Introduction

The Florida coast is one of the most species-rich ecosystems in the world. This paper focuses on the sensitivity of the habitat of threatened and endangered shorebirds to sea level rise induced by climate change, and on the relationship of the habitat with the coastline evolution. We consider the resident Snowy Plover (Charadrius alexandrinus nivosus), and the migrant Piping Plover (Charadrius melodus) and Red Knot (Calidris canutus) along the Gulf Coast of Mexico in Florida.

Methods

We analyze and model the coupled dynamics of habitat patches of these imperiled shorebirds and of the shoreline geomorphology dictated by land cover change with consideration of the coastal wetlands. The land cover is modeled from 2006 to 2100 as a function of the A1B sea level rise scenario rescaled to 2 m. Using a maximum-entropy habitat suitability model and a set of macroecological criteria we delineate breeding and wintering patches for each year simulated.

Results

Evidence of coupled ecogeomorphological dynamics was found by considering the fractal dimension of shorebird occurrence patterns and of the coastline. A scaling relationship between the fractal dimensions of the species patches and of the coastline was detected. The predicted power law of the patch size emerged from scale-free habitat patterns and was validated against 9 years of observations. We predict an overall 16% loss of the coastal landforms from inundation. Despite the changes in the coastline that cause habitat loss, fragmentation, and variations of patch connectivity, shorebirds self-organize by preserving a power-law distribution of the patch size in time. Yet, the probability of finding large patches is predicted to be smaller in 2100 than in 2006. The Piping Plover showed the highest fluctuation in the patch fractal dimension; thus, it is the species at greatest risk of decline.

Conclusions

We propose a parsimonious modeling framework to capture macroscale ecogeomorphological patterns of coastal ecosystems. Our results suggest the potential use of the fractal dimension of a coastline as a fingerprint of climatic change effects on shoreline-dependent species. Thus, the fractal dimension is a potential metric to aid decision-makers in conservation interventions of species subjected to sea level rise or other anthropic stressors that affect their coastline habitat.

Florida coastline-dependent species are characterized by one of the highest extirpation risks in the world because of sea level rise and increase in tropical cyclone activity (Convertino et al. 2010;2011c) due to climate change. The Snowy Plover (Charadrius alexandrinus nivosus; SNPL hereafter) is a residential shorebird of Florida listed as threatened at the state level. The Piping Plover (Charadrius melodus; PIPL hereafter) is federally designated as threatened, and it migrates mostly from the North Atlantic coasts of the USA and Canada to Florida where it winters for 3 months on average (Elliott Smith and Haig 2004). The Red Knot (Calidris canutus; REKN hereafter) is designated as threatened in New Jersey and is federally listed as a potential “at risk” species. REKN uses the Florida Gulf beaches as stop-over areas for about 3 weeks during its migration between South America and North America’s Big Lakes region and Atlantic coast (Harrington 2001). This is considered as the wintering period of the REKN in Florida. An understanding of the spatial distribution of the suitable habitat patches for these shorebirds, their controlling factors, and how these factors are affected by sea level rise is fundamentally important for adopting efficient conservation strategies. An understanding of linkages between the coupled evolution of landforms and ecological patterns is a crucial topic due to the evidence that these patterns are tightly linked. Biocomplexity approaches (Mandelbrot 1982; Rinaldo et al. 1995; Banavar et al. 2001; Pascual et al. 2002; Schneider and Tella 2002; Buldyrev et al. 2003; del Barrio et al. 2006; Solé and Bascompte 2006; Scanlon et al. 2007), despite being accused of adopting simplified biological models (Paola and Leeder 2011), are capable of reproducing macroscale patterns of complex phenomena and of developing indicators, such as the probability of the patch size (Mandelbrot 1982; Bonabeau et al. 1999; Jovani and Tella 2007; Kéfi et al. 2007; Jovani et al. 2008; Convertino et al. 2012), that are useful for assessing ecosystem health (Kefi et al. 2011). One of ecology’s main goals is to detect from observed patterns, such as species occurrence patterns, the organizational rules of species in stationary and evolving ecosystems. Many theories have been proposed to explain the formation of clustered patterns of species in nature. Conspecific attraction, environmental heterogeneities, and food availability have been claimed—alone or together—to be the motivation for the formation of habitat patches in which individuals of a species coexist in colonies. An optimal search theory, the so-called Lévy-flight foraging hypothesis (or predator-prey-food resource dynamics), predicts that predators should adopt search strategies known as Lévy flights where prey is sparse and distributed unpredictably. However, Humphries et al. (2010) showed that Brownian movement is sufficiently efficient for locating abundant prey. This theory explains the clustered patterns of resources in landscapes that may be different from the pattern of species occurrence. Neither the Lévy-flight foraging hypothesis nor Brownian movement model address the linkages of biota with landscape forms and their evolution, which is, in our opinion, one of the main missing points.

The colony size of seabirds (Schneider and Tella 2002), colonial birds (Jovani and Tella 2007), and many other other animals (Bonabeau et al. 1999) has been found to follow a power-law distribution. Analogous scale-free distributions have been detected for bacteria colonies (Buldyrev et al. 2003), for species in complex ecosystems (Solé and Bascompte 2006; Convertino et al. 2012), and also for man-made systems such as cities (Batty and Longley 1994). The ubiquity of the power-law structure for the probability of the patch size in aggregation phenomena of natural and human systems suggests the existence of universal self-organization principles (Pascual et al. 2002; Solé and Bascompte 2006). The scaling exponent of the power-law distribution of the aggregate size was proven to be the fractal dimension of the pattern analyzed (Mandelbrot 1982; Convertino et al. 2012). The word “aggregate” is a general word for indicating the assemblage of individuals with similar or identical features in a landscape. In the presence of a power law for the probability distribution of the aggregate size, the occurrence patterns are scale-free, indicating that the patterns are invariant at different scales of observations (Convertino et al. 2012). The concept of fractal dimension was introduced by Mandelbrot analyzing the coastline of Britain at different scales (Mandelbrot 1967). The work disseminated the use of fractal analysis first in geomorphology (Morais et al. 2011; Baldassarri et al. 2012) and later to a variety of sciences from biology to engineering (Bak 1999). Nonetheless, all these theories, models, and empirical findings have rarely considered any potential effect of slow or abrupt change in the exogenous factors on the heterogenous habitat in which species live. Only recently it was proven quantitatively that ecosystems exhibit variations in the probability distribution of the patch size due to anthropically and naturally driven changes in the environmental variables (Kefi et al. 2011). For example, desertification of water-controlled ecosystems produces a decrease in the fractal dimension of vegetation patches, or in extreme cases, a shift from the power law to exponential distribution of the patch size (Kéfi et al. 2007; Scanlon et al. 2007; Kefi et al. 2011). Climate change scenarios tested in temperate/continental regions depicted an overall decrease in the fractal dimension of patches in time for many different taxa (Barrio et al. 2006). For colonial birds the variation in the fractal dimension of the patches was clearly related to the fluctuations in the population abundance due to interspecies competition (Jovani et al. 2008). In geomorphology the variation in the fractal dimension was used as the signature of the persisting climate over landscapes. For example, the association between landscape evolution and climate has been assessed for river basin ecosystems in Rinaldo et al. (1995). However, none of the previous studies linked the fractal dimension of two ecosystems’ patterns in time (e.g., of geomorphological and ecological patterns) resulting from linked processes. Here we verify for the first time, to the best of our knowledge, that the fractality of the coastline is clearly linked to the habitat patches of shoreline-dependent birds in their breeding and wintering seasons.

