A Methodological details of species distribution modelling

(1)

A P P E N D I X A – Methodological details of

s p e c i e s d i s t r i b u t i o n m o d e l l i n g

(2)

(3)

A Methodological details of species distribution modelling

A.1 Method description

General setup and input data

Predictive distribution modelling was applied for the two most prevalent seaduck species in the offshore waters of Polish EEZ, i.e. Long-tailed Duck and Velvet Scoter (Meissner 2015). The input data were collected by the Pomarinus group using ship surveys between June 2012 and June 2013 on OWF BŚ II, Słupsk Bank and other offshore areas in the Polish EEZ and between November 2013 – February 2014 on OWF BŚ II and Słupsk Bank (Figure 1). Each of these areas was surveyed twice per month (with very few exceptions when only one monthly survey was done) during the first survey period (June 2012 – June 2013) and are reported in detail in the Pomarinus report (Meissner 2015). During the second survey period (November 2013 – February 2014) BŚ II area and Słupsk Bank were surveyed once per month except February 2014, when two surveys were conducted. The survey approach generally followed standard ship-based survey methodology for monitoring seabirds at sea (Camphuysen et al. 2004).

Inclusion of data not only collected at BŚ II but also from other offshore areas in the Polish EEZ allows for better characterization of species distribution offshore, as the sample size is increased and more environmental gradients are covered. There is no reason to assume that benthic feeding species have different habitat use preferences on the different survey sites. Survey effort was also balanced and temporally equal between different areas. Therefore, it was assumed that pooling all simultaneously collected data is justifiable and allows achieving more reliable results than analysing different survey areas separately.

Figure 1 The ship survey coverage, representing geographic distribution of the data used for species distribution modelling.

(4)

Due to different spatial coverage during the first and second survey periods, species distribution modelling was done separately for these periods. Also, it included only Long-tailed Duck for the second survey period, as very few Velvet Scoters were registered then when covering only BŚ II area and Słupsk Bank in November 2013 – February 2014. Low registrations of the species did not allow for distribution modelling.

Long-tailed Duck was the most numerous species, and separate models were fitted for each month during the period when the species were present, November 2012 through April 2013;

species observations were too few in October and May for modelling and no long-tailed Ducks were recorded in other months (Meissner 2015). Period of presence of Velvet Scoters is longer than of Long-tailed Ducks and scoter distribution was modelled using data collected between October 2012 and April 2013. Velvet Scoters, however, were less numerous than Long-tailed Ducks, and only one distribution model was fitted for the entire wintering period using month as a factor variable to account for possible differences between months.

Survey data available for modelling the second survey period (Nov 2013 – Feb 2014) was of much smaller volume both spatially and temporally, therefore fitting separate monthly

distribution models for Long-tailed Duck was not possible. Instead, all survey data were used in a single model setup, where survey month was used as a categorical variable.

Based on experiences in modelling waterbird distribution in other parts of the Baltic Sea, seven potentially important predictor (environmental) variables were included in this study; current speed (from DHI hydrodynamic models, DHI 2015), water temperature (from DHI hydrodynamic models, DHI 2015), salinity (from DHI hydrodynamic models, DHI 2015), a filter feeder index (Skov et al. 2011), water depth (Marine Institute, Poland), bottom slope (derived from

bathymetry raster), distance to shipping lanes (Euclidean distance to main shipping lanes in the Baltic Sea) (Figure 2 – Figure 8). AIS data on ship density representing a typical month (August 2010) was received from the Danish Maritime Authority.

Figure 2 Map representing environmental variable “Bathymetry” that was used among potential predictors in Long-tailed Duck and Velvet Scoter distribution models.

(5)

Figure 3 Map representing environmental variable “Bottom slope” that was used among potential predictors in Long-tailed Duck and Velvet Scoter distribution models.

Figure 4 Map representing environmental variable “Distance to shipping lanes” (km) that was used among potential predictors in Long-tailed Duck and Velvet Scoter distribution models.

(6)

Figure 5 Map representing environmental variable “Filter-feeder index” that was used among potential predictors in Long-tailed Duck and Velvet Scoter distribution models.

Figure 6 Map representing environmental variable “Current speed” that was used among potential predictors in Long-tailed Duck and Velvet Scoter distribution models.

(7)

Figure 7 Map representing environmental variable “Salinity” that was used among potential predictors in Long-tailed Duck and Velvet Scoter distribution models.

Figure 8 Map representing environmental variable “Water temperature” that was used among potential predictors in Long-tailed Duck and Velvet Scoter distribution models.

(8)

The filter feeder index represents modelled filter-feeder carrying capacity, which describes the average carrying capacity using an arbitrary scale based on DHI’s hydrodynamic and geo- biochemical model complex BANSAI 3. The carrying capacity is used as a proxy for biomass of mussels and combines a physiology-based growth model for a standard individual with an advection term that replenishes the food ingested by filter-feeders. On a large scale the index depends on the local primary production and on a smaller scale current speed plays an increasing role. Water depth and bottom slope are two parameters of the seabed which both typically have a relatively large influence on the distribution of benthivorous waterbirds. Distance to shipping lanes reflect important human pressure offshore in terms of disturbance from ship traffic, and hence is expected to influence distribution of waterbirds negatively. Hydrodynamic characteristic of water masses may be as a proxy for productive marine areas (current strength) or reflect water mass characteristics that may influence bird (or their prey) physiology and energetics (water temperature and salinity).

The bird observation data were associated with environmental parameters using custom-made data integration tool that has been specifically designed for extracting the data from numerical hydrodynamic model files to data records based on spatial and temporal attributes (coordinates and date-time fields). Values from static environmental data layers (e.g. bathymetry, bottom slope, distance to shipping lanes) were extracted using Interpolate Raster Values to Point tool available in the Marine Geospatial Ecology Tools package for ArcGIS (Roberts et al. 2010).

Bird observations were recalculated to densities and corrected for distance detection bias, which is a standard approach for treating observational data collected from ships (Thomas et al. 2010).

We used correction factors based on effective strip width (ESW) estimated for the same dataset by Meissner (2015): ESW used for Long-tailed Ducks was 169.29 m and for Velvet Scoters 152.76 m. Considering ship transect width of 300 meters, these ESW values yielded correction coefficients of 1.77 and 1.96 respectively. Effective strip width (ESW) defines the distance so that number of animals detected outside of ESW equal the number of animals missed inside the ESW. Correction factors were applied to birds recorded sitting on water inside the transect line and further birds flying in transect recorded during snapshot counts were added to density calculations (see Meissner 2015 for details of ship survey technique).

