Impact of climate change on seabird species off the east coast of Scotland and potential implications for environmental assessments: study

Appendix C. Survival modelling

Models

The structure of the model for adult survival was similar to that of the binomial model for productivity; in this case the binomial numerator is equal to the number of breeding adults in year that were already breeding adults in year , and the denominator is equal to the number of breeding adults in year . The transformed annual adult survival probability is assumed to be a linear function of environmental variables (relating to climate, but, depending on the biology of the species, not necessarily the same variables that are assumed to impact productivity), and to also depend on random effects for "site", "year" and the interaction of "site" and "year". We initially considered a logit transformation (as the survival probability must lie between zero and one), but the resulting model encountered irresolvable estimation problems due to incompatibility between the model and data, and we therefore used a log transformation instead. The use of a log transformation means the survival probabilities are assumed to be positive, but not necessarily below one – although it is logically impossible to have survival probabilities above one, using a model that allows this (e.g. a log rather than logit transformation) was possible whereas a model that prevents this was not. This is probably because of inconsistencies between the productivity data and count data, which can only be resolved by allowing the survival probability to be greater than one in some instances.

The number of breeding adults in year that were already breeding adults in year is not something that can usually be observed directly, but it will be equal to:

Breeding adults in year - Number of new recruits in year

The former quantity is observable; under a standard Leslie matrix model (with three age classes) the latter quantity can be assumed to follow a Binomial distribution where the denominator is the number of chicks born in year , where denotes the age at first breeding, and the probability of survival is equal to .

We considered two possible models for juvenile survival:

Model A: juvenile survival is a single unknown value, constant across time and space, whose value is derived from the species-level estimate in Horswill & Robinson (2015), so that juvenile survival is independent of climate and other environmental factors;

Model B: juvenile survival varies over time and space (i.e. depends on climate, and other environmental factors) in the same way as adult survival, but has a potentially different mean value, so that:

where denotes the mean adult survival rate from Horswill & Robinson (2015).

Neither model is likely to be entirely plausible, but they offer two scenarios for how variations in juvenile survival are linked to adult survival (completely in Model B, not at all in Model A); applying both models and comparing the results will often be the most useful strategy, as without direct data on juvenile survival (which are rarely available) it is impossible to empirically select between the different models.

A full count of the number of chicks born at the colony in year was estimated from a sample, rather than observed directly from a census. We assumed that the productivity rate for sampled and unsampled birds was identical, and that the actual number of fledged chicks from unsampled nests followed the same distribution as chicks from sampled nests (i.e. binomial).

Non-Bayesian inference

Model A can be fitted using an approximate form of non-Bayesian inference. Specifically, if the number of newly recruited breeding adults at each colony in each year were known then the number of surviving breeding adults from year can be calculated, and the model can be fitted directly to abundance data as a Binomial GLMM. For each year at each colony we therefore estimated the number of new recruits deterministically, using the nest productivity data and juvenile survival rate, via the following steps:

Step 1. Estimated number of chicks born in year