We hypothesize that sea level rise may increase the complexity of the coastline and that such complexity determines fragmentation of the habitat of species. We assume scale invariance of the patches, which is also detectable by the analysis of the shorebird occurrences. We consider a breeding shorebird (Snowy Plover) and wintering shorebirds (Piping Plover and Red Knot) in Florida to quantify the potential effect of sea level rise on resident and migrant species. For the Snowy Plover the nesting season is usually considered part of the breeding season; thus, our model’s input considers the SNPL breeding and nesting occurrences simultaneously. Furthermore, observations indicate that nesting, breeding, and wintering areas for SNPL fall within the same range (Convertino et al. 2011a). Wintering occurrences of SNPL are thus considered together with breeding occurrences.

An integrated ecogeomorphological modeling approach is adopted to predict the viability from 2006 to 2100 of threatened, endangered, and at risk (TER) shorebirds (SNPL, PIPL, and REKN) along the Gulf Coast of Florida as a function of the increasing sea level rise due to climate change. We rescale to 2 m the Intergovernmental Panel on Climate Change (IPCC A1B) scenario described in Chu-Agor et al. (2011) and model the ecosystem at a 120 m spatial resolution. We predict land cover change with the Sea Level Affecting Marshes Model SLAMM (Clough 2010)] which is a geomorphological model at low-medium level of complexity. SLAMM considers coastal wetland types such as swamp, cypress swamp, mangrove, and salt marsh (Additional file 1: Figure S1). The habitat model predicts the habitat suitability for breeding and wintering through a maximum entropy principle approach (MaxEnt) (Phillips and Miroslav 2008) as a function of the recorded species occurrences in the breeding and wintering season, the predicted land cover, and a geology layer. MaxEnt is an ecological model at low level of complexity. The land cover and habitat simulations are produced in Aiello-Lammens et al. (2011). Finally, in this paper a patch-delineation model is introduced to predict the yearly habitat patches for a set of biological constraints imposed on the habitat suitability maps. We assume the stationarity of the habitat patterns at the year scale and absence of biological adaptation of species to climate change. The fractal dimension of the patches is derived by three independent methods: (i) box-counting for the observed occurrences; (ii) probability distribution of the patch size [“Korčak’s law” (Korcak 1940; Mandelbrot 1982)]; and (iii) perimeter-area relationship for the predicted patches. We assume that these three methods produce very close estimates of the fractal dimension of the whole mosaic of patches as shown in Convertino et al. (2012).

The power-law distribution of the patch size is verified by almost a decade (2002–2010) of historical observations of the species. Thus, the patch-delineation model is validated against these observations from 2002 to 2010. The coupled ecogeomorphological organization is shown by the correspondence in time of the fractal dimensions of the habitat-specific coastline and of the predicted patches. The fractal dimension of the habitat-specific coastline, along with habitat loss and population abundance, is demonstrated to greatly influence the number and size of the patches, which are related to habitat loss and population abundance. Although the fragmentation of the habitat (which is proportional to the fractal dimension of the patches) is predicted to fluctuate considerably in this century, the risk of extirpation of the species analyzed is not drastically increased because the connectivity of the patches is predicted to increase. The Piping Plover is the species with the largest fluctuation in the number and size of patches. We believe the research presented in this paper constitutes a contribution to the emerging field of biogeosciences, which explores the interface between biology and the geosciences and attempts to understand the interrelated functions of landscapes and biological systems across multiple spatial and temporal scales. We are aware of the existence of many other complex ecogeomorphological processes that are not included in our modeling effort. However, parsimonious models such as the model presented here can capture large-scale patterns while bypassing small-scale details (Ehrlich and Levin 2005; Pascual et al. 2011). These models can be tested against other more biologically realistic models to fully explore the linkages among various environmental changes, geomorphological dynamics, and biodiversity patterns. We anticipate that further research will explore this issue of process complexity versus model complexity, model relevance, and model uncertainty, which can be synthesized as a “modeling trilemma” (Muller et al. 2010).

This paper is organized as follows. The “Methods” section describes the shorebird data and the study site and explains the models used in this study and the theoretical characterization of patches. The “Results and discussion” section reports the main results with a broad discussion of figures and how these results are interpreted considering our assumptions. The “Conclusions” section reports the most important conclusions, implications for management, and further research efforts. Additional files 1 and 2 are provided to support our main result.

Site description and biogeographical variables

The white fine-sand beaches of the Florida coast of the Gulf of Mexico constitute the habitat of the whole Florida SNPL population. The SNPL population in Florida is distributed along about 80% of the Florida Panhandle and along about 20% of the Florida Peninsula (Lamonte and Douglass 2002; Himes et al. 2006; Burney 2009; Pruner 2010) (Figure 1a). The Florida Peninsula and the Atlantic coasts are the main wintering grounds for the migratory PIPL and REKN, which seem less constrained than the SNPL by the mineralogical properties of the beach substrate captured by the geology layer (Convertino et al. 2010;2011b). The land cover, which includes many wetland types from C-CAP (2009) is represented in Figure S1 of the Additional file 1, and the geology (F-DEP 2001) characterizes the mineralogical substrate of each land cover class (Additional file 1: Figure S6) (Convertino et al. 2011b). In 2006 the PIPL Panhandle-Peninsula and Atlantic populations were 38 and 33%, respectively, of the total migrant PIPL population in Florida. The REKN Panhandle-Peninsula and Atlantic populations were 55 and 20%, respectively, of the total migrant population in Florida. The International Piping Plover Census in 2006 supported the field sampling of SNPL, PIPL, and REKN (USGS-FWS 2009; FWC 2010; Alliance 2010). The 2006 wintering occurrences in Florida are the data used in this study for PIPL and REKN. For the SNPL, data of breeding and nesting occurrences are also available from 2002 to 2010 and are provided year by year by the Florida Wildlife Commission. These occurrences are used to verify the assumption of scale-invariance of SNPL occurrence patterns over time with the box-counting. However, despite the availability of SNPL data from 2002 to 2010, we construct the habitat suitability model with the 2006 SNPL occurrences alone in order to be consistent with the 2006 NOAA land cover (C-CAP 2009) and the 2006 PIPL and REKN occurrences. The geology and the elevation from USGS (USGS 2010; Convertino et al. 2010;2011b) are used in the habitat suitability model and in the land cover model, respectively (Aiello-Lammens et al. 2011; Convertino et al. 2010;2011b).

