Developing a population model for Rum Manx Shearwaters for assessing offshore wind farm impacts and conservation measures

3 Methods

3.1 Data availability

Our survey of available data prioritised use of peer-reviewed results. Since we aimed to understand how the population functions, rather than just producing a current snapshot of its composition and size, we considered historical information to be just as valuable as recent data. Information extracted from publications was of a numerical nature (estimates, confidence intervals, timing and extent of effort), but also descriptive. Descriptive information provided clues about the reliability of these numbers and the possibility of their integration into the study. For example, any indication from the methods description that not all sources of uncertainty had been included in statistical estimation would lead us to treat the confidence interval in the estimate as potentially too narrow. Unpublished reports, working papers, consultations, published government reports or raw data were used as sources of the population’s most recent status. It is particularly relevant that several of the observations provided for the same biological variables have been obtained via different methods, so they are not directly comparable without some form of scaling or calibration.

3.2 Survey data

Manx shearwaters are burrow-nesting or cavity-nesting seabirds that only return to the colony at night. This makes them notoriously difficult to count (Murray, et al. 2003, Bolton and Thompson 2023). Rum has been difficult to survey accurately and so the population size and the direction and shape of the population trend are rather uncertain. Ringing-recapture methods are generally considered ineffective for estimating population size (Wormell 1976) and therefore the population size is estimated by a calculation involving observations on apparently occupied burrows (AOB) per unit area surveyed and an estimation of the overall area of the colony. The density of burrows appears to remain consistent across surveys but the area of ground occupied is highly uncertain and estimates differ considerably so that the total numbers depend primarily on the estimated occupied area of the colony.

It is relevant to note that Manx shearwater burrows at Rum tend to be long (often about 2m from entrance to nest chamber). The ground in which the burrows are dug is very hard gravelly sand and rock. Burrows are often dug under or around loose rocks that are too large for shearwaters to move. Excavating these burrows must take a lot of time, and probably is the effort from several years of digging before a burrow is fit to use for breeding. This means that growth in breeding numbers may be constrained by the number of burrows, and especially by the number of high-quality burrows that provide protection from flooding and from predators/disturbance. Similarly, if the population were to decline, burrows would not disappear quickly, but would remain obvious in the terrain for many years. Therefore, a decline in breeding numbers would result in a higher proportion of unoccupied burrows rather than a large decrease in the number of burrows. These points are relevant to consider in relation to the large variations in estimates of colony size from surveys carried out between the 1960s and 2020s.

Relevant data identified included:

Population monitoring and area mapping study from 1960s-1970s (Wormell (1976)): Carried out detailed mapping and characterisation of the vegetation associated with suitable burrow habitat (the shearwater greens). An estimate of the total area of shearwater greens was made by aggregating all the greens from vegetation surveys carried out between 1961 and 1964 with the help of aerial photographs taken in earlier years. The density of burrows was surveyed in different compartments on the island in 1965 to 1969.

Population monitoring study from 1979 (Thompson and Thompson (1980)): Used the same overall approach and calculation of colony area as (Wormell 1976) with a different set of random quadrats for collecting data on burrow density and occupancy. They found a mean burrow density 21% higher than in any other survey, a result that may easily have occurred if their quadrats tended to be placed in the more densely populated shearwater greens (Murray et al. 2003).

Population monitoring study from 1982 (Philips 1982): Used the same calculation of colony area as (Wormell 1976) but with a different set of random quadrats for collecting data on burrow density and occupancy, again following essentially the same survey protocol. Burrow densities estimated by (Philips 1982) were lower than those in (Wormell 1976) and (Thompson and Thompson 1980).

Population monitoring studies from 1990 & 1995 (Furness 1990, 1997): Two surveys conducted using the same methodology (105 quadrats, 10 mx10 m size). Burrows were classified into used and unused. The proportion of burrows assessed as unused in the 1990 and 1995 surveys was higher than seems to have been the case in (Wormell 1976). The estimated area of the colony was similar to the area derived by (Wormell 1976). The confidence intervals of the later surveys were considerably larger because they included uncertainty in the estimation of areas (as well as sampling uncertainty in the estimation of burrow density).

Population monitoring study from 2001 (Murray et al. 2003): Estimated occupancy rates of burrows based on visual cues at burrow entrances (droppings, feathers etc.) and on responses to tape playback of Manx shearwater calls. An estimate of colony extent was based on the total area of ground above the 457 m altitude contour. The area considered was dramatically larger (4.5 times) than the previous estimates of the colony extent. Estimates of population size were consequently also higher, but not proportionately so (2 times). The estimate was derived from a stratified survey distinguishing between high-density core areas and lower density peripheral areas. A parallel method of data collection using tape playback was undertaken and a correction factor of approximately 1.9 was applied to account for non-responding birds. This may have been a low correction since a similar study using visual and acoustic methods carried out on Welsh populations (Smith et al. 2001) found that response rates are between 0.3 and 0.53 (implying corrections in the range 1.9-3.3).