= (Number of breeding pairs in year

* (Number of fledged chicks in sampled nests in year )

/ (Number of sampled nests in year )

Step 2. Estimated number of new recruits in year

= Estimated number of chicks born in year *

Bayesian inference

Within the context of the survival models there are two key advantages in using a Bayesian approach, even though the non-Bayesian approach is computationally much faster to use, and can be applied in a more automated way (requiring less manual intervention and checking) than the Bayesian approach:

1. The Bayesian approach is more flexible, and so can be used in situations where the non-Bayesian approach is infeasible. Specifically, the Bayesian approach can be used to fit models in which juvenile survival is either assumed to either be linked to adult survival (Model B) or independent of adult survival (Model A), whereas the non-Bayesian approach can only be used to fit the latter. Model B is biologically more plausible than Model A, so it is a substantial restriction of the non-Bayesian approach that it can only be applied to Model A.

2. The Bayesian approach also accounts more fully for uncertainty. This is a general advantage of Bayesian methods, but is particularly important here, because the non-Bayesian analysis requires the strong assumption that the number of new recruits can be calculated deterministically - the number of new recruits can be treated as a stochastic random variable within the Bayesian analysis, avoiding this assumption.

The Bayesian approach that we used here is similar to that used by Freeman et al. (2014) and Jitlal et al. (2017). The key difference is that they used the model to estimate juvenile survival, in a situation where no environmental/climate effects were considered and where a good independent estimate of adult survival was available. We focused instead upon estimating adult survival, and the relationship between adult survival and climate/environment, and assumed juvenile survival was known. Note that it is impossible to estimate both juvenile and adult survival empirically simultaneously, using abundance data, because the model is unidentifiable – i.e. the same observed series of abundance values can be obtained through a range of different possible combinations for juvenile and adult survival rates.

Model selection

Model selection within the context of Model A was performed using the non-Bayesian approach, as in the analyses for productivity and growth rate. The best supported model for adult survival under Model A was re-fitted in the Bayesian framework - we fitted both possible models for juvenile survival (Model A and Model B), and compared these using deviance information criterion (DIC).

Results

Having identified the best supported set of climatic variables and seasonal definitions for adult survival for each species using glmer models, we then assessed how support in the data differed when juvenile survival was assumed to be constant over time and space (Model A, Methods), derived from species-level estimates (Horswill & Robinson, 2015), versus models in which juvenile survival was allowed to vary over time and space in the same way as estimates for adult survival, derived from the model fitting (Model B, Methods). By necessity, these models could only be fitted using the Bayesian framework. We found that in all species, support for a constant value for juvenile survival vastly outweighed support in the data for varying juvenile survival rates, in line with variation in adult survival (Table 20).

We detected some differences in the strength, and occasional differences in the direction of effects when juvenile survival was allowed to vary in comparison to models where it was kept constant (Table 21). Relationships between climate variables and adult survival tended to be less strong in models in which juvenile survival was allowed to vary (Table 21).

Models with juvenile survival as a constant

Here, we focus our report on the main results from models of adult survival in which juvenile survival was assumed to be constant over space and time, and was based on species level estimates from Horswill and Robinson (2015).

Models of survival would generally model the logit of the survival probability as a linear function of covariates, but in our analyses of adult survival we use a log, rather than logit, transformation. This was because initial attempts to fit the models using a logit transformation suggested that such models could not be fit to the SMP data - attempting to do so resulted in major issues of non-convergence that we could not resolve. This issue is likely to arise because of an irresolvable inconsistency between the assumptions of the model and the empirical characteristics of the data - e.g. the data suggest that the survival probability is sometimes above one, which is logically impossible. There are a number of possible reasons that such an inconsistency may occur:

a) the SMP counts and/or productivity data contain a higher level of observation error than our models can account for;

b) the productivity data do not always represent a representative sample from the colony, and so the productivity rates estimated from these data may provide biased estimates of the year-specific productivity rates for the entire colony;

c) the model makes biologically unrealistic modelling assumptions – in particular, both models (Model A and Model B) necessarily make strong assumptions about the relationship between juvenile survival and climate, which may not be realistic. More specifically, they either assume that juvenile survival is unrelated to climate [Model A], or assuming that the relationship between juvenile survival and climate is the same as that between adult survival and climate [Model B].

Using a log transformation provides a technical way of avoiding this issue, but introduces a degree of biological implausibility in to the model, since it implies that it is possible for adult survival to be greater than one. The results obtained with these models should therefore be treated with some caution, particularly when extrapolating (e.g., extending estimates of adult survival into the future).

Black-legged kittiwake

The best-supported model for survival in this species included effects of marine variables and terrestrial wind over the whole year (Table 22 ). This model included strong negative correlations with the North Atlantic Oscillation (NAO; posterior mean: -1.336, 95%CI: -3.647, 0.034) and sea surface salinity (SLM: posterior mean: -0.282, 95%CI: -0.321, -0.238), as well as strong positive correlations with the Atlantic Multidecadal Oscillation (AMO; posterior mean: 1.271, 95%CI: -0.108, 3.591, >92% posterior density positive), sea surface temperature (SST; posterior mean: 0.142, 95%CI: 0.133, 0.152) and terrestrial wind speed (WS; posterior mean: 0.052, 95%CI: 0.044, 0.059) (Table 22). There was some support in the data for models containing alternative combinations of climate variables (only marine variables ΔAIC 2.4; marine plus terrestrial plus wind ΔAIC 2.7; and marine plus terrestrial ΔAIC 3.1; Table 22 ). The model containing marine and terrestrial wind variables defined over the non-breeding period also received similar support to the best supported model (ΔAIC 1.19), with the null model containing no climate variables also receiving some support (ΔAIC 3.64) (Table 22 ).

Common guillemot