We consider PIPL and REKN in the same geographic domain where the full range of the SNPL occurs in order to perform a simultaneous interspecies assessment of the habitat use and extirpation risk of the three species (Figure 1a). Thus, only the Panhandle-Peninsula region was considered in this study. The SNPL is our main interest because its year-round presence in the Florida coastal ecosystem makes this species potentially more vulnerable than PIPL and REKN. Dispersal among the Panhandle and Peninsula SNPL populations has been observed but not quantified. Population subdivision of the SNPL has not been observed; thus, we can adopt the same habitat and dispersal criteria for the whole population. Population subdivision, for example, can be caused by geographic barriers or disturbances [e.g., renourishment (Convertino et al. 2011a)] that interfere with the dispersal. The reduction in dispersal is reported to reduce gene flow and increase genetic drift of independent subpopulations in the long-term. However, this is not the case for the SNPL population in Florida despite the weak interchange of individuals between Panhandle and Peninsula (Aiello-Lammens et al. 2011).

Habitat area and dispersal data for SNPL are mostly from Aiello-Lammens et al. (2011) but also from Page et al. (2009), Patons and Edwards (1996), Stenze et al. (1994; 2007), Warriner et al. (1986). Aiello-Lammens et al. 2011synthesized the biological data and the metapopulation modeling effort of this research for the SNPL. Information is gathered also from field ecologists working on this project [i.e., Dr. R.A. Fischer (Engineering Research and Development Center, US Army Corps of Engineers) and Mrs. A. Pruner (Florida Park Service)]. For PIPL, habitat and dispersal data are from Audubon (2006), Seavey et al. 2010, and USFWS (2009), and for REKN, data are from Fallon (2005) and Leyrer et al. (2006). For a more detailed description of the site under study we refer the reader to Convertino et al. (2011b).

Box-counting algorithm

The characterization of the occurrence patterns of breeding and wintering occurrences and of the coastline is performed using the “box-counting” method. For the SNPL the occurrence pattern of nesting occurrences was observed to be a self-similar pattern (Convertino et al. 2012); thus, the box-counting method is suitable to predict how this pattern changes with the scale of analysis. The box-counting analysis consists of calculating, for grids of different box-side lengths, the number of boxes that contain the object under study. Adjacent boxes constitute an approximation of the real patches at each resolution. The algorithm can be applied to both point and line patterns. The box-counting is performed over eight orders of magnitude in a logarithmic scale of the box-side length, from lo(1)=565km (which corresponds to the box “B1” in Figure 1), which is approximatively the width of Florida, to lo(5000)=11.3×10−3 km (Figure 1). We indicate with lo(i) the length of the box side at resolution “o=i” from i=1,…,5000 where the increment from one resolution to another is 11.3×10−2km, which is slightly smaller than the average home range of the SNPL (Table 1). The order of magnitude is relative to the scales of analysis (extent) investigated by the box-counting, while the resolution is related to the grids chosen for the box-counting. The number N(l)of boxes of size l needed to cover the pattern of occurrences (which is generally a fractal set) follows a power law,

N(l)=N0l−D,

(1)

where D≤d, and d is the dimension of the space (usually d=1,2,3). D is also known as the Minkowski-Bouligand dimension, Kolmogorov capacity, or Kolmogorov dimension, or simply the box-counting dimension and is an estimate of the Hausdorff dimension (Mandelbrot 1982). The fractal dimension for 1−d objects is associated with the Hurst exponent H such that D=2−H(Mandelbrot 1982; Bak 1999). The values of the Hurst exponent vary between 0 and 1, with higher values indicating a smoother trend, less volatility, and less roughness of the analyzed pattern (Mandelbrot 1982). We indicate the fractal dimension of the breeding and wintering occurrences with Db and the fractal dimension of the coastline with Dfderived from box-counting analysis. Both fractal dimensions are determined by the box-counting method. The fractal dimension of the coastline is calculated also for each land cover class that is a species-specific habitat for the species considered (Figure 1b). Many land cover classes are coastal wetland types (Additional file 1: Figure S1).

Figure 1

Box-counting algorithm. (a) Representation of the the box-counting algorithm applied to the 2006 occurrences of the Snowy Plover (SNPL), Piping Plover (PIPL), and Red Knot (REKN), for eight orders of magnitude (in a logarithmic scale), which corresponds to 5000 resolutions of the box-counting grid. In this example at the resolution of box B5 the number of boxes in which there is at least one occurrence is N(B 5)=6. (b) Box-counting example applied to the whole coastline, to the habitat-specific coastline (e.g., beach, salt-marsh), and to other land cover classes as in Convertino et al. (2011b). Many coastal wetland types are included in the land cover, such as swamp, cypress swamp, mangrove, and salt marsh. The shaded grid cells in (a) and (b) have at least one species occurrence or a coastline segment at the represented resolution. Two coastline configurations are presented: the first for high values of Dfand DK(c), the second for low values of Dfand DK(d). The patches presented in green are connected because their neighboring distance is lower than the maximum dispersal length dl.

Table 1

Macroecological parameters of the patch-delineation model, and biological data estimated from the literature

SNPL

PIPL

REKN

Model parameters

Sp (km2)

0.03

0.04

0.06

Sb/w (km2)

0.01

0.02

0.04

dl (km)

8

12

20

Data

hr (km2)

0.016

2.20

10÷15

hrd (km)

0.12

1.48÷2.40

3.10

〈nd〉 (km)

0.72

0.96

1.44

m (g)

38÷40

50÷60

180÷200

Spand Sb/ware the minimum population and breeding/wintering area, respectively. dlis the estimated maximum dispersal length. 〈nd〉is the neighborhood distance in the breeding season for Snowy Plover (SNPL), and in the wintering season for Piping Plover (PIPL) and Red Knot (REKN). hr and hrd are the home range and the home range distance estimated considering the breeding regions for SNPL, PIPL, and REKN. m is the average body mass. Data are taken from Burney (2009), Himes et al. (2006), Lamonte and Douglass (2002), Pruner (2010), Aiello-Lammens et al. (2011). Model parameters are assumed considering data and calibrating the patch-delineation model on the observed patch size in 2006 derived from the box-counting.

Land-cover model

