1 Introduction

Making inferences on population trends of highly mobile animals inhabiting variable environments is challenging. Accurate assessment of population trends in long-lived, nomadic waterbirds capable of moving long distances requires assessment at a continental scale (Johnston et al. 2015). The Eastern Australian Waterbird Survey (EAWS), consisting of 10 survey bands sampling c. 11% of wetlands within a c. \(2.7 \times 10^6\) km\(^2\) survey region of eastern Australia, theoretically makes such estimation possible. The EAWS was designed by the late Graeme CaughleyFootnote 1, to enable the estimation of total abundance across the sample area (methods as described by Caughley 1977) (see Sect. 2.1 for details).

To date, however, the EAWS has not been used to estimate the total abundance of any species of waterbird within the sampled region of eastern Australia, despite some analyses suggesting this to be the case (e.g. see Figs. 2–4 in Kingsford and Porter 2009). This is in contrast to North America, where estimates of breeding populations derived from large-scale surveys date from 1955 for some species (e.g. mallard Anas platyrhynchos) (Nichols et al. 1995). In Australia, the most recent and detailed analysis of the EAWS data examined trends in counts between the largely unregulated Lake Eyre Basin and the more highly regulated Murray-Darling Basin (Kingsford et al. 2017), despite populations being known to be highly connected between these river basins (Roshier et al. 2002; Kingsford et al. 2010). Counts from wetlands within survey bands outside these basins, including coastal wetlands, were ignored. Most Australian waterbirds are essentially nomadic and use wetlands across eastern Australia and beyond. For example, depending on prevailing conditions, very large populations of many species may occur in Western Australia (Halse et al. 1998, 2005), the Alligator Rivers regions (e.g. Morton et al. 1993) and southern Papua New Guinea (Halse et al. 1996). To analyse population trends, there is a need to make inferences based on the estimated total population, including the unobserved component. This requires analysis of the EAWS survey data that accommodates vagaries such as population movement.

Fig. 1
figure 1

Estimated population size (solid points and solid lines) and associated 95% confidence intervals either from the Ratio Method (yellow shading) or bootstrapped (grey shading) of glossy ibis (Plegadis falcinellus) calculated using the Ratio Method

Full size image

Sampling-based methods applicable to estimating total abundance from the EAWS data are detailed in Caughley (1977), who noted the key role of the book “Sampling Techniques” by Cochran (1963). In particular, the “Ratio Estimate” described by Caughley (1977) is at first glance appropriate for the transect design of the EAWS (see Eq. 7 and 8 in Appendix A for details). However, a preliminary application of the method points to a possible reason as to why the approach has not been used. Figure 1 illustrates a problem in relation to applying the method to EAWS counts for the glossy ibis (Plegadis falcinellus) (see Sect. 2.1 for survey details), whereby the 95% design-based confidence intervals for the estimated population are less than zero in a high proportion of survey years. This is not surprising, given the variance formula is a large-sample approximation (noted by Caughley 1977), and the small sample size (\(n=10\) bands) (Fig. 1). Applying the Ratio Estimator to bootstrapped counts (now easily applicable, but not during the early 1980s; see Manly 2007) produces better-behaved confidence intervals in that negative values are avoided, but the coverage is poor (too narrow), particularly in years with low counts (Fig. 1). These problems repeat for a substantial proportion of the other species. There is also a need for regularization, as some of the changes in the estimated population size are biologically implausible, such as a 50-fold increase in population size between 1989 and 1991 (Fig. 1).

Sampling variation, amplified by the ability of waterbirds to move in response to variable wetland condition, is one explanation for the high variation in counts, but there is also a structural design issue that needs addressing. The EAWS sampling design essentially involves two strata (“wetland” and “non-wetland”) (see Sect. 2.1), but with the important difference from the usual approach to stratified sampling being that the non-wetland stratum is not sampled. There are arguably good logistical reasons for this, on the basis of the vastness of the non-wetland stratum and expected sparseness of the counts, and the sampling design would seem appropriate for species that are strongly tied to wetlands, such as fish-eaters (e.g. the Australian pelican Pelecanus conspicillatus), or aquatic herbivores (e.g. the black swan Cygnus atratus). However, such a sampling design creates unique challenges for the estimation of abundance for the many waterbird species that are not obligate users of the wetland stratum (as defined by the EAWS). For example, straw-necked ibis (Threskiornis spinicollis) spend considerable periods of time away from wetlands foraging on terrestrial prey, though predominantly use wetlands to breed (typically in colonies), while grey teal (Anas gracilis) are noted for their use of shallow ephemeral wetlands (Roshier 2009). Grazing species of waterfowl such as wood ducks (Chenonetta jubata) are increasingly using small farm dams as breeding and roosting habitat (Kingsford 1992) that do not meet the minimum one hectare threshold to be considered as part of the wetland stratum. Clearly, methods for estimating the total abundance of these more problematic species need to also model their choice of habitat, and in the current context, how this relates to surveyed and non-surveyed strata. This is somewhat analogous to the “temporary availability” issue within a survey plot as modelled by Chandler et al. (2011), who also use the term “superpopulation” (after Kendall et al. 1997) to describe the total number of individuals that could possibly be present on the plot surveyed.