Finally, bird observation data with associated environmental attributes was aggregated to 500×500 meter grid for separate months by averaging values of bird observations and environmental data if more than one point fell inside the grid cell. The resulting dataset was used in the analyses and included observed bird densities as well as zero values when no birds of a given species was recorded on a transect segment along with the extracted environmental parameters for all data points.

Distribution modelling

A semi-parametric two-step generalized additive model (GAM) was used to account for the zero inflation (a disproportionate large number of zeros) found in the survey data and the potential nonlinear relationships to the environmental variables (Stefansson 1993). As the first step a binomial (presence/absence) model was fitted (Table 1). In the second step all zeros were removed and a gamma model (with a log link) was fitted with the densities of Long-tailed Ducks and Velvet Scoters as the response variables and the environmental variables listed above as the predictor variables (Table 1). In addition to the environmental predictor variables, an interaction term between coordinates X and Y was also included as a predictor in order to account for the variance unexplained by the other variables and to improve the predictive power of the models (Table 1). The predictions from both parts of the models were combined

(multiplied) to yield the final density predictions. The models were fitted using the “mgcv” R package. We defined the maximum degree of smoothing to 5 (k=5) for predictor variables to reduce the risk of over fitting the response curves and no smooth limitation for the interaction term between X and Y coordinates (Table 1).

(9)

Table 1 Main functions in R software code used for modelling seaduck distributions.

require(mgcv) require(MuMIn)

#fitting presence-absence model for the first survey season

PAmodel<-gam(reponse~s(CurrentSpeed,k=5)+s(Temperature,k=5)+s(Salinity, k=5)+s(Distance2ships,k=5) +s(Bathymetry,k=5)+s(BottomSlope,k=5)+s(X,Y),family=binomial,data=data)

#fitting presence-absence model for the second survey season (Long-tailed Duck only) PAmodel<-gam(reponse~s(CurrentSpeed,k=5)+s(Temperature,k=5)+s(Salinity,

k=5)+s(Distance2ships,k=5)+s(Bathymetry,k=5)+s(BottomSlope,k=5)+s(X,Y)+month,family=binomial,data=data)

#fitting positive model

POSmodel<-gam(reponse~s(CurrentSpeed,k=5)+s(Temperature,k=5)+s(Salinity, k=5)+s(Distance2ships,k=5) +s(Bathymetry,k=5)+s(BottomSlope,k=5)+s(X,Y),family=Gamma(link=log),data= data)

#fitting positive model for the second survey season (Long-tailed Duck only)

POSmodel<-gam(reponse~s(CurrentSpeed,k=5)+s(Temperature,k=5)+s(Salinity, k=5)+s(Distance2ships,k=5)+

s(Bathymetry,k=5)+s(BottomSlope,k=5)+s(X,Y)+month,family=Gamma(link=log),data= data)

#automatic model selection based on AICc values

PA_aic<- dredge(PAmodel, beta = FALSE, evaluate = TRUE, rank = "AICc", extra = "adjR^2") POS_aic<- dredge(POSmodel, beta = FALSE, evaluate = TRUE, rank = "AICc", extra = "adjR^2")

#choosing models with AICc <2

topPA_dAIC2<-get.models(PA_aic, subset = delta < 2) topPOS_dAIC2<-get.models(POS_aic, subset = delta < 2)

#averaging best fitting models

avg_top_PA <- model.avg(topPA_dAIC2) avg_top_POS <- model.avg(topPOS_dAIC2)

#Making predictions from each of the models in a set with averaged coefficients

pred_PA<- data.frame(

model = sapply(topPA_dAIC2, predict, newdata = deploy_data, backtransform = TRUE),

averaged.subset = predict(avg_top_PA, newdata = deploy_data, backtransform = TRUE, full = FALSE), averaged.full = predict(avg_top_PA, newdata = deploy_data, backtransform = TRUE, full = TRUE) )

pred_POS<- data.frame(

model = sapply(topPOS_dAIC2, predict, newdata = deploy_data, backtransform = TRUE),

averaged.subset = predict(avg_top_POS, newdata = deploy_data, backtransform = TRUE, full = FALSE), averaged.full = predict(avg_top_POS, newdata = deploy_data, backtransform = TRUE, full = TRUE) )

#Combining predictions from both model parts

Combined_predictions<-pred_PA$averaged.full * pred_POS$averaged.full

Before model fitting we checked the collinearity between the predictor variables, as strong correlation between them can result in inaccurate model parameterization and decreased predictive accuracy (Dormann et al. 2012). High correlation (r = 0.9) was found between depth and filter feeder index (Figure 9), therefore these two variables were not included into the same

(10)

Figure 9 Scatterplot of predictor variables showing correlations between them. The highest correlation of 0.9 was found between filter feeder index (ff) and bathymetry (bath_pl).

All uncorrelated variables were included in the initial global models and further non-significant variables were excluded by employing automatic model selection. Information-theoretical approach was used for model selection, followed by multi-model inference (Burnham and Anderson 2002). AICc – Akaike Information Criterion (AIC) corrected for small sample size was used as the main model selection criterion. Models with ΔAICc (difference in AICc value between the model in question and the best model with the smallest AICc value) less than 2 were retained in the confidence set of models. Their coefficients were also used for model averaging, while their AICc weights (ω) – for the comparison of models and the assessment of the relative importance of predictor variables. Model selection was done using model dredging in

“MuMIn” R package (Barton 2013) and the same procedure was applied for both parts of the model, binomial and gamma as described above (also see Table 1).

Model fit was assessed based on adjusted R² values of separate model parts and standard diagnostic plots used for evaluating GAM performance. We also assessed the model residuals for spatial autocorrelation using a Moran’s I autocorrelogram (defining the nearest

neighbourhood as 1,500 m in the survey data). If strong spatial autocorrelation was found in the model residuals, this would give rise for considering different modelling procedure by refitting the model as a generalized additive mixed model GAMM, which allows for inclusion of a correlation structure and thus is capable for dealing with the spatial autocorrelation.

(11)

Model averaging

To obtain seasonal values of modelled monthly Long-tailed Duck and Velvet Scoter numbers and distribution, model predictions were averaged for months of December, January and February to represent wintering a season. And finally wintering seasons of 2012/2013 and 2013/2014 were averaged together to represent a generic winter season of these species.

Averaging of prediction rasters was done by using a standard ArcGIS tool “raster calculator”.

(12)

A.2 Results: Long-tailed Duck

Due to high over-dispersion in the Long-tailed Duck input data a hurdle modelling approach was chosen, which is a proven method for modelling data with characteristics such as zero-inflation and auto-correlation (Potts & Elith, 2006; Stefánsson, 1996; Heinänen et al., 2008).