The land cover is predicted year by year by using the Sea Level Affecting Marshes Model (SLAMM) (Clough 2006; Chu-Agor et al. 2011) starting from the year 2006 to 2100. These simulations are performed in Aiello-Lammens et al. (2011) and Convertino et al. (2010) to which we refer the reader for more details. The domain of the model is extended inland for about 10 km from the coastline (Convertino et al. 2011a;2011b) (black region along the coast in Figure 1, box B1). We consider the predicted inundation distance in 2100 (∼ 9 km) for a range of [1, 2] m sea level rise (SLR) adding 1 km to consider the uncertainty in the estimation of the flooding distance. The initial condition is the 2006 land cover from NOAA (Klemas et al. 1993). The NOAA land cover classes are changed into SLAMM land cover classes for modeling purposes. SLAMM requires us to group the classes of land cover into model classes. The conversion is reported in Convertino et al. (2011b). The SLAMM model also requires the elevation and slope as input variables. The modeled domain is divided into seven regions (Additional file 1: Figure S1) with distinct historical tidal and SLR trends. Each region is characterized by a unique set of values for the 26 input parameters (Additional file 1: Table S1) related to tide, accretion, sedimentation, and erosion processes. The value of the parameters is derived from the available literature and previous efforts of this research (Chu-Agor et al. 2011). In this effort of modeling the land cover, we do not consider any geomorphological feedback between landforms and climate change that is expected to occur with global warming. All our assumptions are the same as those in Chu-Agor et al. (2011) and Convertino et al. (2010). Also we do not consider any possible barrier island shifting because that is reported to occur over a time period much longer than our predictions (Masetti et al. 2008).

Habitat suitability model

The employed habitat suitability model is MaxEnt (Phillips et al. 2006; Phillips and Miroslav 2008), which is one of the most diffused models in species distribution modeling. MaxEnt is a model based on the principle of maximum entropy that predicts continuous habitat suitability maps of potential species occurrence under a set of selected environmental variables. The environmental variables that are necessary and sufficient for calculating the habitat suitability are the land cover translated into SLAMM classes (Chu-Agor et al. 2011; Convertino et al. 2011a;2011b) and the USGS geology layer (Convertino et al. 2011a;2011b) at a resolution of 120 m. The resolution 120 m is the home-range distance of the SNPL (Table 1). Such distance is sufficient to capture not only the spatial variability of habitat preferences of SNPL, but also that of PIPL and REKN, whose home-range distance are much larger than that of SNPL. The habitat suitability at-a-point (i.e., for each pixel of the modeled domain) can be considered as a proxy to find SNPL, PIPL, and REKN in the breeding and wintering season. The prior probabilities of occurrence are calculated in MaxEnt using the recorded shorebird occurrences constrained to the environmental variables. The occurrences are nest and breeding occurrences for SNPL, and the adult occurrences for PIPL and REKN. Thus, for PIPL and REKN the habitat suitability refers to the suitability for wintering as in Convertino et al. (2011a). No absences are required in MaxEnt. Then the posterior probabilities of occurrence are based on the prior probabilities given the change in the land cover modeled year by year by the land cover model. A regularization parameter that controls the fit of the predicted suitability to the real occurrence data is assumed to be equal to one. Non-randomly placed pseudoabsences are used to improve the predictions, and 25% of the occurrences are taken as a training sample (Convertino et al. 2011a;2011b). The predicted habitat suitability maps represent the average of over 30 replicates for each year to reduce the uncertainty of the predictions. The habitat suitability is calculated with 10,000 random background points. Background points are a subset of points of the domain over which the Bayesian inference between the recorded species occurrences, pseudoabsences, and environmental layers is determined.

We assign a biological interpretation to the predicted habitat suitability score, P(hs), which is the probability at-a-point of finding a breeding and/or a wintering ground. Breeding and wintering grounds are suitable sites for the SNPL as a function of the season considered, and wintering grounds are suitable sites for the PIPL and REKN. We define the suitability index (SI) as a metric from 0 to 100 that captures the quality of the breeding and/or wintering habitat for the species. The higher the SI the larger the biological spectrum of functions performed by the species in that habitat. Hence, P(hs) is also a surrogate of habitat use during the breeding and wintering seasons of the species considered. In fact, it is legitimate to assume that habitat use increases with habitat quality. Every pixel of the HS maps is classified into five SI categories: SI=100 [for 0.8≤P(hs)≤1] is considered the best habitat with the highest survival and/or reproductive success; SI=80 [for 0.6≤P(hs)<0.8] is typically associated with successful breeding and/or wintering; SI=60 [for 0.2≤P(hs)<0.6] is associated with consistent use for breeding and wintering; SI=30 [for 0.2<P(hs)] is associated with occasional use for non-breeding, feeding activities, and wintering; all values less than SI=30 indicate habitat avoided both for breeding and wintering; and SI=0 for completely unsuitable habitat. We refer the reader to Convertino et al. (2011a) for additional details about MaxEnt runs for the SNPL, PIPL and REKN.

Patch-delineation model

Below we define a criterion for delineating breeding and wintering patches for SNPL, PIPL and REKN, respectively. A patch is defined when the following criteria simultaneously hold:

The species-dependent values for the three parameters required for the patch identification are reported in Table 1. The values of biological data in Table 1 are used only to support the choice of model parameters. The model parameters are calibrated to reproduce a patch-size distribution as close as possible to the box-counting distribution of occurrences in 2006. The model with this set of parameters was validated against the patch-size distributions from 2002 to 2010 estimated by the box-counting. We define breeding patch as an area large enough to at least occasionally support a single breeding pair through courtship and rearing of young to dispersal age (Majka et al. 2007). A population patch is defined as an area large enough to support breeding for 10 years or more, even if the patch is isolated from interaction with other populations of the species (Majka et al. 2007). Since population-wide data are lacking for these breeding and population area requirements, we assumed that a population patch is at least two times larger than a breeding patch. For the SNPL these patches contain certain nesting patches. The minimum population and breeding/wintering patch areas are estimated from the literature available and by expert knowledge of the field biologists involved in the sampling campaigns performed for this study (see Burney 2009; Himes 2006; Lamonte and Douglass 2002; Pruner 2010). Sp and Sb/w are the minimum population and breeding/wintering area, respectively, and are proportional to the estimated home range. The minimum breeding/wintering area is the minimum area that will support breeding and wintering activity of the shorebirds. The home range hr and the home-range distance hrd (the square root of the hr) are values estimated considering the breeding regions for SNPL, PIPL, and REKN. We assume that Sp and Sb/w for PIPL, and REKN are much smaller than hr because they refer to the wintering period of these shorebirds in Florida. For REKN, Spis also reduced due to the habitat limitation and the close coexistence with SNPL in the same habitat. Patches are considered connected if their neighboring distance is equal to or smaller than dl, which is the maximum dispersal length. Figure 1c,d shows an example of patches that are connected because their reciprocal distance is lower than dl. These plots also represent our assumption that coastline complexity affects patch distribution. The average neighborhood distance 〈nd〉 is the average dispersal of the species. 〈nd〉 is higher than hrd for the SNPL due to the higher local dispersal ability estimated from recent surveys (Himes et al. 2006; Pruner 2010). For PIPL and REKN, 〈nd〉is smaller than hrd because the reported hrd refers to their breeding range in northern states in the USA and Canada. In the winter season PIPL and REKN migrate to Florida, and their dispersal distance is observed to be smaller. Within the neighborhood distance a subpopulation can be assumed to be panmictic. A panmictic population is one in which all individuals are potential partners. It is usually estimated from the foraging distance of an animal species. In a more abstract way the neighborhood distance is the glue of all the suitable patches. In a particle physics analogy, it describes the Brownian motion of individuals within a larger species group. Thus, by using dl, which is the maximum dispersal, as a criterion in the model, foraging is certain to be considered within patches. Our model considers an upper estimate of the patch size for all the shorebirds considered. m is the average body mass, which is used to discuss some results. We assume the same biological parameters for the SNPL Panhandle and Peninsula as in Aiello-Lammens et al. (2011).