In this species, the best supported model included effects of marine climate variables and terrestrial wind defined over the non-breeding period (Table 22 ). This model included strong negative correlations with the North Atlantic Oscillation (NAO; posterior mean: -0.749, 95%CI: -1.227, -0.013) and sea surface salinity (SLM: posterior mean: -0.114, 95%CI: -0.157, -0.070), as well as strong positive correlations with the Atlantic Multidecadal Oscillation (AMO; posterior mean: 0.748, 95%CI: 0.018, 1.213), sea surface temperature (SST; posterior mean: 0.018, 95%CI: -0.003, 0.038) and terrestrial wind speed (WS; posterior mean: 0.023, 95%CI: 0.014, 0.032) (Table 22). Models containing alternative combinations of climate variables received little support in the data compared to the best supported model (ΔAIC: 3.8 for all other combinations, Table 22 ), as did the null model containing no climate effects (ΔAIC: 3.78). Similarly, when marine climate variables and wind were defined over the whole year, the model received little support in the data comparted to the model where climate influence was defined over the non-breeding period (ΔAIC 3.66; Table 22).

Great black-backed gull

Results for this species indicated poor model fit (Table 22), suggesting the information in the data was insufficient to estimate model parameters and effects robustly. We therefore do not report further on the model for this species.

Herring gull

The best supported model for this species included effects of marine climate variables and terrestrial wind speed defined over the entire year (Table 22). The model included strong negative relationships between adult survival and NAO (posterior mean: -6.195, 95%CI: -13.545, -0.818) and sea surface temperature (posterior mean -0.375, 95%CI: -0.787, -0.070), and a strong positive relationships with AMO (posterior mean: 6.705, 95%CI: 1.359, 13.782) (Table 22). A model containing the same climate variables, but defined over the non-breeding season only, received essentially identical support in the data (ΔAIC 0.007; Table 22). Models containing other combinations of climate variables received little support in the data (ΔAIC 3.1 for all combinations; Table 22), as did the null model with no climate effects (ΔAIC 3.16; Table 22).

Razorbill

Adult survival in this species was best explained by a model including marine climate variables and terrestrial wind defined over the non-breeding period (Table 22). This model included strong negative relationships with NAO (posterior mean: -0.440, 95%CI: -0.609, -0.271) and sea surface salinity (posterior mean: -0.443, 95%CI: -0.567, -0.319), and a strong positive association with terrestrial wind speed (posterior mean: 0.095, 95%CI: 0.078, 0.116) (Table 22). Other models received considerably less support in the data, with the next best model including marine, terrestrial and wind climate variables (ΔAIC 3.5), and other models receiving even less support in the data (marine variables only ΔAIC 6.8; marine and terrestrial variables ΔAIC 7.0), as did the null model with no climatic effects (ΔAIC 6.96). A model containing the same climate variables as the best supported model, but defined over the entire year also received little support in the data (ΔAIC 6.01).

Discussion

We fitted two alternative broad model formulations to estimate relationships between adult survival and climate variables. In the first of these (Model A), juvenile survival was assumed to be constant over time and space, with the value of juvenile survival derived from species-level estimates (Horswill & Robinson, 2015). In the second formulation (Model B), we allowed juvenile survival to vary over time and space in the same way in which estimates for adult survival were estimated to vary over time and space. Overwhelmingly, support in the data for models with constant juvenile survival outweighed that for models with varying juvenile survival. This result was surprising, given that ecologically we would expect juvenile survival to vary over time and space, and the model parameterisation forcing juvenile survival to vary in a similar pattern to that estimated for adult survival also follows from current ecological understanding of a degree of correlation in variation in survival rates across different age classes. We therefore focused our more detailed assessments of the influence of climatic factors on adult survival to those models in which juvenile survival remained constant over time and space.

We were unable to fit survival models for Atlantic puffins due to insufficient data on both abundance and productivity. Similarly, for two additional species, European shags and northern gannet, the best supported survival model contained no climatic effects, suggesting available data for these species was insufficient to detect the influence of climate on adult survival. Finally, one more species, great black-backed gull, demonstrated support in the data for a model including marine and terrestrial climate variables, but model fit was exceedingly poor and model estimates were implausible, again suggesting that in this species there was insufficient data to properly detect and estimate relationships of climate with adult survival.

Black-legged kittiwake and common guillemot showed very similar responses of adult survival to climate. In both species, adult survival was related to marine and terrestrial wind variables, although they differed in the seasonal definition used to aggregate climate impacts – with black-legged kittiwakes responding most to climate defined over the whole year, and common guillemot responding most to climate defined over the non-breeding period only. Adult survival in both species was lower when NAO and sea surface salinity were higher, and was greater when AMO, sea surface temperature and wind speed were higher. Projected future climate predicted higher future adult survival rates for both black-legged kittiwake and common guillemot, although there will little evidence that either species could substantially alter its adult survival rate through extending foraging range around breeding colonies, either now or under future climate scenarios.