Due to high correlation between two predictor variables, ‘bathymetry’ and ‘filter feeder index’, two alternative model sets were tested with these variables included into different sets and keeping all other predictor variables and settings the same. Initial testing made it clear that models with bathymetry performed better based on AIC, R² and other model performance parameters. Therefore, use ‘filter feeder index’ was abandoned and further model selection and predictions were done considering the ‘bathymetry’ variable along with other potential predictors.

As Long-tailed Duck observations were numerous, separate models were fitted for each month during the season of species presence of the first investigation period (June 2012 – June 2013).

The fit of the models representing migration seasons, October and April, however was very poor and therefore these months were excluded from the modelling exercise, as observed birds likely represented not only staging but also migrating individuals, which did not always occur within typical species habitats at sea.

Long-tailed Duck model for November 2012: presence-absence part

All possible presence-absence part models (GAMs) with 6 variable combinations were tested and the best-fitting models with ΔAICc<2 were automatically selected using “MuMIn” package for R (Barton 2013). This procedure allowed selection of the 4 top-ranking models, which differed in predictor variables included, though model performance was relatively similar (Table 2). ‘Bathymetry’, ‘current speed’, ‘distance to ships’ and ‘coordinates XY’ ‘distance to ships’ were consistently occurring in top-ranking models (Table 2). Relative importance of these variables was also high according to AIC weights (Table 3).

Table 2 Structure and fit of the top-ranked binomial models (ΔAICc<2) in the presence/absence part of Long-tailed Duck distribution modelling for November 2012.

Intercept Bathymetry Current speed Distance to ships Salinity Water temperature Coordinates X * Y

R² AICc ΔAICc model weight

-1.14 + + + + + 0.52 301.8 0.00 0.39

-1.14 + + + + + 0.52 302.6 0.81 0.26

-1.08 + + + + 0.50 303.0 1.21 0.21

-1.14 + + + + + + 0.52 303.8 1.97 0.14

(13)

Table 3 Relative variable importance in both parts, presence-absence and positive, of Long-tailed Duck distribution model for November 2012.

Variable Presence-absence model part

Positive model part

Bathymetry 1.00 1.00

Current speed 1.00 0.18

X Y 1.00 ---

Distance to ships 1.00 0.20

Temperature 0.40 0.17

Salinity 0.53 0.12

Response curves of predictor variables in the presence-absence model part were plotted for full model (with all variables included) and represent generic characterisation between the response and predictor variables. The response curves show that probability of bird presence increases with decreasing depth, increasing current speeds, increasing distance from major shipping lanes and with decreasing salinity (Figure 10). Interaction term between coordinates XY suggested that probability of Long-tailed Duck occurrence decreases when moving eastwards along the X axis and peaks at high Y axis values (representing Southern Midsjöbank) considering the input data (Figure 10). The shape representing temperature revealed poor performance of this variable in the presence-absence model part (Figure 10).

Diagnostic plots of the best fitting presence-absence model showed good fit (Figure 11).

Assessment of presence-absence model using Moran’s I test indicated that the model accounted for all auto-correlation that was present in the original input data, which suggests good model performance (Figure 12).

Figure 10 Response curves of the presence/absence GAM model of Long-tailed Duck distribution modelling for November 2012 when all predictor variables were considered.

(14)

Figure 11 Standard diagnostic plots for assessing GAM performance for the best fitting presence- absence part of the Long-tailed Duck distribution model for November 2012.

Figure 12 Moran’s I showing autocorrelation in the original dataset (left) and removed autocorrelation among the residuals of the presence-absence part (right) of the Long-tailed Duck distribution model for November.

Original (input) data Residuals of presence- absence GAM

(15)

Long-tailed Duck model for November 2012: positive part

All possible models (GAMs) of the positive part with 5 variable combinations were tested and the best-fitting models with ΔAICc<2 were automatically selected using “MuMIn” package for R (Barton 2013). This procedure allowed selection of the 5 top-ranking models, which differed in predictor variables included and model performance was also rather variable (Table 4). Only bathymetry consistently occurred in all top-ranking models (Table 4). Relative importance of the variables also indicated that bathymetry had the highest weight among the predictors (Table 3).

Table 4 Structure and fit of the top-ranked binomial models (ΔAICc<2) in the positive part of Long- tailed Duck distribution modelling for November 2012.

Intercept Bathymetry Current speed Distance to ships Salinity Water temperature

R² AICc ΔAICc ΔAIC weight

4.55 + 0.41 1178.0 0.00 0.328

4.53 + + 0.42 1179.1 1.02 0.197

4.52 + + 0.43 1179.2 1.16 0.183

4.51 + + 0.44 1179.3 1.29 0.172

4.55 + + 0.41 1180.0 1.99 0.121

Response curves of predictor variables in the positive model part were plotted for the full model (with all variables included) and represent generic characterisation between the response and predictor variables. The response curves showed that bird densities increase with increasing current speed, declining salinity, increasing distance to shipping major lanes and at shallower depths (Figure 13). The shape representing ‘temperature’ indicated poor performance of this variable in the positive model part (Figure 13).

Diagnostic plots of the best fitting positive model showed moderate fit (Figure 14). Assessment of the positive model using Moran’s I test indicated that there was no auto-correlation present in the original input data (Figure 15), therefore no testing for autocorrelation in the model residuals was needed.

(16)

Figure 13 Response curves of the positive GAM model of Long-tailed Duck distribution modelling for November 2012 when all predictor variables were considered.

Figure 14 Standard diagnostic plots for assessing GAM performance for the best fitting positive part of the Long-tailed Duck distribution model for November 2012.

(17)

Figure 15 Moran’s I showing absence of autocorrelation in the original input dataset for the positive part of the Long-tailed Duck distribution model for November 2012.

The highest ranked presence-absence and positive models were averaged and used for predicting Long-tailed Duck distribution for November in the entire area that was considered.

The results show that high densities of this species occurred on the Słupsk Bank and Southern Midsjöbank, and moderate to low densities in the surrounding areas (Figure 16). The predicted Long-tailed Duck densities in the BŚ II area were low to high ranging from less than 1 bird/km² to 80 birds/km²(Figure 16). The highest densities of this species with the BŚ II area occurred in the south-western part which is closest to the Słupsk Bank.

Figure 16 Modelled Long-tailed Duck distribution in Polish offshore waters for November 2012.

Original (input) data

(18)

Long-tailed Duck model for December 2012: presence-absence part

All possible presence-absence part models (GAMs) with 6 variable combinations were tested and the best-fitting models with ΔAICc<2 were automatically selected using “MuMIn” package for R (Barton 2013). This procedure allowed selection of the 7 top-ranking models, which differed in predictor variables included, though model performance was similar (Table 5).

‘current speed’ and ‘coordinates XY’ appeared in all top models and ‘salinity was also frequently occurring (Table 5). Relative importance of these variables was also high according to AIC weights (Table 6).