Probability distribution of the patch size

The probability of exceedance of the patch size is known in literature as Korčak’s law (Korcak 1940; Nikora et al. 1999), which is expressed by:

P(S≥s)=cS−εFssc,

(3)

where c is a constant, F is a homogeneity function that depends on a characteristic size sc, and ε=DK/2 is the scaling exponent (Korcak 1940; Mandelbrot 1982). DKis the fractal dimension of the patches. The probability of exceedance exhibits a power-law behavior. The probability distribution of the patch size for the predicted patches was used to validate the patch-delineation model against the box-counting estimates on the real occurrences from 2002 to 2010. The fit of the predicted distribution of patches is performed using a Maximum Likelihood Estimation technique (MLE), which is described in the Additional file 1.

Perimeter-area relationship

The scaling relationship between the perimeter p and the size S of the patches:

p=kSDc/2,

(4)

determines the fractal dimension of the mosaic of patches, which considers the fractality of the patch edge. Here we indicate the fractal dimension Dc, which is derived from the same predicted patches of the introduced patch model (see the “Patch-delineation model” section) but also considers their perimeters. Because Korčak’s law (Korcak 1940) considers only the size of the patches, the perimeter-area scaling law has been considered as a more precise tool for measuring the fractal dimension. In literature the ratio p/S is adopted to measure the quality of the patches for population survivability, that is, the likelihood of surviving in a suitable patch (Helzer and Jelinski 1999; Airoldi 2003; Imre and Bogaert 2004). In general the higher the ratio p/S the less suitable the patch area for the species, and the higher the ratio p/S, the higher the fractal dimension Dc.

The relationship between the number of cells occupied by shorebird occurrences, N(l), and the length of the side of the box, l, at each scale of analysis is shown in Figure 2. The relationship is a power-law function, N(l)∼lDb, whose exponent Db is the fractal dimension of the shorebird occurrence pattern. Figure 2a reports the power-law relationship for PIPL and REKN breeding occurrences in 2006, and Figure 2b for the SNPL breeding occurrences from 2002 to 2010. The results confirm the supposed scale-free distribution of the shorebird occurrences. The fractal distribution of the predicted patches is captured by Korčak’s law (Figure 3). The box-counting overestimates the fractal dimension with respect to the fractal dimension of Korčak’s law as shown in Convertino et al. (2012). The fractal dimension of the box-counting (Db) is 1.63, 1.85, and 1.53, and the fractal dimension of Korčak’s law (DK) is 1.47, 1.70, and 1.42 for SNPL, PIPL, and REKN in 2006, respectively (Additional file 1: Tables S2 and S3). The box-counting envisions a more pessimistic scenario for the patch size of shorebirds. However, as in Convertino et al. (2012) we believe that in the absence of any modeling effort box-counting constitutes a valid technique to calculate the fractal dimension of the mosaic of patches. For the SNPL occurrences, box-counting allows us to detect the fluctuation over time of the fractal dimension of the recorded nest occurrences and of the coastline. The insets in Figure 2b show the empirical evidence of the correlation between Df and Db, and Additional file 1: Table S3 reports the values of the fractal dimensions. The analysis raises the question of whether the variation in Dbis caused by natural fluctuations of the species range or by changes in external forcing such as natural or anthropogenic stressors. We observe that in 2004 and 2005 the fractal dimension showed a jump possibly due to the exceptional hurricane season in those years, which altered the positive feedback between tropical cyclones and SNPL nest abundance (Convertino 2011c). This supposition is confirmed by the results of Convertino et al. (2011c).

Figure 2

Box-counting scaling-law in time. (a) Power-law N(l)=N0l−Dbderived from the box-counting algorithm applied to the occurrences of PIPL (black dots) and REKN (green) in 2006 and to the whole Florida Gulf coastline. In the inset the schematized Florida coastline is evaluated at different box sizes. (b) Box-counting algorithm applied to the 2002-2010 occurrence of the SNPL. The fractal dimension derived from the analysis of the breeding and nesting occurrences is Db=1.63, 1.62, 1.75, 1.74, 1.63, 1.64, 1.66, 1.68, 1.70 for the years from 2002 to 2010, respectively. In the inset Dband Dfare reported for each year. Values of Dfare reported (Additional file 1: Table S3).

Figure 3

Korčak’s law of the predicted suitable patches. The fractal dimension of the patches is derived from the scaling exponent, ε=DK/2, of the probability of exceedance of the patch size (Equation 3) for SNPL, PIPL, and REKN. The probability of exceedence of the patch size is represented for the years 2006, 2020, 2040, 2060, 2080, and 2100. The probability of exceedance is compared against the box-counting scaling laws for the years 2002-2010. Additional file 1: Table S2 reports the values of DK. The insets represent the probability density functions (pdfs) of the patch size that show a heavy-tailed behavior.

The potential effect of sea level rise, one of the main controlling factors of land cover of coastal habitats, is studied here. The simulated variation in land cover classes over time is performed in SLAMM (Clough 2010) for the Gulf Coast of Florida (Additional file 1: Figures S1 and S2). We predict by 2100 a decrease in the salt-marsh and estuarine beach classes, which are crucial habitats for PIPL, SNPL, and REKN. We also predict a net decrease in swamp and inland fresh marsh habitats. Following flooding predicted to occur after 2060, undeveloped drylands will change mostly into tidal flats, which may shift into estuarine open water (Additional file 1: Figure S2). We estimate a 6% increase in estuarine open water and a 10% increase in ocean open water from 2006 to 2100. We expect global land-loss independent of the land cover class of about 16% with respect to the 2 m sea level rise. A video in Additional file 2 and Figure S2 in Additional file 1 show the evolution of land cover and of the coastline geomorphology over time. Additional file 1: Figures S3, S4, and S5 report the suitability index derived from the predicted habitat suitability maps using MaxEnt corresponding to the yearly land cover maps. The patches are then calculated using the patch-delineation model introduced in the “Patch-delineation model” section and the habitat suitability maps.