MarPAMM from 2021 (Inger et al. 2022): The report by Inger et al. (2022) estimates a much higher number of Manx shearwaters at Rum than previous counts. The estimate is based on a very similar mean density of burrows within quadrats but a very much larger estimate of colony area (characterisation of the colony area was based on habitat suitability modelling, rather than direct survey). Statistical concerns about this estimate are the lack of validation/ground-truthing of the habitat model, and the partial propagation of uncertainty to the end results.

3.3 Demographic data

We considered the following aspects of demography:

Life history: The earliest age to reproduction, as inferred by early studies (Harris 1966, Brooke 1977) on the rates of return by young shearwaters is approx. 5 years. These data indicate the possibility of birds being found for the first time in later years and this could be an artefact of detectability, or a result of postponed breeding due to density-dependent limitations (see section on density dependence below). A discussion on the limitations of this estimate (Horswill and Robinson 2015) also includes the possibility that some of the birds in these early studies had actually bred undetected at age 4.

Dispersal: There is some evidence that the levels of breeding dispersal in the species are low (Harris 1966, Brooke 1978, Horswill and Robinson 2015). This permits the assumption of a closed population, allowing us to attribute all adult losses to mortality, and also to assume that all new recruitment comes from recruitment of juveniles. It is less clear whether cross-subsidies of juvenile recruits from nearby colonies are substantial. However, levels of productivity in the most proximate colony at Canna are lower than at Rum (Thompson and Furness 1991), making Rum a source population (hence, one donating surplus individuals that are unlikely to have an impact on a model of a closed population).