Table 5 Structure and fit of the top-ranked binomial models (ΔAICc<2) in the presence/absence part of Long-tailed Duck distribution modelling for December 2012.

-0.99 + + + 0.64 237.1 0.00 0.25

-0.99 + + + + 0.64 237.6 0.54 0.19

-0.99 + + + + 0.65 238.4 1.31 0.13

-0.93 + + 0.64 238.5 1.37 0.13

-0.99 + + + + + 0.65 238.9 1.80 0.10

-0.99 + + + + 0.64 238.9 1.80 0.10

-0.93 + + + 0.64 238.9 1.80 0.10

Table 6 Relative variable importance in both parts, presence-absence and positive, of Long-tailed Duck distribution model for December 2012.

Variable Presence-absence model

part

Positive model part

X Y 1.00 1.00

Salinity 0.77 0.59

Response curves of predictor variables in the presence-absence model part were plotted for full model (with all variables included) and represent generic characterisation between the response and predictor variables. The response curves show that probability of bird presence increases with increasing current speeds, decreasing salinity and shallower depths (Figure 17).

Interactions with water temperature and distance to ships show weak importance of these variables on species occurrence in December 2012 (Figure 17). Interaction term between

(19)

coordinates XY suggested that probability of Long-tailed Duck occurrence was highest in the eastern half along the X axis considering the input data (Figure 17).

Figure 17 Response curves of the presence/absence GAM model of Long-tailed Duck distribution modelling for December 2012 when all predictor variables were considered.

(20)

Figure 18 Standard diagnostic plots for assessing GAM performance for the best fitting presence- absence part of the Long-tailed Duck distribution model for December 2012.

Figure 19 Moran’s I showing autocorrelation in the original dataset (left) and removed autocorrelation among the residuals of the presence-absence part (right) of the Long-tailed Duck distribution model for December 2012.

Long-tailed Duck model for December 2012: positive part

All possible models (GAMs) of the positive part with 6 variable combinations were tested and the best-fitting models with ΔAICc<2 were automatically selected using “MuMIn” package for R (Barton 2013). This procedure allowed selection of the 8 top-ranking models, which differed in

(21)

predictor variables included but model performance was similar (Table 7). Bathymetry and coordinates XY consistently occurred in all top-ranking models, but current speed and salinity were also common (Table 7). Relative importance of the variables also indicated that

bathymetry and coordinates XY had the highest weight among the predictors, followed by current speed and salinity (Table 6).

Table 7 Structure and fit of the top-ranked binomial models (ΔAICc<2) in the positive part of Long- tailed Duck distribution modelling for December 2012.

Intercept Bathymetry Current speed Distance to ships Salinity Water temperature XY R² AICc ΔAICc ΔAIC weight

4.73 + + + + 0.49 1198.0 0.00 0.23

4.73 + + + + + 0.50 1199.2 1.14 0.13

4.77 + + + 0.48 1199.2 1.16 0.13

4.83 + + 0.43 1199.2 1.23 0.13

4.78 + + + + 0.47 1199.6 1.58 0.11

4.74 + + + + + 0.49 1199.7 1.72 0.10

4.81 + + + 0.44 1199.8 1.79 0.10

4.81 + + + 0.44 1200.0 1.97 0.09

Response curves of predictor variables in the positive model part were plotted for the full model (with all variables included) and represent generic characterisation between the response and predictor variables. The response curves showed that bird densities were higher at lower current speed, higher water temperature, lower salinity, and at shallower depths (Figure 20).

Diagnostic plots of the best fitting positive model showed reasonable model fit (Figure 21).

Assessment of the positive model using Moran’s I test indicated that there was no autocorrelation present in the original input data of positive observations (Figure 22), therefore no testing for autocorrelation in the model residuals was needed.

(22)

Figure 20 Response curves of the positive GAM model of Long-tailed Duck distribution modelling for December 2012 when all predictor variables were considered.

Figure 21 Standard diagnostic plots for assessing GAM performance for the best fitting positive part of the Long-tailed Duck distribution model for December 2012.

(23)

Figure 22 Moran’s I showing absence of autocorrelation in the original input dataset for the positive part of the Long-tailed Duck distribution model for December 2012.

The highest ranked presence-absence and positive models were averaged and used for predicting Long-tailed Duck distribution for December 2012. The predicted Long-tailed Duck densities in the BŚ II area were very high, up to 420 birds/km²(Figure 23).

Figure 23 Modelled Long-tailed Duck distribution in BŚ II area for December 2012.

(24)

Long-tailed Duck model for January 2013: presence-absence part

Table 8 Structure and fit of the top-ranked binomial models (ΔAICc<2) in the presence/absence part of Long-tailed Duck distribution modelling for January 2013.

-0.49 + + + + 0.31 318.1 0.00 0.22

-0.54 + + + + + 0.34 318.2 0.07 0.22

-0.53 + + + 0.35 318.4 0.34 0.19

-0.53 + + + + 0.35 318.5 0.38 0.19

-0.53 + + + + 0.34 318.5 0.38 0.19

Table 9 Relative variable importance in presence-absence and positive parts of best (<ΔAICc) Long- tailed Duck distribution models for January 2013.

Variable Presence-absence model part

Positive model part

X Y 0.78 1.00

Temperature 1.00 ---

Salinity 0.63 0.27

The response curves show that probability of bird presence increases with decreasing water temperature, increasing distance to shipping lanes and at shallower depths (Figure 24). Current speed and salinity showed weak effect on species occurrence in January 2013 (Figure 24).

Interaction term between coordinates XY suggested that probability of Long-tailed Duck occurrence was highest in the eastern half along the X axis and northern part along Y axis considering the input data (Figure 24).

(25)

Figure 24 Response curves of the presence/absence GAM model of Long-tailed Duck distribution modelling for January 2013 when all predictor variables were considered.

Figure 25 Standard diagnostic plots for assessing GAM performance for the best fitting presence- absence part of the Long-tailed Duck distribution model for January 2013.

(26)

Figure 26 Moran’s I showing autocorrelation in the original dataset (left) and removed autocorrelation among the residuals of the presence-absence part (right) of the Long-tailed Duck distribution model for January 2013.

Long-tailed Duck model for January 2013: positive part

All possible models (GAMs) of the positive part with 6 variable combinations were tested and the best-fitting models with ΔAICc<2 were automatically selected using “MuMIn” package for R (Barton 2013). This procedure allowed selection of the 13 top-ranking models, which differed in predictor variables included but model performance was similar (Table 10). Coordinates XY consistently occurred in all top-ranking models, and bathymetry and current speed were also common (Table 10). Relative importance of the variables also indicated that coordinates XY had the highest weight among the predictors, but weight of other variables were relatively low (Table 9).

(27)