The power-law structure of the patch size holds for every year simulated (Figure 3), which proves the scale-invariance of the suitable habitat over time. By using the maximum likelihood estimation (MLE) criteria, we found that the Pareto-Lévy probability function has the best fit for the predicted distribution of the patch size (Additional file 1). Korčak’s law exhibits some finite-size effects before the upper truncation and a potential lower-cutoff in the power-law behavior. However, these variations from the power law are quite common in natural systems due to the finiteness of the variable sampled. Thus, we can claim an overall scale-invariance of the patch size. Additional file 1: Table S2 reports the fractal dimension derived from Korčak’s law for 2006, 2020, 2040, 2060, 2080, and 2100. The scale-invariance of the habitat patterns of the SNPL was shown in Convertino et al. (2011b) for the prediction of the habitat suitability in 2006. Here we show that, given the scale-invariance of the patch size, fluctuations in the scaling exponent ε=DK/2 of Korčak’s law occur. We believe that these fluctuations are related to variations in the land cover, which changes the coastline fractality. The higher the fractal dimension, the higher the fragmentation of the shorebird habitat. The fragmentation of the habitat creates smaller patches for wintering and breeding for PIPL and REKN, and for SNPL, respectively. Brownian-Lévy movements of shorebirds might be the cause for the scale-invariance of the occurrence patterns that can be detected by the box-counting. This has been proven for other marine animals (Humphries et al. 2010) and colonial birds (Jovani et al. 2008). However, in this study we do not reproduce any movement of species as we believe that the size and number of patches is affected by the geomorphological evolution of the coastline, which in turn affects the movement of shorebirds.

The worst scenario for the vulnerability of shorebirds is predicted considering the fractal dimension Db from the box-counting. Moreover, the box-counting suffers from the risk of potentially unsampled occurrences. The Korčak’s law fractal dimension (Figure 3) is based on the size of the predicted suitable patches (potential habitat range), while the box-counting (Figure 2) is an approximation that only captures the recorded occurrences (realized range). The fact that DK≅Db for the 2006-2010 period in which SNPL nest occurrences are available confirms the good estimation of MaxEnt of the realized range as previously found in Convertino et al. (2011b). A more accurate estimation of the fractal dimension that is intermediate between DKand Db is given by the patch perimeter-size scaling relationship (Figure 4). The perimeter-size relationship captures the edge effects of patches on species. In general shorebird species prefer to live in patches whose shapes are as regular as possible versus highly irregularly shaped patches such as the patches determined by a very complex coastline. The survivability of the species is higher for those inhabiting patches with large perimeters and simple shapes than for those inhabiting patches of equivalent area but complex shape. The larger the edge effect determined by the complexity of the patch parameter, the lower the probability of survival for the individuals within the species within the patch. However, there are some cases of “edge species” for whom irregular shapes are preferred. In our case it was observed that DK≤Dc≤Db. Hence, the estimation of the fractal dimension by using Korčak’s law forecasts the best scenario predicting the least amount of fragmentation due to sea level rise. Dc predicts greater fragmentation than DKbecause the fractality of the patch’s perimeter is considered, but overall Dc seems the best estimate of the fractal dimension between Korčak’s law and the box-counting estimates.

Figure 4

Perimeter-size relationship for SNPL, PIPL, and REKN. Perimeter-size relationship (p=kSDc∗/2) for the predicted suitable patches of the SNPL (a), PIPL (b), and REKN (c), in 2006 and 2100. The exponent Dcfor the SNPL is listed in Additional file 1: Table S3.

Figure 5a,b respectively, show the time series of the fractal dimension of the species-dependent habitat coastline Df (mostly beach for SNPL, PIPL, and REKN, but also salt marsh for the PIPL), and of the fractal dimension of the patches DK (from Equation 3) computed with the patch-delineation model. Additional file 1: Figures S3, S4, and S5 show the patches for SNPL, PIPL, and REKN in the years 2006, 2020, 2040, 2060, 2080, and 2100. The majority of patches are along the barrier islands and particularly in the Panhandle region. After 2060, when sea levels start to rapidly rise, a consistent portion of the patches will be found along the shore as barrier islands gradually disappear. Figure 5a also shows the variation in the fractal dimension of the whole coastline independently of the land cover class. The probability of finding a patch of size S is lower in 2100 than in 2006. DK values are similar for SNPL and REKN and are higher for PIPL (Figure 5b). Thus, on average the relationship DcREKN≤DKSNPL≤DKPIPL holds for the modeled period. Big variations in DKPIPL are observed particularly in correspondence with big variations in the salt-marsh habitat (Figure 5a,b), which confirms the likelihood of finding a breeding ground of PIPL in the salt-marsh habitat (class 8 contained in Additional file 1: Figure S6) as reported in literature (Convertino et al. 2011a) and as found by our results (Additional file 1: Figure S6). In 2100 the fractal dimension of REKN is very similar to the fractal dimension of SNPL, while the fractal dimension of PIPL is the highest. The PIPL shows the lowest probability of large patches with respect to the other shorebirds considered because DK is the highest. The area under the power-law distribution of patches for 2100 in Figure 3 has a 5, 3 and 8% negative variation with respect to the area for 2006 for SNPL, PIPL and REKN. The area under the curve is the overall probability of finding patches of any size in a given year. Just comparing 2006 with 2100 is not enough to derive any conclusion about the species with the highest potential risk of decline. The location and size of the patches are determined by the habitat suitability at-a-point and by a combination of dispersal and area criteria (see the “Patch-delineation model” section). The Piping Plover, despite having a larger spectrum of habitat preferences than SNPL and REKN [transitional marsh and salt marsh areas are favorable classes as shown in Additional file 1: Figure S6 and as reported in Convertino et al. (2011a)] seems to be at risk due to the high fragmentation of its habitat. This is evidenced by the larger fluctuation of DKPIPL than DKSNPL and DKREKN.

Figure 5

Fractal dimension time series of the shorebirds patches and of the coastline. (a) Time series of the fractal dimension Dfof the entire coastline (blue line), of the salt-marsh (red), and of the beach (green) habitat coastlines, determined by the box-counting algorithm. (b) Fractal dimension DKover time for the patches for SNPL (blue dots), PIPL (red), and REKN (green) derived from Korčak’s law. The dashed gray lines (a, b) represent the 95% confidence interval of the estimated Dfand DK. (c) Scaling relationship among the fractal dimension of the patches for the threatened, endangered, and potentially at-risk shorebird species (TER-s) and the fractal dimension of the favorable habitat coastlines (salt marsh for PIPL, and beach for SNPL and REKN). The average species-independent scaling exponent is γ= 1.67. The gray cloud (c) represents the 95% confidence interval for the linear regression between DKand Df.

We believe that it is important to observe the fluctuations of DK over time for each species. DKvalues of SNPL and REKN are on average steady and increasing over time, respectively; thus, the probability of finding large patches for these shorebirds decreases over time with respect to 2006. DK of PIPL has the largest fluctuations, but most of these fluctuations imply an increase in the probability of finding large patches with respect to 2006. Nonetheless we believe the frequent and large variation in patches is not a good scenario for species.