The EAWS is an “observational experiment”, subject to variation and trend in rainfall, the most important driver of waterbird dynamics and movements through the provision of wetlands (creation and expansion, filling and drying) and ultimately food resources (Crome 1986, 1988). Rainfall needs to be accounted for, such that we can estimate population growth rates under specified (e.g. average) conditions (recognising the uncertainty that a warming climate may bring), rather than marginalised over the conditions during the study period. This enables inferences on whether species are undergoing population decline driven by underlying negative growth rates. This is particularly pertinent in the Australian context where the high inter-annual variability of rainfall (Morton et al. 2011) makes it unlikely that starting and finishing conditions over the survey period will be equivalent. Indeed, the early counts within the East Australian Waterbird Survey were notably high (see Sect. 3 and Supplementary Material), and this has been postulated to result from the preceding period (1970s) being notably wet, leading to increases in both wetland areas in inland Australia (e.g. Roshier et al. 2001), and waterbird populations.

In this paper we analyse the EAWS data as was originally intended: to make inference on trends in abundance of the waterbird superpopulations that the survey transects sample. Our aim is to develop a method that is generic enough to be applied to those species with sufficient non-zero yearly observations. This requires a model-based approach (see Borchers et al. 2002) that accounts for drivers and sources of variation in: (1) The underlying “state” of the population including demographic processes; and (2) The observation (counting) process and the survey design (spatial sampling process). Results of the analyses are considered in light of published information on population trends of several species, both those using the data in the EAWS (Kingsford et al. 2017) and independent sources (e.g. Australian white ibis Threskiornis molucca: Martin et al. 2010; Smith et al. 2013).

2 Methods

2.1 Data sources: the East Australian Waterbird Survey