Breeding success: Baseline information on the rates of breeding can be derived from estimates of breeding propensity (Newman et al. 2021, Wood et al. 2021) as well as studies carried out elsewhere (Brooke 1977, 1978, Mavor et al. 2006, Booker et al. 2008). Time series of breeding success (counts of nests with successful reproductions) come from multiple reports for Rum (Thompson and Furness 1991, Wood et al. 2021) and the neighbouring colony at Canna (Thompson and Furness 1991). The time series on breeding success on Rum spans the years 1958-69, 1971, 1973, 1985-86, 1994-2018. Sampling effort varies from ~30 nests in the early years, up to 221 in 2014. Data on breeding success of Manx shearwaters on Canna ((Thompson and Furness 1991) and JNCC Seabird Monitoring Programme database, now hosted by BTO at https://app.bto.org/seabirds/public/index.jsp) exist for the period 1973-81, 1982-85. These can be used as a relative index of productivity for the missing years of data on Rum, but not a simple proxy since productivity on Canna is generally lower.

Adult survival: There is good baseline information on adult survival from various sources (Perrins et al. 1973, Brooke 1977, Horswill and Robinson 2015, Horswill et al. 2016, Wood et al. 2021). A previous analysis of ringing data (Duff 2011) has yielded a continuous time series of estimates with associated confidence intervals for the period 1994-2010. This is used as a product (i.e., treated as data for the state-space model, with associated uncertainty).

Juvenile survival: There is only baseline information for juvenile survival from three publications (Harris 1966, Perrins 2014, Horswill and Robinson 2015). The estimates are derived from different viewpoints. The Perrins (2014) study provides an overall survival estimate for the fledgling and juvenile years, whereas the Harris (1966) study provides an estimate for survival of 4y old juveniles. These estimates were used jointly by the model to inform a prior for 1^st year survival. There is no information on annual variations in juvenile survival.

3.4 Covariate data

A key consideration regarding the availability of covariate data is that we require as complete a time series of covariates as possible. Some intermittency can be addressed via data imputation by the model, but we should weigh the relative value of information contributed by a fragmentary covariate against the inferential cost of the additional degrees of freedom introduced by imputation. In this work, we examined an autocovariate (density dependence) and an external covariate (rainfall). In the discussion we examine future possibilities for inclusion.

Density dependence: Thompson (1987) discusses possible mechanisms of density-dependent regulation. Her findings suggest that the availability of high-quality burrows (as defined with respect to vulnerability to flooding, see below) may act as a temporary or permanent barrier to population growth, mainly via its regulation of breeding success. It is not generally clear if density dependence in seabirds operates at the colony (i.e., through competition for nest sites e.g. burrows) or at sea (i.e., through competition for food), but either of these may have a suppressing effect on reproduction.

Rainfall: There is evidence (Thompson 1987, Thompson and Furness 1991) that heavy rainfall events during the incubation period can impact the breeding success of birds at Rum. A monthly rain gauge at the colony was maintained for a number of years by the National Nature Reserve (NNR) team. Additional data from a nearby area (Tiree; the closest continuously monitored site but about 60 km from Rum, but a low-lying island with much lower rainfall than in the hills of Rum) were used as complementary information to help with imputation of missing data on Rum.

In summary therefore, the available data for the model were:

1. A total of seven surveys designed to estimate population size (1976, 1979, 1982, 1990, 1995, 2001 & 2021), each with a reported 95% confidence interval, mainly relating to sampling variation. Each survey has a different method to estimate the size of the colony and the estimated colony area is also known (but subject to debate).

2. A single, additional tape-playback survey in 2001 to obtain an estimate of bird abundance via a radically different method to the survey inspection.

3. Direct or indirect prior knowledge about the four demographic baseline processes of interest (survival of adults, juveniles and fledglings as well as fecundity).

4. Several continuous years of data on productivity of sampled burrows. Different studies have observed different sample sizes and have applied slightly different definitions of breeding success estimated by various methods (e.g., frequency of nest burrow checks and means of assessing burrow contents). The period 1973-1993 is relatively sparse in such data from the Isle of Rum.

5. A contiguous period of estimates for adult survival (1994-2010), with associated standard errors.

6. A period of years between 1973 and 1985 where productivity was recorded for the neighbouring island of Canna.

7. Time series of rainfall data, available for different times of the year. Only a limited number of years exists for Rum itself, and more extensive uninterrupted data from other locations, such as Tiree. An integrated cross-calibration exercise of the rainfall data is therefore needed, to impute the missing Rum rainfall data.

8. Very short (5y) period of rat abundance index, not overlapping with data on breeding success.

3.5 Modelling

The key modelling challenge of this study is in distinguishing between true ecological fluctuations and artefacts created by the changes in survey methodology and calculation. Strong points of the data set are the auxiliary measures of demographic rates (survival and fecundity) and body of expert knowledge from field ecologists who have been dedicatedly working on this population for several decades.

The first priority was to develop an integrative model that could accommodate all of the available demographic data, have capacity to evaluate multiple covariates and be supported by expert views. We adopted a Bayesian Integrated Population Modelling (IPM) approach (Buckland et al. 2007, Fieberg et al. 2010, Maunder and Punt 2013, Matthiopoulos et al. 2014, Zipkin and Saunders 2018). In particular, given the dynamical nature of the questions involved, we used a state-space model, which couples a model for the biological (i.e., generative) process together with a model for the observation (i.e., data collection) process.

Models are not solely informed by data inputs and expert opinions, but also by biological constraints. Here, these originate from the key principles of population dynamics. A key priority of our modelling therefore was to recognise such mechanistic principles and use an approach that explicitly incorporates constraints.

A further priority was to correctly address uncertainty in reconstructed parameters, system states, forecasts and counterfactuals. This requires us to tread a fine balance between acknowledging uncertainty of different forms at the information-input stage, but without being so agnostic that we end up ignoring basic biological principles and cause the model to fail (by e.g., not converging to an answer) because of lack of information.

Achieving the first modelling objective (data integration) can be done via a frequentist approach, but it is generally easier to achieve on a Bayesian framework. The Bayesian approach also affords a more flexible treatment of uncertainty in model predictions and is, by design, suited to the incorporation of expert opinion in the form of priors.

The model's state variable was the number of breeding females, or, alternatively, the number of active burrows in the population. We modelled population dynamics in discrete time at the annual scale. In principle, next year's population size will comprise surviving adults from this year and the new recruits (the survivors from the cohort of chicks produced 5 years ago). This was expressed as a stochastic model to allow for demographic stochasticity and random effects to absorb the effects of unknown environmental covariates or errors in the observation process.

The process of breeding success was, in later versions of the model written as a function of rainfall, and its effect on population size was lagged by five years, to allow for the age-to-maturity in Manx shearwaters. Recruitment of maturing juveniles was constrained to be negatively dependent on current population density (compensatory density dependence (Miller et al. 2019)). In the project’s Technical Report, we investigated theoretical extensions of this rudimentary form of density dependence to allow for the imbalances in processes of burrow-building and natural degradation.

We implemented stochastic (overdispersion) terms in the linear predictors of the model's demographic rates and allowed the model to estimate the extent of overdispersion in each demographic process. This is a good indicator of which process is currently most vulnerable to extrinsic variations, and therefore which process should be getting investigated in the future for inclusion of candidate covariates.

The recorded data are as much a result of biological processes, as they are of observation artefacts and effort. In this sense, state-space models require equivalent amounts of modelling effort for the specification of the observation component. A large part of integrated population modelling, particularly when the time series available are not complete, or overlapping, is the imputation of missing data. In some cases (e.g., reconstructing the time series of demographic rates in this model) the results of imputation are interesting in their own right. In other cases (e.g., filling-in gaps in covariate data), imputation is necessary to deal with nuisance variables or parameters.

Population surveys: All surveys recorded the density of burrows in selected regions and then scaled up to total population size by calculating the total area of the colony on the island. Different surveys took different approaches in the field, but also implemented a differing scaling-up calculation to derive population size. Two types of error are built into the given estimates. First, bias in the point estimate. Underestimates could result from consistently selecting low-density areas for the sample surveys. Overestimates could result from assuming that more areas are suitable, or currently contain burrows. Second, bias in the precision. Accounting for only some sources of uncertainty in the estimate could have led to a smaller-than-necessary CV. Being too precautionary about the sources of error could have the opposite effect. We addressed such issues with a critical evaluation of the field and estimation limitations of different surveys and an extensive sensitivity analysis of model performance on different levels of bias (see below).

Demographic data: Observed proportions of burrows surveyed that reproduced successfully were modelled while accounting for variable observation effort. Mark-recapture estimates of adult survival (Duff 2011) over the period 1994-2010 were used, together with their associated standard errors. There were no survival estimates for the fledgling or juvenile stages in the model, so the only information provided directly on those demographic rates was in the form of baseline priors. The full trajectories for all four demographic rates were reconstructed as latent variables by the model.

Prior specification: The key parameters of the process model that were supported from prior knowledge are the baseline demographic rates (baseline adult survival and juvenile survival). We used a precautionarily broad range for the priors of these parameters. Although we don't have direct information for first-year survival, an informative joint prior was constructed from the available information on compounded survival from birth to recruitment. Breeding success was well informed by the available data, so it was not advisable to constrain the intercept arbitrarily, especially since all the information available for the baseline fecundity comes from the fecundity data that are already used by the model.

Addition of data from playback surveys: During the 2001 survey (Murray et al. 2003), a parallel method of data collection using tape playback was implemented at a time of year when non-breeding burrow occupiers would not be an issue. A correction factor of approximately 1.9 was applied to account for non-responding birds. This may have been a low correction since a similar study using visual and acoustic methods carried out on Welsh populations (Smith et al. 2001) found that response rates are between 0.3 and 0.53 (implying corrections in the range 1.9-3.3). To make use of this survey datum, we inflated the estimated confidence interval for the call-back population estimate by accounting for the uncertainty in the response rate. This led to a mean estimate of 97,945 and a CV of 0.175. It is worth noting that this calculation is over the same estimate for the area of the colony provided by (Murray et al. 2003) for the visual part of their survey.

Addition of breeding data from Canna: The data from Canna have two distinct advantages. First, they cover a time period characterised by data-sparsity on Rum. Second, they have a one-year overlap with breeding data on Rum, allowing at least a rudimentary calibration. Despite its proximity to Rum, Canna offers a different geomorphology and vulnerability to risks such as rats and flooding. We therefore treated the Canna breeding success data as a relative index of breeding success on Rum. The model estimated the scaling relationship between Canna and Rum, as a nuisance parameter.

3.6 Fitting to real data

The first version of the model (Model 1) was fitted to the real data, and henceforth treated as a baseline for further development. Model 1 contained no covariates (just generic random effects) and assumed that the coefficients of variation in survey counts were exactly as reported in the corresponding studies.

3.7 Validation with simulated data

Simulated data offer a sanity check for the performance of the model. General practice is to generate synthetic data from the model to be used for fitting to the real data, using plausible parameter values, and then proceed to fit the model to these data. Given that there is, by construction, no model misspecification in this procedure, assuming that the synthetic data provided are adequate, the model fitting should be able to retrieve the underlying parameters and reconstruct the true demographic trajectories. The specifications of the key model parameters were selected to be the median values retrieved from the real data in model fitting. Variations in data collection effort were set to be identical to the real data.

3.8 Sensitivity on population estimates

A particularly difficult aspect of the above analysis is that it relies heavily on estimates of population size (and associated confidence intervals) that, as well as being irregular and intermittent, have been generated by different field methodologies and estimation methods. An apparently strong deviation from the established estimation approaches is the most recent population estimate obtained by Inger et al. (2022). This led to an effective doubling of the population estimate due to the analytical approach used to derive apparently suitable habitat for nesting. The resulting number is quite influential for the modelling undertaken here (particularly for the forecasts and future scenarios) and therefore the analytical approach of Inger et al. (2022) was investigated closely in the Technical Report. Particular strengths of the approach are its novel use of model-based components that employ species-habitat modelling frameworks to derive suitable areas for the shearwater burrows. The decision to quantify uncertainty at parts of the analysis is good and the logical sequence of the analysis steps is unambiguous and well justified. The approach also has certain limitations primarily to do with underrepresentation of uncertainty and the lack of a validation test.

Due to these shortcomings, we carried out a sensitivity analysis of the results from Model 1. Given the (numerical and methodological) divergence of the 2021 population estimate from historical surveys, and also, considering that its high value might be generating over-optimistically increasing trends for the population, we were interested to see if the increasing tendency predicted by Model 1 would remain if more modest numbers were used. We were also interested to see whether any of the results from Model 1 were an artefact of overconfidence in the complete time series of population estimates. As part of the sensitivity analysis we examined versions of the model that either adjusted the 2021 estimate downwards, or increased the uncertainty in all historical population estimates.

The major difference between the Inger et al. (2022) estimate and previous efforts is in the total area considered to make up the breeding colony (209 ha). Murray et al. (2003) considered the extent of the colony to be 148 ha. By scaling their burrow density and response rate estimates to a breeding area of 148 ha, Inger et al. (2022) estimated that the number of AOBs would decrease to a total of 134,514 (95% HDI: 85,122–212,886), a comparable estimate to the 2001 estimate of 119,950 (95% CI: 106,730–133,500) AOBs. Earlier surveys, before the 2000/2001 survey, considered the colony to occupy a much smaller area – or at least focused only on the obvious shearwater greens.

Model 2: To explore the sensitivity of our results to the reported 2021 estimate, we refitted the model by using the alternative, low estimate calculated by Inger et al. (2022). This experiment was aimed at capturing the effect of bias in the estimate and was considered a worst-case scenario in terms of future forecasts.

Model 3: As an additional experiment, we considered the possibility that the true population size is in-between the best and worst-case population estimates for 2021. To capture that possibility, we used the average population estimate (180,262 AOB) from the two values reported by Inger et al. (2022). To acknowledge the fact that Inger et al. (2022). had not propagated several sources of uncertainty into their final confidence interval for population size, and therefore they are likely to have underestimated uncertainty, we took the precautionary approach of deriving a CV based on the upper CI extreme of the high population estimate (403,915 AOB) and the lower CI extreme of the low population estimate (85,122 AOB).

Model 4: A doubling of CVs. This was comparable to the inflation of the 2021 CV between Model 1 (CV=0.154) and Model 3 (CV=0.442, an inflation of 2.87).

Model 5: A scaling of CVs by an order of magnitude.

Comparison between these models (see Results) led us to select Model 3, for use in examination of covariates and generation of forecasts.

3.9 Examination of covariates

In addition to the density dependence autocovariate, which was included in all models, we developed a detailed examination of the effect of rainfall on breeding success. We used the aggregated rainfall data for the period May-June from Tiree. To account for differences in rainfall between Tiree and Rum, we also used a partial and intermittent time series of equivalent rainfall measurements for Rum. Fecundity was directly linked to Rum rainfall, and the missing values from Rum were imputed simultaneously by linking Rum to Tiree rainfall by a stochastic process.

3.10 Forecasts

A population simulation using model 3 without covariates was iterated as many times as the number of parameterisations (5,000 iterations were used to create stable credible intervals for the predictions), recording the population projections for the three component classes (Adults, Recruits and Fledglings). We recorded results for a time horizon of 100 years and derived metrics of risk at 25 and 100 years into the future. Risk was defined as the probability that the population will be below its current size in a pre-determined time horizon of 25 or 100 years.

To facilitate current and future investigation, any impact was characterised by the three key characteristics of 1) Intensity (i.e., proportional losses in the ability to survive or reproduce), 2) Temporal pattern (Here, distinguishing between press and pulse perturbations (Bender et al. 1984)) and 3) Demographic effect (i.e., the demographic rate on which the impact was felt).

The code generating these projections has been optimised for speed and ease of use, so that multiple scenarios can be examined rapidly. To illustrate the type of result that can be generated for future populations, we examined scenarios of impact on adult survival and fecundity.