Table 10 Structure and fit of the top-ranked binomial models (ΔAICc<2) in the positive part of Long- tailed Duck distribution modelling for January 2013.

4.59 + 0.35 1480.3 0.00 0.13

4.64 + + 0.30 1480.4 0.11 0.12

4.59 + + 0.35 1480.8 0.44 0.10

4.65 + + + 0.31 1451.0 0.68 0.09

4.57 + + 0.37 1481.1 0.79 0.09

4.56 + + + 0.37 1481.5 1.20 0.07

4.67 + + + 0.28 1481.6 1.25 0.07

4.61 + + + 0.33 1481.6 1.25 0.07

4.64 + + + + 0.31 1482.1 1.81 0.05

4.61 + + + + 0.34 1482.1 1.82 0.05

4.23 + + 0.32 1482.2 1.87 0.05

4.66 + + + + 0.29 1482.2 1.90 0.05

4.61 + + + 0.33 1482.3 1.97 0.05

The response curves showed that bird densities were increasing with increasing current speed, lower salinity, and at shallower depths (Figure 27). Temperature and distance to ships were not important variables defining bird densities in January 2013 (Figure 27).

Diagnostic plots of the best fitting positive model showed moderate model fit (Figure 28).

(28)

Figure 27 Response curves of the positive GAM model of Long-tailed Duck distribution modelling for January when all predictor variables were considered.

Figure 28 Standard diagnostic plots for assessing GAM performance for the best fitting positive part of the Long-tailed Duck distribution model for January 2013.

(29)

Figure 29 Moran’s I showing absence of autocorrelation in the original input dataset for the positive part of the Long-tailed Duck distribution model for January 2013.

The highest ranked presence-absence and positive models were averaged and used for predicting Long-tailed Duck distribution for January 2013. The predicted Long-tailed Duck densities in the BŚ II area were very high, 6-180 birds/km²(Figure 30).

Figure 30 Modelled Long-tailed Duck distribution in BŚ II area for January 2013.

(30)

Long-tailed Duck model for February 2013: presence-absence part

Table 11 Structure and fit of the top-ranked binomial models (ΔAICc<2) in the presence/absence part of Long-tailed Duck distribution modelling for February 2013.

-0.53 + + + + 0.25 430.1 0.00 0.38

-0.55 + + + 0.25 430.8 0.69 0.27

-0.54 + + + 0.22 431.4 1.28 0.20

-0.53 + + + + 0.24 431.9 1.77 0.16

Table 12 Relative variable importance in presence-absence and positive parts of best (<ΔAICc) Long- tailed Duck distribution models for February 2013.

Variable Presence-absence model part Positive model part

Current speed 1.00 ---

X Y 0.64 1.00

Temperature 0.16 ---

Salinity --- ---

The response curves show that probability of bird presence increases with decreasing current speed and at shallower depths (Figure 31). Current speed and temperature showed weak effect on species occurrence in February 2013 (Figure 31). Bird probability of occurrence in relation to distance to shipping lanes is rather difficult to interpret as it shows decreasing tendency with increasing distance to shipping lanes (Figure 31).

Diagnostic plots of the best fitting presence-absence model showed moderate fit (Figure 32).

(31)

Figure 31 Response curves of the presence/absence GAM model of Long-tailed Duck distribution modelling for February 2013 when all predictor variables were considered.

Figure 32 Standard diagnostic plots for assessing GAM performance for the best fitting presence- absence part of the Long-tailed Duck distribution model for February 2013.

(32)

Figure 33 Moran’s I showing autocorrelation in the original dataset (left) and removed autocorrelation among the residuals of the presence-absence part (right) of the Long-tailed Duck distribution model for February 2013.

Long-tailed Duck model for February 2013: positive part

All possible models (GAMs) of the positive part with 6 variable combinations were tested and the best-fitting models with ΔAICc<2 were automatically selected using “MuMIn” package for R (Barton 2013). This procedure allowed selection of the 2 top-ranking models, which were similar in predictor variables and performance (Table 13). Coordinates XY and bathymetry occurred in both top-ranking models, and distance to shipping lanes occurred in one (Table 13). Relative importance of the variables also indicated that coordinates XY and bathymetry had the highest weight (Table 12).

Table 13 Structure and fit of the top-ranked binomial models (ΔAICc<2) in the positive part of Long- tailed Duck distribution modelling for February.

3.64 + + 0.52 1586.8 0.00 0.58

3.64 + + + 0.52 1587.4 0.62 0.42

The response curves showed that bird densities were increasing with decreasing depth and highest densities occurred in Słupsk Bank area according to the interaction term XY (Figure 34).

Other variables were of little importance in determining bird densities in February 2013 (Figure 34).

Diagnostic plots of the best fitting positive model showed reasonably good model fit (Figure 35).

(33)

Figure 34 Response curves of the positive GAM model of Long-tailed Duck distribution modelling for February 2013 when all predictor variables were considered.

Figure 35 Standard diagnostic plots for assessing GAM performance for the best fitting positive part of the Long-tailed Duck distribution model for February 2013.

(34)

Figure 36 Moran’s I showing absence of autocorrelation in the original input dataset for the positive part of the Long-tailed Duck distribution model for February 2013.

The highest ranked presence-absence and positive models were averaged and used for predicting Long-tailed Duck distribution for February 2013. The predicted Long-tailed Duck densities in the BŚ II area were low to high, 0.2-115 birds/km²(Figure 37).

Figure 37 Modelled Long-tailed Duck distribution in the BŚ II area for February 2013.

(35)

Long-tailed Duck model for March 2013: presence-absence part

All possible presence-absence part models (GAMs) with 6 variable combinations were tested and the best-fitting models with ΔAICc<2 were automatically selected using “MuMIn” package for R (Barton 2013). This procedure allowed selection of the 4 top-ranking models, which

differed in predictor variables included, though model performance was rather similar (Table 14).

Bathymetry, water temperature and ‘coordinates XY’ appeared in all top models (Table 14).

Relative importance of these variables was also the highest according to AIC weights (Table 15).

Table 14 Structure and fit of the top-ranked binomial models (ΔAICc<2) in the presence/absence part of Long-tailed Duck distribution modelling for March 2013.

-1.68 + + + + 0.67 271.3 0.00 0.34

-0.21 + + + 0.57 271.3 0.02 0.34

-0.21 + + + + 0.57 272.6 1.32 0.18

-0.17 + + + + + 0.57 273.1 1.81 0.14

Table 15 Relative variable importance in presence-absence and positive parts of best (<ΔAICc) Long- tailed Duck distribution models for March 2013.

X Y 1.00 1.00

Salinity 0.14 ---

The response curves show that probability of bird presence increases with decreasing water depths, closer to shipping lanes and peaks at water temperature of 2.5°C (Figure 38). Current speed and salinity showed weak effect on species occurrence in March 2013 (Figure 38).