In Figure 5c we propose a scaling relationship between the fractal dimension of patches and the fractal dimension of the habitat-specific coastline, DK∼Dfγ. The relationship holds over at least two orders of magnitude, from the smallest patches (∼ 0.01km2) and short coastline segments to the largest patches and the whole Florida Gulf coastline. The same scaling exponent is observed for SNPL, PIPL, and REKN, underlining a possible common ecogeomorphological organization of the landscape under sea level rise pressure. In Figure 5cDf is characteristic of the portion of coastline in which there is a suitable habitat for SNPL, PIPL, and REKN, which is evidenced in Additional file 1: Figure S6. The coupled evolution of the land cover and habitat patterns may hold clues about the linkage of geomorphological and ecological processes. The scaling relationship between the fractal dimensions of patches and coastline can be a potential tool to measure the vulnerability of the species in the future. The higher the exponent γ, the higher the potential risk of decline of the species. For small changes in the configuration of the coastline, a large fragmentation of the suitable habitat would potentially be observed. For species with comparable values of γ, which is the case for SNPL, PIPL, and REKN, the range of values of DK and Df is important for detecting which species may be subjected to the most significant change in the suitable habitat patches. The lower DK, the higher the likelihood of having large patches. To the best of our knowledge this is the first scaling relationship to be identified between fractal dimensions of landscape and ecological patterns. In this respect this relationship brings insights into the field of “landscape allometry,” which is the study of the possible scaling of landscape and ecological patterns and processes. The relationship is between fractal dimensions, which are indicators that focus on how measured quantities vary as a power of measurement scale, but at the same time the relationship has an allometric focus, between the coastline complexity and the magnitude of habitat fragmentation.

However, fragmentation per se does not directly imply loss of connectivity among patches. Figure 6 shows how the average size of the patches 〈s〉for SNPL, PIPL, and REKN decreases with the increase in the fractal dimension of the patches. Here we consider DK of Korčak’s law for the fractal dimension. At the same time we observe an increase in the number of patches Np. Thus, the variation in the coastline produces fragmentation, rather than shrinking, of the suitable habitat. The former does not imply the latter as erroneously assumed by many theoretical models in the ecological literature. The average size of the PIPL patches is lower than that for SNPL and REKN, and the habitat for the PIPL is the most fragmented (Np is the highest on average). This is related to the high value of DKfor the PIPL with respect to SNPL and REKN. Thus, although the variations in DKPIPL would predict bigger patches, the fragmentation of the PIPL habitat is the greatest. In 2100 the number of suitable patches for SNPL, PIPL, and REKN is predicted to be higher than in 2006, but the average size of the patches is predicted to be smaller (Additional file 1: Table S2). As sea level rise (SLR) increases the complexity of the coastline, habitat patches moderately shrink and split. On the contrary when the coastline complexity decreases, habitat patches enlarge and coalesce (Figure 1c) as in our assumption depicted in Figure 1b. The PIPL seems to be the shorebird most affected by the changes in its breeding habitat due to sea level rise.

Figure 6

Relationships among patch number, size, and connectivity, and fractal dimension of the habitat-specific coastline. 〈s〉vs DK(a), Npvs DK(b), Npvs 〈s〉(c), and 〈c〉vs DK(d) for the threatened, endangered, and at-risk shorebirds (TERs) considered. The dots are the bin averages over 30 simulations for each year for the period 2006-2100. The dashed lines represent the 95% confidence intervals for the dependent variables considered.

The average size and the number of the patches are inversely proportional given the relationship in Figure 6a,b and as shown in Additional file 1: Figure S7. The average patch size 〈s〉 for the shorebirds is not proportional to the average body mass m as possibly expected (Table 1), although the latter scales with the average dispersal length. The 〈s〉is for the PIPL, while it is larger for SNPL and REKN. This emphasizes the controlling role of habitat geomorphology in shaping the patch distribution. The PIPL also depends on the salt-marsh habitat, which t is one of the classes more seriously compromised by SLR. We consider dl, the estimated maximum dispersal length, in order to determine the average number of connected patches 〈c〉. dlconsiders rare “Lévy flights” of individuals of the species in the ecosystem. Lévy flights are a special class of random walk with movement displacements drawn from a probability distribution with a power-law tail (the so-called Pareto-Lévy distribution), and they give rise to stochastic processes closely linked to fractal geometry and anomalous diffusion phenomena. Because it has the largest maximum dispersal distance, the REKN has the highest number of connected patches. However, for the three shorebird species 〈c〉increases with the fractal dimension of the patches, indicating a measure of the habitat fragmentation. Because we find that climate change is responsible for the splitting of the patches, rather than their shrinking, and because the dispersal capability of species is not expected to change consistently in the modeled period, the result seems justifiable. The increase in the number of connected patches is explainable because Npincreases without a drastic reduction in the habitat. The average connectivity of the predicted breeding and wintering patches is an increasing function of the fractal dimension of the patches. The increasing roughness of the Florida coastline due to climate change produces a larger number of patches with smaller dimensions. The increased connectivity would potentially enhance the survivability of the shorebirds despite the decrease in the average size of suitable patches. Thus, the predicted patch patterns for the Florida shorebirds are not the worst case scenario in which both the connectivity and the dimension of the patches are reduced. Further explanation of the land cover, habitat, and patch dynamics is provided in Additional file 1.

Sea level rise due to climate change, beyond being a human-population threat, is shown to strongly affect biodiversity such as residential and migrant shorebird populations in Florida. The integrated patch-prediction modeling framework proposed in this paper constitutes a parsimonious but useful risk assessment tool for species decline with respect to more accurate metapopulation models. In our opinion, the understanding of ecogeomorphological processes at any scale of analysis together with the detection of useful indicators of such dynamics is one of the primary goals to protect biodiversity against the anticipated changes in the landscape due to climate change. On the one hand, it is impossible to consider, or to estimate with low uncertainty, all the factors affecting the processes that govern the distribution of species (e.g., conspecific attractions, interspecific competition, density dependence, sex structure, life history, phenotypic plasticity, and phenological changes in dispersal ability and in breeding/wintering area requirements), the geomorphological processes, and the links and feedbacks among these processes. On the other hand, we believe that a top-down approach of biocomplexity is useful to detect the fundamental drivers of the observed patterns of interest (Schwimmer 2008; National Research Council 2009; Reinhardt et al. 2010). We are aware that many geomorphological and biological processes are not incorporated in the presented model; however, the uncertainty in the quantification of these processes and the interaction of these uncertainties may produce erroneous results in the predictions. The integrated model is capable of providing valuable macroscale predictions with relatively few data and variables. Thus, the model is useful for evaluating conservation actions for increasing the survivability of shorebirds in Florida. We are also confident that the proposed model, properly tuned, can be applied to many different species in coastal ecosystems worldwide that are threatened by sea level rise. We anticipate further development of this model at higher levels of complexity and also for inland sites. The following conclusions are worth mentioning.