Razorbill showed similar responses of adult survival to climate as was detected for black-legged kittiwakes and common guillemots, with adult survival being most strongly correlated with marine climate variables and terrestrial wind, as defined over the non-breeding period. This species also showed lower adult survival rates when NAO and sea surface salinity were higher, and higher adult survival rates when terrestrial wind speed was greater. Projected future climate also resulted in higher predicted adult survival rates for this species, again with no opportunity for this species to augment adult survival rates by extending foraging range around breeding colonies, either under current or future climate conditions.

Finally, adult survival in herring gull was also most strongly related to marine climate variables and terrestrial wind, as defined over the entire year. This species also showed lower adult survival when NAO was greater, and lower adult survival with higher sea surface temperatures. However, in contrast to the other species, herring gull showed higher adult survival when AMO was greater. Adult survival in herring gull was predicted to decline under future climate conditions, although there was considerable uncertainty around these estimates. As with other species, there was little to no predicted opportunity for this species to increase adult survival by extending foraging range around breeding colonies.

The fitted models, when applied to future climate conditions, often predicted adult survival rates of greater than one, which is clearly biologically unfeasible. Typically statistical models of probabilities using logit functions to ensure estimated rates lie between zero and one, however when applied to the SMP count and productivity data, these models did not converge, likely due to irreconcilable inconsistencies between the assumptions of the model and the empirical characteristics of the count and productivity data. For instance, there are various combinations of counts of breeding pairs and fledged young in some years which imply the adult survival rate is sometimes greater than one. In all survival models, we therefore had to use a log transformation to achieve model convergence. This parameterisation avoids non-convergence issues by relaxing the constraint that adult survival estimates must always lie between zero and one, allowing inconsistencies between data and modelled processes to be avoided. It does, of course, however introduce a degree of biological implausibility into the model. As such, the results obtained from these models, particularly the extrapolation of survival rates under future climate, should be treated with caution.

The underlying mechanism explaining the positive effect of wind on survival is not clear, and therefore warrants further investigation. Many climate models predict future changes in wind speeds. A key question revolves around the importance of effects across the range of speeds experienced by the birds which we considered here, and the frequency and severity of extreme events, which are predicted to increase in future in many regions (Rahmstorf and Coumou, 2011; IPCC 2018), because the effect of wind speed may be non-linear such that there is a positive effect overall but a negative effect of high winds (Frederiksen et al. 2008). Our results also show that temperature is having a strong positive effect on survival of black-legged kittiwakes and common guillemots. This relationship is in stark contrast to Frederiksen et al. (2004) who demonstrated a negative effect of temperature on survival in kittiwakes. Given the widespread evidence that warming is having a negative impact on prey of seabirds, this appears counterintuitive, but one possibility is that these potential impacts are more than compensated for by the reduced energetic costs of key activities such as foraging and resting that seabirds, as endotherms, experience at higher temperatures (Amelineau et al. 2018). The predicted large declines in productivity under future climate projections may also in part explain why these models predicted counter-intuitive increases in survival in relation to climate variables such as temperature.

We used count and productivity data to estimate survival because it was the only approach that would allow us to study multiple populations and species. Mark-recapture data are only available for a few studies, and therefore do not allow for a UK-wide multi-species assessment of relationships between adult survival and environmental variables. However, survival estimates derived from mark-recapture data are a significantly more powerful and reliable approach for estimating drivers of change in survival. For example, using that approach, Frederiksen et al (2008) demonstrated the effect of wind on survival of European shags, whereas we could not detect any climatic effects in this species. Of more concern are cases where we found opposing results than published survival studies using mark recapture, such as the effects of warming on survival of black-legged kittiwakes. This, together with predictions of survival in our models that exceeded one, leads us to conclude that estimating survival from counts and productivity is not reliable where there are significant gaps as is the case with the SMP.

A further challenge is that for several species in this study – notably black-legged kittiwake and northern gannet – a considerable proportion of the adult population spends the winter outside UK waters. As such, the environmental variables used here may not be particularly relevant to survival prospects, since most adult mortality occurs at this time. Thus, incorporating environmental drivers at wintering grounds would potentially have provided important insights (Reiertsen et al. 2014).