Diagnostic plots of the best fitting presence-absence model showed good model fit (Figure 39).

(36)

Figure 38 Response curves of the presence/absence GAM model of Long-tailed Duck distribution modelling for March 2013 when all predictor variables were considered.

Figure 39 Standard diagnostic plots for assessing GAM performance for the best fitting presence- absence part of the Long-tailed Duck distribution model for March 2013.

(37)

Figure 40 Moran’s I showing autocorrelation in the original dataset (left) and removed autocorrelation among the residuals of the presence-absence part (right) of the Long-tailed Duck distribution model for March 2013.

Long-tailed Duck model for March 2013: positive part

All possible models (GAMs) of the positive part with 6 variable combinations were tested and the best-fitting models with ΔAICc<2 were automatically selected using “MuMIn” package for R (Barton 2013). This procedure allowed selection of 8 top-ranking models, which varied in predictor variables but were similar in general model performance (Table 16). Current speed and coordinates XY occurred in all top-ranking models (Table 16). Relative importance of the variables also indicated that coordinates XY and current speed had the highest weight (Table 15).

(38)

Table 16 Structure and fit of the top-ranked binomial models (ΔAICc<2) in the positive part of Long- tailed Duck distribution modelling for March 2013.

Intercept Bathymetry Current speed Distance to ships Water temperature XY R² AICc ΔAICc ΔAIC weight

3.88 + + + 0.57 1404.8 0.00 0.21

3.91 + + 0.55 1405.6 0.81 0.14

3.88 + + + + 0.57 1405.6 0.81 0.14

3.88 + + + + 0.5 1405.8 0.96 0.13

3.91 + + + 0.55 1406.4 1.53 0.10

3.91 + + + 0.55 1406.4 1.59 0.10

3.88 + + + + + 0.57 1406.6 1.73 0.09

3.95 + + 0.51 1406.7 1.89 0.08

The response curves showed that bird densities were increasing with decreasing depth and temperature and the highest densities occurred in Slupsk Bank area according to the interaction term XY (Figure 41). Other variables were of little importance in determining bird densities in March 2013 (Figure 41).

Assessment of the positive model using Moran’s I test indicated that slight auto-correlation that was present in the first lag of the original input data was successfully accounted for by the model so that no autocorrelation was left in the model residuals (Figure 43).

(39)

Figure 41 Response curves of the positive GAM model of Long-tailed Duck distribution modelling for March 2013 when all predictor variables were considered.

Figure 42 Standard diagnostic plots for assessing GAM performance for the best fitting positive part of the Long-tailed Duck distribution model for March 2013.

(40)

Figure 43 Moran’s I showing absence of autocorrelation in the original input dataset for the positive part of the Long-tailed Duck distribution model for March 2013.

The highest ranked presence-absence and positive models were averaged and used for predicting Long-tailed Duck distribution for March 2013. The predicted Long-tailed Duck densities in the BŚ II area were high, 14-110 birds/km²(Figure 44).

Figure 44 Modelled Long-tailed Duck distribution in the BŚ II area for March 2013.

(41)

Long-tailed Duck model for November 2013 – February 2014: presence-absence part

All possible presence-absence part models (GAMs) with 8 variable combinations were tested and the best-fitting models with ΔAICc<2 were automatically selected using “MuMIn” package for R (Barton 2013). This procedure allowed selection of the 9 top-ranking models, which

differed in predictor variables included, though model performance was rather similar (Table 17).

Bathymetry and month in all top models and salinity, bottom slope and water temperature were included in many of them (Table 17). Relative importance of these variables was also the highest according to AIC weights (Table 18).

Table 17. Structure and fit of the top-ranked binomial models (ΔAICc<2) in the presence/absence part of Long-tailed Duck distribution modelling for November 2013 – February 2014.

Intercept Bathymetry Bottom slope Current speed Distance to ships Salinity Water temperature Coordinates X * Y Month

-0.99 + ⁺ + ^0.38 ^538.5 ^0.00 ^0.162

-1.78 + ⁺ ⁺ + ^0.39 ^538.6 ^0.09 ^0.154

-0.99 + ⁺ ⁺ + ^0.38 ^538.8 ^0.28 ^0.141

-1.68 + ⁺ ⁺ ⁺ + ^0.39 ^539.0 ^0.44 ^0.130

-2.22 + ⁺ + ^0.38 ^539.1 ^0.53 ^0.124

-2.23 + ⁺ ⁺ + ^0.38 ^539.8 ^1.26 ^0.086

-0.97 + ⁺ ⁺ + ^0.38 ^540.1 ^1.53 ^0.075

-0.98 + ⁺ ⁺ ⁺ + ^0.38 ^540.3 ^1.78 ^0.067

-0.99 + ⁺ + ^0.38 ^540.5 ^1.98 ^0.060

Table 18 Relative variable importance in presence-absence and positive parts of best (<ΔAICc) Long- tailed Duck distribution models for November 2013 – February 2014.

Bathymetry 1.0 1.00

Bottom slope 0.42 --

Distance to ships 0.06 --

Salinity 0.79 0.21

X Y 0 1.00

Month 1.0 0

The response curves show that probability of bird presence increases with decreasing water depths, in more saline water, at low water temperature, and increasing bottom slope (Figure 45).

Other variables showed weak effect on species occurrence (Figure 45).

Diagnostic plots of the best fitting presence-absence model showed reasonable model fit (Figure 46). Assessment of presence-absence model using Moran’s I test indicated that the model

(42)

Figure 45 Response curves of the presence/absence GAM model of Long-tailed Duck distribution modelling for November 2013 – February 2014 when all predictor variables were considered.

(43)

Figure 46 Standard diagnostic plots for assessing GAM performance for the best fitting presence- absence part of the Long-tailed Duck distribution model for November 2013 – February 2014.

Figure 47 Moran’s I showing autocorrelation in the original dataset (left) and removed autocorrelation among the residuals of the presence-absence part (right) of the Long-tailed Duck distribution model for November 2013 – February 2014.

Long-tailed Duck model for November 2013 – February 2014: positive part

All possible models (GAMs) of the positive part with 6 variable combinations were tested and the best-fitting models with ΔAICc<2 were automatically selected using “MuMIn” package for R (Barton 2013). This procedure allowed selection of 3 top-ranking models, which varied in

(44)

current speed and coordinates XY occurred in all top-ranking models (Table 19). Relative importance of the variables also indicated that these had the highest weight (Table 18).

Table 19 Structure and fit of the top-ranked binomial models (ΔAICc<2) in the positive part of Long- tailed Duck distribution modelling for November 2013 – February 2014.

Intercept Bathymetry Current speed Salinity Water temperature XY Month