Background: The East Australian Waterbird Survey (EAWS) began in 1983 (Braithwaite et al. 1986), and has continued annually since (Kingsford and Porter 2009; Kingsford et al. 2017). An early aim of the EAWS was to gather information on the distribution and abundance of waterfowl to better inform controls on harvesting (“duck hunting”), which at the time was increasing rapidly. The ultimate objective had been envisaged as, “the provision of data enabling the measurement of the relative effects of weather and hunting on waterfowl population and the estimate of optimal yields and appropriate levels of hunting pressure.”(Braithwaite et al. 1985). Since then, concerns regarding hunting pressure have somewhat abated following restrictions on recreational duck hunting (particularly in New South Wales), and the EAWS counts have been increasingly used as environmental indicators of wetland health, within a broader debate regarding water use (e.g. Kingsford 1999; Kingsford et al. 2017). We analysed the first 36 years of survey data, 1983–2018, accessed from the publicly available National Waterbird Survey Database (https://aws.ecosystem.unsw.edu.au/) (Kingsford et al. 2020).

Fig. 2
figure 2

Locations (blue dots) of surveyed wetlands within bands surveyed by the East Australian Waterbird Survey, with the extent of the Lake Eyre Basin (pink shading) and Murray-Darling Basin (yellow shading) shown

Full size image

Sampling design: The sampling design is described in detail by Braithwaite et al. (1986). The survey design involves ten survey bands, each 30 km wide, spaced at intervals of \(2^{\circ }\) latitude from \(20^{\circ }30'\)S to \(38^{\circ }30'\)S, starting from \(136^{\circ }\)E to the eastern coastline of Australia (Fig. 2). Within each survey band, wetlands and waterways \(>1\) hectare were identified using 1 : 250, 0000 topographic maps, for which the underlying aerial photography was dated between 1944 and 1980 (Braithwaite et al. 1986). Only this wetland stratum within bands is surveyed, which contains c. 11% of wetlands within a c. \(2.7 \times 10^6\) km\(^2\) survey region of eastern Australia. Counting is undertaken from a fixed-wing (overhead) aircraft flown at a height of 30–46 m and a speed of 167–204 km hr\(^{-1}\) (Kingsford and Porter 2009). The counting method depends on the nature of each wetland, with transect methods used for sub-sampling within extensive wetlands (Braithwaite et al. 1986). Not all waterbirds have visual characteristics and behaviours that allow reliable identification to species from the air. We acknowledge there are methods that can account for systematic variation in the counted population, but addressing this is not central to our arguments. Furthermore, the data are not readily available in disaggregated (i.e. within wetland) form.

Species selection: More than 50 waterbird species have been recorded during the EAWS. We restricted the data set to 45 individual species for which at least two-thirds of the yearly total counts were non-zero. Species that were pooled together within the dataset (e.g. “Small Waders”), due to species discrimination being difficult from the air, were also excluded from the analysis. For reporting purposes, species were classified into the following functional groups: “Ducks”, “Gulls and Terns”, “Herbivores”, “Large Waders”, “Pelicans, Darters and Cormorants”, and “Small Waders” with group membership detailed in Table 2, Appendix B.

Covariate selection: Rates of increase of game species of waterbirds in eastern Australia have been shown to be significantly related to changes in antecedent rainfall (Briggs and Holmes 1988). We used the deviation from average rainfall over the antecedent 12-month period (October-September) as a covariate to model the rate of population increase with yearly time-steps, and the deviation from the average rainfall over the antecedent 6-month period (April–September) to model habitat selection with respect to the wetland stratum. In doing so, we sought to illustrate the use of covariates that were broadly applicable to the 45 species involved, without engaging in species-level optimisation. It would be possible to tailor covariates for particular species based on their predominant geographic distribution and critical life history needs, but that is beyond the scope of this analysis, given such a large number of species. Antecedent rainfall data were calculated from 5 km rainfall grids (Bureau of Meteorology: http://www.bom.gov.au/climate/maps/rainfall) for areas of the Australian mainland bounded to the west by \(134^{\circ }0'\hbox {E}\) and to the north by \(18^{\circ }0'\hbox {S}\) using the period 1969–2018 (50 years) as a baseline. The rainfall data features runs of above and below average conditions in response to La Niña and El Niño conditions (see Fig. 1 within Supplementary Material). Prior to the waterbird surveys starting, the 1970s featured 7 consecutive years of above average conditions (including the 1974 floods), whilst the 2000s featured 7 consecutive years of below average conditions (“The Millennium Drought”, amongst the most severe on record). We return to the impact of different choices of baseline on the results in Sect. 4.

2.2 Model detail

Considering each species as a separate analysis, let:

  • \(N_t\) = the size of the superpopulation that is sampled by the survey in year t,

  • \(n_t\) = the number of individuals counted on all bands in year t (we return to this in Sect. 4),

  • \(\phi _t\) = the proportion of the superpopulation within the surveyed bands in year t (following the nomenclature of, Chandler et al. 2011),

  • \(r_t\) = the exponential rate of superpopulation increase in year t,

  • \(x_{1t}\) = deviation of antecedent 12-month rainfall for eastern Australia in year t from the mean over the period 1969–2018.

  • \(x_{2t}\) = deviation of antecedent 6-month rainfall for eastern Australia in year t from the mean over the period 1969–2018.

State model: The unknown (latent) superpopulation \(N_{t}\), grows exponentially with log-normal process noise. The exponential rate of increase in year t is linearly related to the anomaly in antecedent rainfall in the current year (Eq. 1), where \(r_0\) can be interpreted as the rate of increase of the population under average rainfall conditions (i.e. \(x_{1t}=0\)). For a previous example of this approach to rescaling the covariate see Hobbs and Hooten (2015). Process noise (W) around population growth is log-normally distributed (as \(r_t\) is already on a logarithmic scale). That is:

$$\begin{aligned} \log (N_{t})= & {} \log (N_{t-1}) + r_t,~~\text {where} \nonumber \\ r_t= & {} r_0 + \beta _1 x_{1t} + W,~~\text {and} \nonumber \\ W\sim & {} N(0,\sigma _W^2). \end{aligned}$$
(1)

We did not fit a growth model that includes carrying capacity (e.g. a Ricker type), as in such a variable environment it is hard to justify a constant carrying capacity (McLeod 1997). We note that there are approaches that accommodate time-varying carrying capacity (e.g. Hone et al. 2007) and Gompertz-style density-dependent population growth. These are not explored here.

Sampling model: The proportion of the superpopulation on surveyed wetlands within bands (\(\phi _t\)) is specified by a prior for population coverage (\(\pi _c\)) based on the sampling fraction, combined with a linear effect of rainfall anomaly in the antecedent 6 months, and a random effect (\(\tau \)) for year, with the effects modelled on a logit scale:

$$\begin{aligned} \text {logit}(\phi _t)= & {} \text {logit}(\pi _c) + \beta _2 x_{2t} + \tau ,~~\text {where} \nonumber \\ \pi _c\sim & {} \text {Beta}(a,b),~~\text {and} \nonumber \\ \tau\sim & {} N(0,\sigma _{\tau }^2). \end{aligned}$$
(2)

Observation model: Total counts of waterbirds on surveyed wetlands within bands are a proportion (\(\phi _t\)) of the superpopulation, with log-normal observation error:

$$\begin{aligned} \log (n_t)= & {} \log (\phi _t N_t) + V,~~\text {where} \nonumber \\ V\sim & {} N(0,\sigma _V^2) \end{aligned}$$
(3)

Zero counts were replaced with a count of one, affecting 14 of the 45 species, typically in a minor way. For 5 species (gull-billed tern Gelochelidon macrotarsa, black-necked stork Ephippiorhynchus asiaticus, musk duck Biziura lobata, pied cormorant Phalacrocorax varius, purple swamphen Porphyrio porphyrio) this substitution was only needed for a single year (median non-zero counts 260, 9, 30, 1149, 163 respectively). For three species (chestnut teal Anas castanea, magpie goose (Anseranas semipalmata), plumed whistling duck Dendrocygna eytoni) the substitution was made in 2 years (median non-zero count 229, 2762, 1033). For two species (banded stilt Cladorhynchus leucocephalus, great crested grebe Podiceps cristatus) it was 4 years (median non-zero counts 955 and 31 respectively). For 1 species (dusky moorhen Gallinula tenebrosa) it was 8 years (median non-zero counts 20), and for 2 species (Cape Barren goose (Cereopsis novaehollandiae), rajah shelduck (Radjah radjah)) is was 9 years (median non-zero counts 74 and 22). In the case of the banded lapwing (Vanellus tricolour) there were zero counts on 11 years (median non-zero count 17). Given the low percentage of zero counts, and the typical count being an order of magnitude greater than the substitution value, we were confident of not having introduced significant bias.

Finally, implicit in this model formulation (Eq. 3) is that posterior samples of the superpopulation \(\widehat{N_t}\) can be estimated using the derived quantity:

$$\begin{aligned} \widehat{N_t}=\frac{n_t}{\widehat{\phi _t}}. \end{aligned}$$
(4)

This is essentially a model-based implementation of the Ratio Method (Caughley 1977), where the sampling fraction is treated as a parameter coming from a distribution (c.f. a fixed scalar) within an observation model (see Borchers et al. 2002, p.44). Also implicit in this estimator is that the probability of detection within a surveyed wetland is certain (i.e. \(p=1\)). This is revisited in Sect. 4.

2.3 Model priors

The priors chosen for model parameters along with their rationale are given in Table 1. For priors on variances, following Gelman (2006) we used a uniform prior on the scale of the standard deviation. The truncated normal prior for \(r_0\) (rate of increase under average antecedent rainfall) may initially appear highly informative, but we note it provides considerable support well beyond the range (– 0.18 to 0.29 year\(^{-1}\)) of empirical observation of long-term rates of increase collated by Saether and Engen (2002).

Table 1 Details of the parameters and prior values of the statistical model fitted to the 45 waterbird species
Full size table

2.4 Model fitting

Models were fitted separately for each species using the “rjags” library (Plummer et al. 2016) for the software R (R Core Team 2018). To account for possible effects of prior conditions on the starting population size, model runs were all started from 1973 (a 10-year lead-in). Each model run consisted of 3 Monte Carlo Markov chains, with a posterior of 100,000 samples drawn following a burn-in period of 20,000 samples and an adaptive phase of 1000 samples. Convergence of chains was assessed using the convergence diagnostic of Gelman and Rubin (1992) that calculates potential scale reduction factors.

The full posterior to be sampled is:

$$\begin{aligned}&[N_t,\beta _1,\beta _2,\phi _t,r_0,\pi _c,\sigma _\tau ^2,\sigma _W^2,\sigma _V^2 |n_t,N_{t-1},x_{1t},x_{2t},a,b]\nonumber \\&\quad \propto [n_t|N_t,\phi _t,\sigma _V^2] \nonumber \\&\qquad \times \; [N_t|N_{t-1},r_t] \times [r_t|r_0,\beta _1,x_{1t},\sigma _W^2] \nonumber \\&\qquad \times \; [\phi _t|\pi _c,\beta _2,x_{2t},\sigma _\tau ^2] \times [\pi _c|a,b] \nonumber \\&\qquad \times \; [r_0][\beta _1][\beta _2][\sigma _V^2][\sigma _W^2][\sigma _\tau ^2] \end{aligned}$$
(5)

Posterior distributions for parameters were summarized by the mean (e.g. \(\bar{r}_0\)) and credibility intervals (C.I.). For posterior plots of the coefficients \(\beta _1\), \(\beta _2\) and \(r_0\) we present 90% credibility intervals, and hence a credibility interval excluding zero implies c. \(>95\)% belief the parameter differs from zero (assuming symmetric posterior distributions).

2.5 Model evaluation

We evaluated the model within the context of its purpose to estimate the population growth rate under average conditions (\(r_0\)). To do this we compared the width of the credibility intervals around \(r_0\) for the full model (Eq. 5) with those from a reduced model that had neither a model-based antecedent rainfall effect, nor yearly random effect for the proportion of the superpopulation available to be counted, while retaining the effect of rainfall on rate of increase (Eq. 6). Given the importance of rainfall to wetland availability, comparisons with a model without rainfall-driven demographics would not be useful.

$$\begin{aligned}{}[N_t,\beta _1,\phi _t,r_0,\pi _c,\sigma _W^2,\sigma _V^2 |n_t,N_{t-1},x_{1t},a,b]\propto & {} [n_t|N_t,\phi _t,\sigma _V^2] \nonumber \\&\times \; {[}N_t{|}N_{t{-}1}{,}r_t{]} \times {[}r_t{|}r_0{,}\beta _1{,}x_{1t}{,}\sigma _W^2{]} \nonumber \\&\times \;[\phi _t|\pi _c] \times [\pi _c|a,b] \nonumber \\&\times \;[r_0][\beta _1][\sigma _V^2][\sigma _W^2] \end{aligned}$$
(6)

As \(\beta _1\) and \(\beta _2\) are both functions of measures of antecedent rainfall that are likely correlated, it could be questioned whether \(\beta _1\) and \(\beta _2\) are identifiable. They are. Briefly, a year-on-year change in population when antecedent conditions are the same will be ascribed to \(r_t\) (and hence informing \(\beta _1\) as in Eq. 1), as the proportion available to be counted will be modelled as remaining the same (as per Eq. 2), and vice versa. In addition, the correlation between \(\beta _1\) and \(\beta _2\) is surprisingly low (\(r=0.46\), see Fig. 1 and Fig. 2 in Supplementary Material). We further checked for the identifiability of \(\beta _1\) and \(\beta _2\) by examining joint posterior plots for the presence of ridges.

3 Results

Graphical summaries of all effects for the 45 species and plots of inferred population trajectories for a limited number of species are presented in Figs. 3, 4, 5, 6, 7, 8, 9, 10. Full summary details of fitted models for all species, including 95% C.I.s for all parameters, plots of counts, inferred proportions on the surveyed stratum, and inferred population trajectories, are contained in the Supplementary Materials (Figs. 3–47).

Fig. 3
figure 3

Estimated coefficient (\(\beta _2\)) for the effect of antecedent 6-month rainfall anomaly from the 50-year mean on the proportion of waterbird populations that were on the wetland stratum within bands sampled by the East Australian Waterbird Surveys. Intervals are 90% credibility intervals

Full size image

3.1 Model evaluation

For 44 of the 45 species assessed, adequate convergence (Gelman & Rubin’s \(R<1.1\)) was reached for the key parameters of interest (\(r_0\), \(\beta _1\) and \(\beta _2\)), the exception being the white-faced heron (Egretta novaehollandiae). Joint posterior plots confirmed that \(\beta _1\) and \(\beta _2\) were uniquely identifiable for all species (see Supplementary Material Fig. 48), and that the prior values for \(r_0\) were strongly informed by the data (see Fig. 49 in Supplementary Material). The full model performed better in estimating \(r_0\) for 35 of the 45 species assessed, with 95% credibility intervals reduced on average by 25% (range [\(-14\%, 68\)%]), compared with the model without the sampling model.

3.2 Proportion of superpopulation on surveyed wetlands

Effect of rainfall: The effect of rainfall on the proportion of waterbirds on surveyed wetlands ranged from strongly negative (e.g. freckled duck Stictonetta naevosa, pink-eared duck Malacorhynchus membranaceus) to strongly positive (e.g. gull-billed tern) (Fig. 3). At the functional group level, ducks tended to have strongly negative values of \(\beta _2\), providing model-based evidence for their dispersal across the landscape (beyond the sampled wetland stratum) in response to recent rain, as did the pelicans, darters and cormorants group (Fig. 3). Large waders exhibited a mixed response, with some of the colonial nesting species (e.g. ibis and spoonbill species) showing large positive effects of rainfall on the proportion of the superpopulation on wetlands within surveyed bands, though others (e.g. glossy ibis) having a strong negative response (Fig. 3).

Unexplained variation: The standard deviation (\(\sigma _\tau \)) for the random effect (\(\tau \)) in Eq. 2 can be considered as a proxy for the unexplained year-to-year variation in the proportion of the superpopulation present on surveyed wetlands. This varied considerably, generally being lowest for the piscivores (e.g. Pelicans, darters and cormorants group) and some of the aquatic herbivores (e.g. black swan), followed by large waders, ducks and small waders (Fig. 4). Two species with high unexplained variation—the banded lapwing and gull-billed tern—forage extensively in dryland environments outside of the surveyed stratum (e.g. Marchant and Higgins 1990, 1993). Species of conservation interest that have a reasonably high unexplained variance include the glossy ibis and freckled duck (Fig. 4).

Fig. 4
figure 4

Estimated standard deviation (\(\sigma _\tau \)) for the \(N(0,\sigma _\tau )\) additive random effect for the proportion (on a logit scale) of the total waterbird population that were on the wetland stratum within bands sampled by the East Australian Waterbird Surveys. Intervals are 95% credibility intervals

Full size image

Overall effect: The estimated year-to-year variation in \(\phi \) is a result of the combined effects of the antecedent rainfall term and the yearly random effect. For some species, such as the black swan, both these effects are small and variation in \(\phi \) from year-to-year is correspondingly small (Fig. 5a). For the glossy ibis, the combination of a strong effect of rainfall-driven absence from the surveyed stratum and unexplained variation leads to moderately high year-to-year variation in \(\phi \) (Fig. 5b). Finally, for the banded stilt, despite there being no clear effect of rainfall on \(\phi \), the estimated year-to-year variation in \(\phi \) is very large due to unexplained variation, including runs of low years (Fig. 5c). We return to this in Sect. 4.

Fig. 5
figure 5

Estimated proportion (\({\hat{\phi }}\)) of the waterbird superpopulation on the wetland stratum within bands sampled by the East Australian Waterbird Survey for a black swan (Cygnus atratus), b glossy ibis (Plegadis falcinellus) and c banded stilt (Cladorhynchus leucocephalus). Note the scale of the y-axes differs between plots

Full size image

3.3 Rates of increase

Rainfall response: The effect of antecedent rainfall on the rate of population increase was estimated to be positive with strong belief (i.e. \(P(\beta _1>0)>0.95\)) for about half the species examined. Of the remaining species, all but two had a positive estimate for \(\beta _1\), though with lesser belief of a positive effect (Fig. 6). The estimated effect sizes are non-trivial. For example, the La Niña conditions of 2010 and 2011 delivered rainfall c. 200 mm above average across the survey region for the October–September period. For a value of \(\beta _1=0.3\) (the typical response for about one third of the species), this equates to an increase in r of c. 0.6 \(\text {yr}^{-1}\) compared with average conditions (\(r_0\)). In linear terms, this corresponds to a finite growth rate of \(\lambda \simeq 2~\text {yr}^{-1}\). The effects of dry conditions are of similar magnitude though in the other direction. For example, the drought conditions of 2018, with a negative anomaly of c. 120 mm compared to the long-term mean, equates to a finite growth rate of \(\lambda \simeq 0.7~\text {yr}^{-1}\) for \(\beta _1=0.3\).

Fig. 6
figure 6

Estimated coefficient (\(\beta _1\)) for the effect of antecedent 12-month rainfall anomaly from the 50-year mean on the rate of increase of waterbird populations sampled by the East Australian Waterbird Surveys. Intervals are 90% credibility intervals

Full size image

Estimates under average rainfall: Under average antecedent rainfall conditions, the rates of increase were negative for most of the species assessed, though for no species was the belief in an ongoing decline (i.e. \(r_0<0\)) greater than 95% (Fig. 7). That said, there were several species for which the estimated rate of increase was substantially negative though imprecise, including the freckled duck, Australasian shoveler (Anas rhynchotis), whiskered tern (Chlidonias hybrida), pink-eared duck and glossy ibis (Fig. 7). The species with the strongest belief of decline under average conditions was the black-winged stilt (Himantopus himantopus) (\(\bar{r}_0=-0.06\) year\(^{-1}\), \(P(r_0<0)=0.92\)) followed by the Australasian shoveler (\(\bar{r}_0 = -0.09\) yr\(^{-1}\), \(P(r_0<0)\) = 0.90). The estimates of \(r_0\) for ducks tended to be at the lower end of the scale (Fig. 7). For some functional groups, such as the piscivorous pelicans, darters and cormorants, the inferred population trajectories suggest fluctuating though stable populations (see Fig. 7, 8 and Supplementary Material). Finally, there was a single species, the Cape Barren goose (Cereopsis novaehollandiae), for which there was strong belief in an increasing population trend under average rainfall conditions (\(\bar{r}_0=0.14\) yr\(^{-1}\), \(P(r_0>0)\) = 0.97).

Fig. 7
figure 7

Estimated median exponential rate of increase (\(r_0\)) for waterbird populations sampled by the East Australian Waterbird Surveys subject to average antecedent rainfall conditions (baseline 1969–2018). Intervals are 90% credibility intervals

Full size image
Fig. 8
figure 8

For the Australian pelican (Pelecanus conspicillatus), a inferred median superpopulation size (solid line) with 95% credibility intervals (grey shading) estimated from b summed counts on wetlands surveyed within the East Australian Waterbird surveys (blue bars)

Full size image

4 Discussion

We have presented a new model-based approach to making inference on population size and trends for 45 waterbird species commonly observed by the EAWS. For most species, the model estimates substantial effects of antecedent rainfall on population growth rate and/or the proportion of the population available to be counted on wetlands contained within survey transects, with random effects additionally used to account for sampling variation. A striking result across the 45 species assessed was the general lack of strong belief in either a decreasing or increasing long-term population trend. This result is somewhat at odds with previous publications that have inferred major declines in waterbird species populations in eastern Australia using either all the EAWS data (Kingsford and Porter 2009), or the Murray-Darling Basin subset of the data (e.g. Kingsford et al. 2017; Bino et al. 2020). We explore reasons for these differing conclusions after first discussing some aspects of our modelling approach and potential improvements in terms of additional covariates and model specifications. We then discuss some of the design issues inherent in the EAWS, before concluding with fruitful lines of future inquiry.

4.1 Model specification

Our decision to ignore band-to-band variation invites criticism as wasting information. We have taken the view that these data are more likely to mislead than inform, typically by underestimating variation when counts are sparse. The greater the spatial aggregation, the higher the chances that yearly estimates of band-to-band variation will be biased downward: the information such variation provides is not currently useful. For species with sparse counts, a higher sampling fraction (which is logistically challenging) would be needed to reduce bias, but in the current context will not alleviate problems arising from waterbirds using non-surveyed strata at the time of survey or moving outside the survey region. Developing methods for using band-to-band variation to help inform the uncertainty around the proportion of waterbirds available to be counted would appear a useful avenue for further research. In the meantime we have focussed on exploring how the proportion of the population that was not available to be counted can be estimated. Although we have presented one approach, we note that we have not explicitly accounted for the inter-annual variability in extent and productivity of wetlands across the whole of Australia and southern Papua New Guinea.

4.2 On explaining variation

We have taken a one-model-fits-all approach (i.e. with the same structure and covariates) for the 45 waterbird species considered here. Even within this simple model structure, it is illuminating how species varied in their utilisation of the surveyed wetland stratum. This finding agrees with Roshier et al. (2002) who found changes in abundance of most functional groups of waterbirds in north-western New South Wales changed in response to changes in wetland distribution, and that some groups responded immediately to wetland distribution and area flooded (e.g. dabbling ducks) whereas others responded to sequences of wetting and drying, both within the focal region and neighbouring river basins. Models tailored for individual species would consider different covariates (possibly including specific threats and pressures), and different spatial scales. This is an area for future research.

Our model does not currently distinguish between uncounted individuals being outside the surveyed stratum though within the survey area, and those outside the survey area altogether. The high unexplained variances for the radjah shelduck, magpie goose and plumed whistling duck are likely caused by the latter, as their typical range straddles the northern boundary of EAWS survey area, with greater movement into the survey area on occasion in response to seasonal conditions. Likewise, recent studies suggest that the entire population of banded stilt can move into and out of the survey region (Pedler et al. 2014, 2018), which explains the high year-to-year variation in \(\phi \), including very low (and sometimes zero) counts in some years (particularly during the Millennium Drought: see Fig. 5c and Fig. 37 in Supplementary Material).

4.3 On inferring trend

Our approach for analysing trend accounts for serial dependence in the counts through the population being sampled, and the effect of time-varying rainfall covariates on population demographics and movements. Rainfall will not necessarily be stationary over the period of the observational experiment. Indeed, our results for a number of species confirm that the above average rainfall during the 1970s could well explain the exceptionally high counts recorded during the early years of the EAWS. The grey teal is one such example, whereby the estimated superpopulation size in the 1970s is higher than at any other time period (Fig. 9). The structure of our model enables estimation of the rate of population increase under average conditions. This is an improvement on simply estimating population trend effectively marginalised over the environmental conditions experienced during the observation period. The resulting estimate, however, is influenced by the time-window over which the average is calculated. The baseline period we have presented here (1969–2018) includes the run of above average rainfall years in the 1970s. If, for example, only the EAWS survey years are included (1983–2018) to calculate baseline conditions, the deviations below the estimated mean are less, and the inferred rate of increase under average conditions (\(r_0\)) is lower. An exploration of what the future average rainfall may be is beyond the scope of this work.

Fig. 9
figure 9

For the grey teal (Anas gracilis), a inferred median superpopulation size (solid line) with 95% credibility intervals (grey shading). Note the logarithmic (base 10) scale for the population estimates. b Summed counts on wetlands surveyed within the East Australian Waterbird surveys (blue bars). Note that the first survey counts start from 1983

Full size image

An early analysis of trends in the EAWS for a selection of speciesFootnote 2 by Kingsford and Porter (2009) ignored serial dependence in counts when estimating exponential rates of increase from log-linear regression, resulting in the precision of estimates being overestimated (with increased risk of Type I errors). This problem is well known (Humbert et al. 2009), with better methods readily available (e.g. Holmes et al. 2012), as has been previously pointed out when analysing East Australian Waterbird Survey counts (Colloff et al. 2015). A subsequent analysis of Kingsford et al. (2017) using the EAWS data to compare rates of increase between Lake Eyre Basin and the Murray-Darling Basin again ignored such serial dependence in counts; both within log-linear regression, and additionally in linear models comparing mean counts between the first and last decades of the EAWS. The risk of not accounting for serial dependence is that the strength of inference (as typically expressed by p-values in relation to null hypothesis tests of significance) is overstated.

In our modelling efforts, the effect of accounting for rainfall, observation error, and serial dependence in the data when making inferences around population trend, has been to greatly reduce the certainty around what appeared to be compelling declines in some species. For example, visualising the trajectory of counts for the Australasian shoveler, and inferred superpopulation size (Fig. 10), would leave little doubt in most minds that the population has decreased “significantly”. Indeed, comparing the posterior distribution for the initial and final inferred population size supports this (Fig. 10b and Fig. 10c). Our model however, after correcting for rainfall effects is more reserved as to whether there is a long-term downtrend underway, with \(P(r_0<0)\) = 0.90. It could be reasonably argued from a management perspective that a 90% belief in a long-term decline provides strong grounds for taking action to reverse the trend (e.g. Field et al. 2004). The Australian white ibis is another species that appears to have undergone a dramatic decline based on visually inspecting the counts (see Fig. 7 and Fig. 24 in Supplementary Material), yet the 90% credibility intervals still contain zero (\(P(r_0<0)=0.88\)). In this case, however, any concern that this is a true decline is tempered by the observation that the ecology of this species has shifted from natural wetlands towards urban habitat over the period of this study (Martin et al. 2010; Smith et al. 2013). These urban habitats are not part of the surveyed stratum.

Fig. 10
figure 10

For the Australasian shoveler (Anas rhynchotis), a inferred median superpopulation size (solid line) with 95% credibility intervals (grey shading), and b inferred median superpopulation size (solid line) and 95% credibility intervals (grey shading) on a logarithmic (base 10) scale, and c summed counts on wetlands surveyed within the East Australian Waterbird surveys (blue bars)

Full size image

Regarding the single species with strong belief to be increasing in numbers, the Cape Barren goose, this increase is explicable on the basis of the population recovering following severe declines by the 1930s that resulted from persecution for grazing crops and fouling pastures, and use as a food source. Under protective measures, the population began to recover in early 1960s (Parker et al. 1985).

There is a difference between an ongoing, genuine population decline (i.e. negative population growth rate being driven by an underlying negative rate of increase) and a population decrease arising from fluctuations in time-varying covariates. The former will eventually lead to a small and endangered population unless the cause of the decline can be identified and mitigated (Caughley and Gunn 1996), whereas the latter may not be of concern. The inference from our results is that there is no strong belief for a declining or increasing trend in any of the 45 species modelled, after correcting for rainfall and if adopting the (arbitrary) 95% credibility interval as “strong evidence”. It appears likely that a component of the previously inferred decline in the abundance of waterbirds over the period of the EAWS is driven by high starting populations arising from a run of wet years during the 1970s. Note that any effect of changed water extraction (e.g. for irrigated agriculture) on wetland availability is subsumed within our estimation of \(r_0\). Increasing rates of water extraction for irrigation over the survey period has been argued (e.g. Kingsford et al. 2017; Bino et al. 2020) to be impacting on the number of waterbirds utilising the Murray-Darling Basin. A concomitant decline in the superpopulation would be expected given that the Murray-Darling Basin is a substantial component of the range for many species. However, our analysis suggests that such declines in the superpopulation are less certain. Some possible reasons for this discrepancy include that the analysis of Bino et al. (2020) did not account for the effect of rainfall in driving fluctuations in waterbird abundance, and that data were pooled across species within functional groups prior to analysis (thus smoothing over much of the sampling variation seen at a species level). Alternatively, and despite our modelling efforts, the high unexplained year-to-year sampling variation is undoubtedly impacting on the ability of the dataset to be used to detect trends for some species.

4.4 Accuracy of population estimates

Aerial surveys of waterbirds (and most other animal taxa) usually underestimate population size, as measured from ground counts (e.g. Kingsford 1999) or aerial photography (e.g. Bayliss and Yeomans 1990). With regard to estimating the total population, the bias associated with the EAWS counting methods has not been calibrated, so we are reluctant to incorporate such a correction for detection probability (e.g. by putting a prior on p). The expectation, however, would be that the estimates are underestimates. Aerial counts of waterbirds have also been found to vary by at least an order of magnitude between different observers and different species (Morton et al. 1993). Such imprecision may limit the spatial resolution at which these data can be analysed, though the sums of the counts across all wetlands and bands (as we have done here) should have considerably better precision. If the accuracy of the population size is of importance, then future surveys should include a calibration component (e.g. Morton et al. 1993), such as taking ground counts at a modest number of discrete wetlands along the aerial survey bands.

Australian pelicans represent one of the better species for assessing how accurate the population estimates are, because they are large, conspicuous, piscivorous (and hence an obligate wetland species), relatively slow-moving and typically do not occur in very large (and possibly mixed) flocks, other than in breeding aggregations (see below), that makes counting extremely challenging. The population of pelicans in Australia is thought to number around 300,000–500,000 (Jaench unpubl., quoted by Reid 2009). The estimate for the population sampled by the EAWS ranges over c. 50,000–250,000 (Fig. 8). There are few independent data to test the accuracy of the population estimates. In April 2002, there were an estimated minimum 80,000 pelicans nesting in a single colony on Lake Machattie (south-west Queensland), within the survey region following flooding in the Lake Eyre Basin (Reid 2009). In comparison, the estimated superpopulation for late 2002 was c. 54,000 (95% C.I. 24,000–121,000).

4.5 Issues with the design of the EAWS

A feature of the data arising from the EAWS is the very high year-to-year variation in counts, of which much is caused by sampling variation. Variation in the proportion of the total population within survey bands could be driven either by rarity, spatial aggregation (even if species are not rare), and the use by waterbirds of non-surveyed strata. There are other, largely irreducible sources of variation such as waterbirds being in transit between wetlands at time of survey, or moving between wetlands during the period (c. 4 weeks) of the survey. This leads to high sampling within- and between-year variation within the surveyed stratum. For some species such as the larger piscivores (e.g. Australian pelican), and herbivores (e.g. black swan), our analyses suggest that the current wetland stratification is probably adequate, in that the year-to-year random effect for the proportion of these populations on the surveyed stratum is small. However, a substantial proportion of waterbird species spend considerable periods outside of the surveyed wetland stratum, leading to additional sampling variation (Fig. 4). Banded stilt, black-tailed native-hen (Tribonyx ventralis) and plumed whistling duck are good examples, being uncommon in some years, subject to clumping behaviour and, in the latter two cases, disposed to terrestrial foraging. For such problematic species, we could simply choose to accept this variation, and continue to gather additional data and make inferences on trend regardless. Indeed, a common response to lack of power in statistical studies arising from sampling variation is to increase the number of samples, and/or in the case of longitudinal monitoring, to monitor for longer (Gerrodette 1987). Increasing the number of survey bands within a given year will certainly help (though would be logistically difficult), but would not solve the problem of species using non-wetland habitat at the time of survey. For many species, this appears to generate such high year-to-year variation that no amount of additional years of surveying will improve the statistical power sufficiently to detect trends. Accounting for sampling variation appears to be a key research need for the statistical power of analyses to be improved.

A possible alternative solution to the movement in and out of the surveyed stratum (a “temporary emigration” problem, see Chandler et al. 2011), and one that we have explored here, is to explicitly model yearly changes in the proportion of the population that is available to be counted within the observation model. The covariate we have used for this (antecedent 6-month rainfall), is a proxy measurement of wetland area and resource availability, and could be improved on. More generally, improvement necessitates making inferences for the considerable area (c. 89% of wetlands) that is not surveyed, both in terms of resource availability and waterbird abundance. Relevant environmental layers are increasingly becoming available that can be used to help model this unobserved component. One useful example would be wetland area availability, at an appropriate spatial and temporal resolution to fit with the EAWS survey design. Such covariates may add considerably to our understanding of waterbird population dynamics at continental scale. Other ancillary data sources of waterbird abundance such as bird lists could also be of assistance (e.g. Wen et al. 2016). A further complication arises due to the change in breeding season of many waterbirds, largely spring breeding in southern Australia to late summer-autumnal breeding in the north (see Marchant and Higgins 1990, 1993).

Extra information gained from additional surveys of waterbird populations outside the current survey coverage (though still within the designated wetland stratum) will be most easily captured and incorporated into inference if these areas are surveyed using the same methods as the EAWS. Band rotation would be one way of achieving this, without greatly altering the cost of the current surveys. For example, rotating 20% of bands each year would mean each band would be surveyed consecutively for five years, with all the survey area covered over c. 50 years. To continue with the status quo and survey the same bands indefinitely would miss the opportunity to reveal what might be occurring in the remaining 90% of the study area, thus providing a means of improving inference. How to estimate waterbird numbers not on the wetland stratum is an open question. The overall density of waterbirds outside of the current wetland stratum is expected to be low (despite the numbers involved potentially being high), so extending the aerial surveying to all these non-wetland areas will probably be inefficient. An alternative approach is to improve on the current model-based estimate of \(\phi \). Independent data on \(\phi \) could be obtained from habitat use studies, with selection of the existing wetland stratum modelled as a function of measurable covariates. Incorporation of such information into the existing Bayesian model would be straightforward, and the animal tracking technology to collect such data already exists (e.g. Roshier and Asmus 2009).

As would be expected with a sample-based survey approach, there are major wetlands not included in the EAWS survey bands. Having large numbers of individuals not on survey bands is not an issue if this phenomenon is occurring at random, however this may not be the case for aggregations of waterbird associated with breeding. These sometimes dense aggregations tend to recur at specific locations. Knowledge of the locations of such breeding locations has built since the original design of the EAWS. For example, in the Murray-Darling Basin, several Ramsar wetlands listed for their importance as waterbird breeding sites are not included, including Narran Lakes, Gwydir Wetlands, Hattah Lakes and Barmah-Millewa Forest. In the Lake Eyre Basin, major wetlands in the Georgina, Diamantina and Cooper Creek catchments are excluded, including all of the Coongie Lakes. In 2001 (a La Niña year), Reid and Jaensch (2004) estimated there to be c. 4 million waterbirds in Lake Eyre Basin wetlands that were not part of the EAWS survey bands. Such breeding locations could qualify as an additional survey stratum if the survey was to be redesigned. Of course, the large number of species involved makes identifying useful strata (e.g. breeding sites) that cover off on all species difficult. Alternatively, incorporating band rotation using the existing survey protocol would be a means to start learning about waterbird populations more widely in a structured way, without precluding comparisons with previous surveys.

4.6 Conclusions

Our analysis is the first to use the EAWS data to provide yearly estimates for the total population of individual waterbird species being sampled; the main purpose for which it was designed. To do this we developed a model-based framework to better explain variation in the observation and state processes using readily available rainfall data. This has enabled us to make inferences on long-term rates of increase that avoid issues of (1) serial dependence in the data, (2) sampling variability arising from waterbirds not exclusively being confined to the surveyed stratum, and (3) fluctuations in population size driven by large fluctuations in resource availability over the course of the study. We expect that generating better covariates (e.g. annual fine-scale mapping of wetland area across Australia and southern Papua New Guinea) to model waterbird habitat selection and how this relates to the surveyed stratum can lead to further model improvements, and hence better inferences from this important dataset.