A scale-free distribution of nesting, breeding, and wintering occurrences was detected for the Snowy Plover in Florida. The scale-free distribution was also found for the wintering occurrences of Piping Plover and Red Knot. The distribution was derived through the box-counting technique applied to the breeding and wintering occurrences, which gives a proxy of the fractal dimension of shorebird patches. Empirical evidence shows that the fractal dimension of the occurrences is strongly positively correlated with the coastline fractal dimension, which underlines an ecogeomorphological organization, i.e., a coupling of ecological and geomorphological patterns. The power law held for any season of the shorebird annual cycle, demonstrating the high importance of the physical habitat on species processes.

We predicted breeding and wintering patches of shorebirds, simulating land cover (which comprises many coastal wetland types) and habitat suitability at the year scale from 2006 to 2100 as a function of sea level rise. Patches were identified by a set of macroecological criteria, such as area, habitat suitability, and neighboring distance, as a function of the maximum dispersal. The distribution of the predicted patch size was Korčak’s law, whose exponent is half of the fractal dimension of the patches. We validated the model by predicting the observed patch-size distribution and patch patterns from 2002 to 2010 where data were available. We also investigated the perimeter-size relationship for estimating the fractal dimension of the patches at a higher level of complexity because of the calculation of the perimeter. The fractal dimension provided by the perimeter-size relationship provided a median estimate between the values derived from Korčak’s law and the box-counting distribution. Korčak’s law provided the most optimistic scenario of fragmentation in which the probability of finding large patches was the highest, while the box-counting provided the most pessimistic scenario. Hence, the perimeter-area relationship is suggested as the best method to calculate the fractal dimension of the mosaic of habitat patches.

The robustness of the Pareto-Lévy distribution of the patch size was verified for predictions of patches from 2006 to 2100. Thus, the scale-invariance of the patch patterns holds in time despite the strong influence of sea level rise. This may be related to a sort of simulated “biological resilience” of species to the external changes (Folke et al. 2004) by assuming invariant habitat area and dispersal requirement. Scale-free habitat patterns have proven to be the most resilient to external stressors in previous studies (Kefi et al. 2011). Thus, the shape of the patch-size probability and the fractal dimension when this probability is a power law can be useful indicators to estimate the “degree of stress” of coastal ecosystems. Further research is anticipated to understand when and how the patch-size probability deviates from a Pareto-Lévy behavior. The fragmentation, which is proportional to the fractal dimension of the habitat-specific coastline, varied considerably over time and in particular for the Piping Plover. However, the risk of extirpation in 2100 for SNPL, PIPL, and REKN was not high with respect to 2006. We note that the comparison between final and initial years’ risk should not be the only comparison in evaluating the risk of decline of a species. The overall trend of the fractal dimension in the modeled period has to be evaluated as well.

A scaling relationship was found between the fractal dimensions of the patches and of the habitat-specific coastline. The scaling exponent of this relationship appears to be species-independent for the shorebirds considered. Further research is needed to explore the conditions of universality (species- and ecosystem-wise) of this relationship, which may be related to the species considered. The fluctuation in the fractal dimension of the coastline can be assumed to be a valuable ecological indicator for assessing variation in patch patterns of breeding and wintering shorebirds.

We demonstrated that habitat loss, fragmentation, and connectivity are three separate concepts. Although these variables are closely linked to each other, their causality is not trivial. For the shorebirds studied, the predicted fragmentation was coupled with habitat loss while the connectivity increased. The fact that the patches, even if smaller, were connected is an extremely positive factor that ensures dispersal and gene flow; thus, the connectivity of patches enhances the survivability of shorebirds. Birth, death, and dispersal processes of a species can overcome the habitat-loss effect and a decrease in the average size of patches. Yet, a lower metapopulation risk of extirpation exists if interpatch migration is allowed (Kindvall and Petersson 2000). However, a decrease in the average patch size can potentially increase intra-species competition for foraging (Ritchie 1998) and decrease carrying capacity. A possible optimal ecogeomorphological state of the coastal ecosystem may be characterized by the smallest fractal dimension of the coastline that maximizes the compactness of the suitable patches. This configuration also minimizes the fractal dimension of the patches. The highest entropy of this configuration may translate into the smallest energy expenditure of the species that inhabit the habitat, for example, for foraging and breeding activities. The entropy of geomorphological landforms (Nieves et al. 2010) may, in fact, be highly correlated with the scale-invariance of ecological patterns such as species-patch patterns.

MC is Research Scientist at the University of Florida, Gainesville, and a Contractor of the Engineering Research and Development Center of the US Army Corps of Engineers at the Risk and Decision Science Team. AB is currently a financial analyst at Frontier Airlines. AB got his B.Sc and M.Sc. from MIT, Civil and Environmental Engineering program. AB performed his research internship at the Risk and Decision Science Team in the summer of 2011. GAK and RMC are Associate and Professor at the University of Florida, Gainesville, respectively. IL is team leader of the Risk and Decision Science Team of the Engineering Research and Development Center of the US Army Corps of Engineers.

Acknowledgements

This research was supported by the US Department of Defense, through the Strategic Environmental Research and Development Program (SERDP), Project SI-1699. M.C. acknowledges the funding of project “Decision and Risk Analysis Applications Environmental Assessment and Supply Chain Risks” for his research at the Risk and Decision Science Team. The computational resources of the University of Florida High-Performance Computing Center (http://hpc.ufl.edu) are kindly acknowledged. The authors cordially thank Dr. RA Fisher (Engineering Research and Development Center of the US Army Corps of Engineers) and the Eglin Air Force Base personnel for their help in obtaining the data and for the useful information about the breeding information of SNPL. Tyndall Air Force Base and Florida Wildlife Commission are also gratefully acknowledged for the assistance with the data. We thank M.L. Chu-Agor (currently at the Center of Environmental Sciences, Department of Biology and Earth and Atmospheric Sciences, Saint Louis University, St. Louis, MO) for her computational effort with SLAMM at the University of Florida. Permission was granted by the USACE Chief of Engineers to publish this material. The views and opinions expressed in this paper are those of the individual authors and not those of the US Army or other sponsor organizations.

Competing interests

The authors declare that they have no competing interests.

Author’s contributions

MC designed the study, managed and analyzed the data, wrote the model (box-counting and patch delineation model), developed the theory, and wrote the manuscript. AB assisted in making the calculations and analysis, and helped in writing the manuscript. GAK and RMC participated in the habitat suitability modeling framework and reviewed the manuscript. IL supervised the whole work, and reviewed the manuscript by providing a practical angle to this research for effective environmental management. All authors read and approved the final manuscript.

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License(http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.