4.03 + + + + 0.39 2059.3 0.00 0.43

4.03 + + + 0.38 2059.7 0.40 0.35

4.03 + + + + 0.39 2060.7 1.41 0.21

The response curves showed that bird densities were increasing with decreasing depth and and the highest densities occurred on Slupsk Bank and southern part of BS II area according to the interaction term XY (Figure 48). Other variables were of little importance in determining bird densities in March 2013 (Figure 48).

Assessment of the positive model using Moran’s I test indicated that slight auto-correlation that was present in the first lag of the original input data was successfully accounted for by the model so that no autocorrelation was left in the model residuals (Figure 50).

Figure 48 Response curves of the positive GAM model of Long-tailed Duck distribution modelling for November 2013 – March 2014 when all predictor variables were considered.

(45)

Figure 49 Standard diagnostic plots for assessing GAM performance for the best fitting positive part of the Long-tailed Duck distribution model for November 2013 – February 2014.

Figure 50 Moran’s I showing absence of autocorrelation in the original input dataset for the positive part of the Long-tailed Duck distribution model for November 2013 – February 2014.

The highest ranked presence-absence and positive models were averaged and used for predicting Long-tailed Duck distribution for separate months of the second survey season (winter 2013-2014).

In contrast to the first survey season, Long-tailed Duck densities in the BŚ II area were quite a bit lower in winter 2013-2014, and ranged from 0.01-4 birds/km²in November to 0.3-

(46)

Figure 51 Modelled Long-tailed Duck distribution in the BŚ II area for November 2013.

Figure 52 Modelled Long-tailed Duck distribution in the BŚ II area for December 2013.

(47)

Figure 53 Modelled Long-tailed Duck distribution in the BŚ II area for January 2013.

Figure 54 Modelled Long-tailed Duck distribution in the BŚ II area for February 2013.

(48)

Seasonal averaging of wintering Long-tailed Duck numbers and distribution

Distribution modelling revealed that abundance of Long-tailed Ducks on the BS II area is highly variable within different months of the same wintering season and even more so between different winters. Aiming to rule out that chance events or unusual season influence the

assessment too much, the average values for winter season were calculated separately for both winter seasons when the surveys have been conducted and overall winter average. The results revealed that Long-tailed Duck densities were considerably higher during the winter 2012/2013 (Figure 55) compared to the following winter 2013/2014 (Figure 56), the distribution patterns being rather similar. The averaged distribution for two winters, obviously represented mean values recorded during the two winter of observations (Figure 57).

Figure 55 Averaged Long-tailed Duck distribution in the BŚ II area for winter 2012/2013 (Dec-Feb).

(49)

Figure 56 Averaged Long-tailed Duck distribution in the BŚ II area for winter 2013/2014 (Dec-Feb).

Figure 57 Averaged Long-tailed Duck distribution in the BŚ II area for winter conditions (Dec-Feb).

(50)

A.3 Results: Velvet Scoter

Due to high over-dispersion in the input data a hurdle modelling approach was chosen, which is a suitable method to model the data with characteristics such as zero-inflation and auto-

correlation (Potts & Elith, 2006; Stefánsson, 1996; Heinänen et al., 2008).

Due to high correlation between two predictor variables, ‘bathymetry’ and ‘filter feeder index’ two alternative model sets were tested with these variables included into different sets and keeping all other predictor variables and settings the same. Initial testing made it clear that models with bathymetry performed better based on AIC, R² and other model performance parameters.

Therefore, use ‘filter feeder index’ was abandoned and further model selection and predictions were done considering ‘bathymetry’.

Presence-absence part of the Velvet Scoter distribution model

All possible presence-absence part models (GAMs) with 7 variable combinations were tested and the best-fitting models with ΔAICc<2 were automatically selected using “MuMIn” package for R (Barton 2013). This procedure allowed selection of the 5 top-ranking models, which differed in predictor variables included, though model performance was relatively similar (Table 20). ‘Bathymetry’, ‘current speed’, ‘coordinates XY’ and ‘distance to ships’ were consistently occurring in top-ranking models (Table 20). Relative importance of these variables was also high according to AIC weights (Table 21).

Table 20 Structure and fit of the top-ranked binomial models (ΔAICc<2) in the presence/absence part of Velvet Scoter distribution modelling.

Intercept Bathymetry Current speed Distance to ships Salinity Bottom slope Water temperature Coordinates X * Y

-4.71 + + + + 0.34 531.1 0.00 0.341

-4.43 + + + 0.34 531.9 0.91 0.228

-4.73 + + + + + 0.34 532.7 1.61 0.152

-4.70 + + + + + 0.34 532.9 1.76 0.142

-4.71 + + + + + 0.34 632.9 1.81 0.138

(51)

Table 21 Relative variable importance in both parts, presence-absence and positive, of Velvet Scoter distribution models.

Current speed 1.00

X Y 1.00 0.62

Bottom slope 0.15 0.71

Salinity 0.14 0.61

Response curves of predictor variables in the presence-absence model part were plotted for full model (with all variables included) and represent generic characterisation between the response and predictor variables. The response curves show that probability of bird presence increases with decreasing depth, at moderate current speeds and increasing distance from major shipping lanes (Figure 58). Interaction term between coordinates XY suggested that probability of Velvet Scoter occurrence increases when moving eastwards along the X axis considering the input data (Figure 58). Shapes representing ‘bottom slope’, ‘temperature’ and ‘salinity’ indicated rather poor performance of these variables in the presence-absence model part (Figure 58).

Diagnostic plots of the best fitting presence-absence model showed moderate fit (Figure 59).

Assessment of the positive model using Moran’s I test indicated that there was no auto-

correlation present in the original input data (Figure 60), therefore no testing for autocorrelation in the model residuals was needed.

(52)

Figure 58 Response curves of the presence/absence GAM model of Velvet Scoter distribution modelling when all predictor variables were considered.

(53)

Figure 59 Standard diagnostic plots for assessing GAM performance for the best fitting presence- absence part of the Velvet Scoter distribution model.

Figure 60 Moran’s I showing autocorrelation in the original dataset (left) and removed autocorrelation among the residuals of the presence-absence part (right) of the Velvet Scoter distribution model.

Positive part of the Velvet Scoter distribution model

All possible models (GAMs) of the positive part with 7 variable combinations were tested and the best-fitting models with ΔAICc<2 were automatically selected using “MuMIn” package for R (Barton 2013). This procedure allowed selection of the 10 top-ranking models, which differed in predictor variables included and model performance was also rather variable (Table 22). All

(54)

variables except current speed occurred in top-ranking models (Table 22). Relative importance of these variables was quite similar (Table 21).

Table 22 Structure and fit of the top-ranked binomial models (ΔAICc<2) in the positive part of Velvet Scoter distribution modelling.

Intercept Bathymetry Current speed Distance to ships Salinity Bottom slope Water temperature Coordinates X * Y

2.55 + + + + 0.40 633.5 0.00 0.159

2.56 + + + + 0.37 633.6 0.11 0.151

2.56 + + + + 0.38 634.0 0.50 0.124

2.55 + + + + + 0.40 634.2 0.71 0.111

2.54 + + + + + + 0.41 634.5 1.02 0.095

2.52 + + + 0.45 634.6 1.19 0.088

2.57 + + + 0.37 634.9 1.47 0.076

2.57 + + + 0.37 635.0 1.54 0.074

2.52 + + + + 0.44 635.3 1.89 0.062

2.53 + + + + + 0.43 635.4 1.94 0.060

Response curves of predictor variables in the positive model part were plotted for full model (with all variables included) and represent generic characterisation between the response and predictor variables. The response curves showed that bird densities increase with decreasing temperature, salinity and depth, and increasing distance to shipping major lanes (Figure 61).

Interaction with bottom slope show higher densities at flat bottoms and also at bottom slopes between 0.1 and 0.17 degrees (Figure 61). Interaction term between coordinates XY suggested that densities of Velvet Scoters increase when moving eastwards along the X axis considering the input data (Figure 61). The shape representing ‘current speed’ indicated poor performance of this variable in the positive model part (Figure 61).

There was no autocorrelation found in the positive part of the original input data according to Moran’s I test therefore there was no need to check for autocorrelation in the model residuals (Figure 63).

(55)

Figure 61 Response curves of the positive GAM model of Velvet Scoter distribution modelling when all predictor variables were considered.

(56)

Figure 62 Standard diagnostic plots for assessing GAM performance for the best fitting positive part of the Velvet Scoter distribution model.

Figure 63 Moran’s I showing absence of autocorrelation in the original input dataset of the positive part of the Velvet Scoter distribution model.

The highest ranked presence-absence and positive models were averaged and used for predicting Velvet Scoter distribution in the entire area which was considered. The results showed that low to moderate densities of this species occurred in the southern part of the study area, Słupsk Bank and towards the Polish mainland coast, and near the Southern Midsjöbank (Figure 64). Predicted Velvet Scoter densities in the BŚ II area were very low and never exceeded 0.5 birds/km² (Figure 65 – Figure 70).

(57)

Figure 64 Modelled Velvet Scoter distribution in Polish offshore waters for November 2012.

Figure 65 Modelled Velvet Scoter distribution in BŚ II area for October 2012.

(58)

Figure 66 Modelled Velvet Scoter distribution in BŚ II area for December 2012.

Figure 67 Modelled Velvet Scoter distribution in BŚ II area for January 2013.

(59)

Figure 68 Modelled Velvet Scoter distribution in BŚ II area for February 2013.

Figure 69 Modelled Velvet Scoter distribution in BŚ II area for March 2013.

(60)

Figure 70 Modelled Velvet Scoter distribution in BŚ II area for April 2013.

Seasonal averaging of wintering Velvet Scoter numbers and distribution

Distribution modelling revealed that abundance of Velvet Scoter abundance is low on the BS II area not very variable among different months. In order to follow the same data processing approach as for the Long-tailed Duck, the average values for winter season were calculated for Velvet Scoters as well (Figure 71).

(61)

Figure 71 Averaged Velvet Scoter distribution in the BŚ II area for winter 2012/2013 (Dec-Feb).

(62)

A.4 References

Barton, K. (2013) MuMIn: Multi-model inference. R package version 1.9.0.

Burnham, K. P. and Anderson, D. R (2002) Model selection and multimodel inference: a practical information-theoretic approach. 2nd ed. New York, Springer-Verlag.

Camphuysen, K.C.J, Fox, A.D., Leopold, M.F. and Petersen, I.K. (2004) Towards standardised seabirds at sea census techniques in connection with environmental impact assessments for offshore wind farms in the U.K. Report COWRIE – BAM- 02-2002.

DHI (2015). Environmental Impact Assessment of Baltyk Środkowy II Offshore Wind Farm Model setup and hydrographic impact assessment.

Dormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré, G., García Marquéz, J. R., Gruber, B., Lafourcade, B., Leitão, P. J., Münkemüller, T., McClean, C., Osborne, P. E., Reineking, B., Schröder, B., Skidmore, A. K., Zurell, D. & Lautenbach, S. (2013) Collinearity: a review of methods to deal with it and a simulation study evaluating their performance.

Ecography, 36, 27-46.

Durinck, J., Skov, H., Jensen, F.P. and Pihl, S. 1994. Important marine areas for wintering birds in the Baltic Sea. EU DG XI research contract no. 2242/90-09-01. Ornis Consult report.

Heinänen, S., Rönkä, M. & Numers, von, M. (2008) Modelling the occurrence and abundance of a colonial species, the arctic tern Sterna paradisaea in the archipelago of SW Finland.

Ecography, 31,601-611.

Meissner, W. (2015). Ornithological monitoring of the area of the planned offshore wind farm

„Bałtyk Środkowy II” Final report and the results of the monitoring. Gdańsk.

Potts, J. M. & Elith, J. (2006) Comparing species abundance models. Ecological Modelling, 199, 153-163.

Roberts, J.J., Best, B.D., Dunn, D.C., Treml, E.A., and Halpin, P.N. (2010) Marine Geospatial Ecology Tools: An integrated framework for ecological geoprocessing with ArcGIS, Python, R, MATLAB, and C++. Environmental Modelling & Software 25: 1197-1207.

Skov, H., Heinänen, S., Žydelis, R., Bellebaum, J., Bzoma, S., Dagys, M., Durinck, J., Garthe, S., Grishanov, G., Hario, M., Kieckbusch, J.J., Kube, J., Kuresoo, A., Larsson, K., Luigujoe, L., Meissner, W., Nehls, H.W., Nilsson, L., Petersen, I.K., Roos, M.M., Pihl, S., Sonntag, N., Stock, A., Stipniece, A., Wahl, J., (2011). Waterbird populations and pressures in the Baltic Sea.

TemaNord 2011:550, Nordic Council of Ministers, Copenhagen, Denmark.

Stefánsson, G. (1996) Analysis of groundfish survey abundance data: combining the GLM and delta approaches. ICES Journal of Marine Science, 53, 577-588.

Thomas, L., Buckland, S.T., Rexstad, E.A., Laake, J.L., Strindberg, S., Hedley, S.L., Bishop, J.R.B., Marques, T.A. and Burnham, K.P. 2010. Distance Software: design and analysis of distance sampling surveys for estimating population size. Journal of Applied Ecology, 47, 5-14.

(63)