J. theor. Biol. (2000) 202, 25-32 Article No. jtbi.1999.1030, available online at http://www.idealibrary.com on IDEM Model Assessments of the Optimal Design of Nature Reserves for Maximizing Species Longevity Jon D. Pelletier* Department of Geosciences, University of Arizona, 1040 E. Fourth St., Tucson, AZ 85710, U.S.A. (Received 2 November 1998, Accepted in revised form on 28 September 1999) Using a computational model for the population growth and dispersal of a model species in a fluctuating environment, we test three nature reserve geometries (one large, many small, and a self-similar distribution of reserve sizes) to determine which geometry maximizes species longevity. The self-similar distribution is a close approximation to the distribution of managed areas in the conterminous United States. We consider models with and without migration from or between reserve fragments and both short- and long-range dispersal mechanisms. The optimal geometry depends on the type of dispersal and on the relative probability of survival in protected and non-protected areas. When no migration is allowed from or between reserve fragments of the three geometries, many small equally sized reserves are the optimal geometry. When migration is allowed, the optimal geometry is a single large reserve when the survivability in non-protected areas is low and a self-similar distribution when the survivability is high. © 2000 Academic Press 1. Introduction Ecologists, biogeographers and conservationists have long debated about the optimal design of nature reserves with respect to species longevity. The principle questions are "How big? How many? How arranged?" (Soule & Simberloff, 1986). Although there is no single answer to these questions since the optimal geometry of a nature reserve always depends on variables unique to the set of species being protected, general patterns that result from simple optimization criteria for reserve networks of a fixed total area and uniform habitat may provide useful guidelines in many cases. The debate on the optimal design of nature reserves was begun by Wilson & Willis (1975) and Diamond & May (1976) who argued that the *E-mail: jon@geo.arizona.edu 0022-5193/00/010025 + 08 $35.00/0 equilibrium theory of biogeography implied that a single compact reserve was optimal. More recently, it has been recognized (Abele & Connor, 1979; Gilpin & Diamond, 1980) that the equilibrium theory is neutral on the questions of how big?, how many? and how arranged? and the argument has shifted to determining minimum viable population sizes within nature reserves (Soule & Simberloff, 1986). Archipelagos have, in most cases, a greater number of species than a single island of the same total area (Soule & Simberloff, 1986; Simberloff & Abele, 1982). This does not necessarily suggest, however, that archipelagos of nature reserves are an optimal design since the total number of species may be an inadequate optimality criterion and because we must consider the dynamic evolution of species abundances in reserves to determine the likelihood of survival (Soule & Simberloff, 1986). There is also no simple answer to the question of © 2000 Academic Press 26 J. d. pelletier which geometry maximizes the longevity of a species or a group of species. Widespread species in large reserves are likely to enjoy a smaller probability of extinction (McKinney, 1997) and lower abundance volatility (Boulinier et al, 1998; Glazier, 1986). However, there are potentially overriding advantages to having several small reserves instead. Although the probability of local extinction will be higher for species in small reserves compared to large ones, it is necessary for the species to become locally extinct in each and every small reserve rather than only the one large reserve for total extinction to occur in the case of many small reserve fragments. In addition, if migration is allowed between reserve fragments, the "rescue effect" may significantly increase the persistence of species (Schoener & Spiller, 1987; Gonzalez et al, 1998; Burkey, 1989; Hastings, 1982). The efficacy of the "rescue effect" will necessarily depend on whether the dispersal is long or short ranged and on the probability of survival in non-protected areas if the reserve fragments are not connected. Thus, to answer the question of which geometry maximizes the time to extinction of a species we must simulate realistic models of the variation of species abundance in space and time in a variety of reserve geometries. An alternate approach to simulating realistic models of the variation in species abundance is the development of algorithms for optimal reserve selection and design based on theoretical models or empirical data on survivability (Pressey et al, 1999; Clemens et al, 1999). 2. Modeling Population Variability Abundance variations in population dynamic models can be driven by three types of random processes: demographic stochasticity, environmental stochasticity, and catastrophic stochasticity. In the case of demographic stochasticity alone, the probability for an individual in a species to reproduce or die is equal for all individuals and constant in space and time. For a very small number of individuals demographic stochasticity can result in large fluctuations of population density. However, the variance becomes very small for large abundances (variance decreases as 1/JV where JV is the effective population size) and will probably be dominated by environmental stochasticity for all but the smallest populations. With environmental stochasticity, the net population growth rate is a random variable in space and time: the net population growth from one generation to the next is a(t + 1) — a(t) = t](x, y, t)a(t), where a(t) is the local abundance at time t and t](x, y, t) is an uncorrelated Gaussian noise in space and time representing many complex factors of the physical and biotic environment including variations in climate, availability of nutrients or prey, and degree of predation by other species (Sousa, 1984). Catastrophic stochasticity occurs when the net population growth rate is spatially correlated over large distances. As an example, the effects of forest fires can be modeled by wiping out all vegetation at a random location within a random area drawn from the probability distribution corresponding to the frequency-size distribution of forest fires. Wright & Hubbell (1983) considered the problem of species persistence in a model with demographic stochasticity. They found that one large reserve was optimal when migration out of reserve fragments was not allowed and that the average species longevity of one large and two small reserves was nearly the same when migration was allowed. Goodman (1987) argued that environmental stochasticity was likely to be dominant over demographic stochasticity in most natural cases and found that the "rescue effect" significantly increased persistence in fragmented reserves for a model with environmental stochasticity. We will consider all three types of stochasticity in our model. 3. Model Geometries We will consider three reserve geometries representing two extreme cases and one intermediate case: (a) one large square (or nearly square: 15 x 14 + 6 cells), (b) many small squares of equal size, and (c) a self-similar distribution formed by a Sierpinski carpet (Mandelbrot, 1983) construction (Fig. 1). Figure 1(a) and (b) deviate slightly from square geometries in order to make their total areas equal to that of (c). Black regions in Fig. 1 (except for the border) are protected areas while the remaining white space of the 60 x 60 cell grid is the non-protected area. The unprotected areas in our models have a smaller model assessments of nature reserves 27 (a) (b) (c) Fig. 1. Geometries of the reserve networks considered: (a) one large nearly square (15 x 14 + 6 cells), (b) many small squares of equal area, and (c) a self-similar distribution formed by an order 3 (33 x 33 cells) Sierpinski carpet. Black areas (except for the border) are protected areas and all other areas are non-protected. The total protected area and non-protected area is the same in (a), (b), and (c). net population growth rate to model the effects of habitat destruction in those areas. The distribution of reserves in Fig. 1(c) may seem contrived. However, we have chosen to model the Sierpinski carpet because a stochastic version of the Sierpinski carpet corresponds closely to the distribution of managed areas in the conterminous United States. The managed areas principally include wildlife refuges, state parks, national parks, and Native American reserves (Fig. 2). We first consider a more realistic, stochastic version of Fig. 1(c) illustrated in Fig. 3. In the deterministic construction of Fig. 1(c) a square of 27 x 27 cells is divided into nine 9x9 squares. The central square is made a protected area and the remaining eight squares are further subdivided as the original 27 x 27 square was. In the stochastic Sierpinski carpet construction of Fig. 3, which of the nine squares at each iteration is made into a continuous reserve area is determined by a random selection. The cumulative frequency-size distribution of connected areas, the number of areas greater than or equal to A, of the stochastic Sierpinski carpet is a power-law function with an exponent of —0.8 [Fig. 4(a)]. The corresponding distribution for all connected managed areas of the conterminous United States is identical [Fig. 4(b)]. The distribution of managed areas was given by McGhie et al. (1996). In the distribution of Fig. 4(b) the managed areas that share a common border were joined to make a single larger area. One major difference between the model geometry of Fig. 3 and the managed areas of the conterminous United States is that the model areas are assumed to be fragments in a uniform habitat whereas the United States contains many different habitats. Nevertheless, since the distribution of managed areas is scale-invariant, we can expect the distribution of reserve fragments to be similar at smaller spatial scales where habitats are more uniform. As we will show, this distribution is an optimal geometry in the case when successful migration between reserve fragments can occur. The opti-mality of self-similar nature reserves in this case is consistent with the self-organization of species occurrence in nature into a power-law frequency-size distribution of areas (Hastings et al, 1982; Kunin, 1998). Cumulative frequency-size distributions of areas with an exponent between —1 and —\ arise in a number of geographical contexts such as the size distribution of lakes and islands (Mandelbrot, 1983). One possible reason for the 28 J. d. pelletier Fig. 3. Geometry of a one realization of a stochastic order 5 (35 x 35 cells) Sierpinski carpet. occurrence of this distribution is that regions enclosed within contour lines of a self-affine fractal surface have a power-law distribution of sizes. In the case of islands and lakes this surface is the self-affine fractal topographic surface of the Earth (Turcotte, 1997). Self-affine fractal surfaces can be quantified with the power spectrum. The power spectrum S(k) of the height of a self-affine fractal surface has a power-law dependence on wave number k with exponent —ß: S(k) cck~ß. For the Earth's topography ß x 2 is observed (Sayles & Thomas, 1978). Other ecologically relevant variables within landscapes besides elevation have also been modeled as self-affine fractals (Gardner et ah, 1987; Krümmel et ah, 1987; Milne, 1991). Nature reserves are likely to be chosen, in part, because within them certain environmental factors are above or below a threshold value that defines a certain habitat. If so, reserves may be analogous to islands in which domains above a threshold elevation are islands model assessments of nature reserves 29 10 r k!' 10 10 10' A (Unit area) 10' probability that a randomly chosen contour loop has a length s) as N(s)ccs~\ where t = 1 + (2 — H)/D and D is the fractal dimension of the contours. The cumulative distribution (the number of contours with length greater than s) is the integral of the non-cumulative distribution, JV( > s) cc S"<2-H>/D. Since the length of a contour is related to the area it encloses by s cc ADjZ (by definition), the cumulative distribution of areas enclosed by contours is iV( > ^4) cc ^4~(2_ff)/2. For a surface with /? = 2, H = j, and we have JV( > A) cc A~314, nearly consistent with the size distribution of Fig. 4(a). 102 1()3 10* L05 106 A (km2) Fig. 4. (a) Cumulative frequency-size distribution of connected areas, the number of connected areas greater than or equal to A, of the model self-similar geometry of Fig. 3. A power-law distribution JV( > A) cc A~°'s is obtained, (b) Cumulative frequency-size distribution of connected managed areas of the conterminous United States (Fig. 2). Areas which share a border are joined together to make a single larger area for this distribution. in an archipelago. If the spatial variability in one or more of those environmental factors can be modeled as a self-affine fractal, then nature reserves may be expected to have a power-law distribution of sizes similar to islands. Kondev & Henley (1995) have related the length distribution of contour lengths of self-affine fractal surfaces to the Hausdorff dimension H. The Hausdorff dimension is related to the power spectral exponent of the fractal surface by the relation /? = 2H + 1 (Mandelbrot, 1983; Hastings & Sugihara, 1993). Kondev & Henley (1995) have shown that the size distribution of contour lengths (the 4. Modeling Migration from or Between Reserves We will consider two types of dispersal in our model calculations: random-walk (short-range) dispersal and long-range dispersal. In random-walk dispersal, each individual moves to one of the four adjacent gridpoints with an equal probability during each time step. Random walk dispersal is diffusive. Random-walk movements of individuals in a species are similar to the random Brownian motion of molecules in a liquid that drives diffusion in liquids. Diffusion provides an accurate model of species migration in many cases (Andow et al, 1990). For many floral species, however, the diffusion approximation breaks down since many plants transport their seeds over long distances. Thus, in order to model plant dispersal properly seeds must be transported according to a "fat-tailed" (slower than exponential) distribution (Clark et al, 1998). In our model with long-range dispersal seeds are distributed over distance r from the parent according to an inverse power-law distribution f(r) cc r_1, where/(r) is the frequency of dispersal over distance r. The two principal features of our model, dispersal and random variations in population growth rate, appear to be necessary features in order to generate the "patchiness" often observed in the spatial distribution of abundance in a wide variety of species (Hastings et al, 1982; Levin & Paine, 1974; Roughgarden, 1977; Fasham, 1978). Each simulation begins with only one individual on each protected cell. The results of our modeling studies are summarized in Table 1. The 30 J. d. pelletier Table 1 Summary of simulation results Migration Survivability Stochasticity Dispersal Fig. 1 (a) Fig. 1 (b) Fig. 1 (c) (a) No N/A Environmental Short 20 100 30 (b) No N/A Demographic Short 10 100 15 (c) No N/A Catastrophic Short 35 40 38 (d) Yes Low Environmental Short 45 15 30 (e) Yes High Environmental Short 80 65 100 (f) N/A Low Environmental Long 65 70 65 Numbers are the average species longevity for 100 simulations, (a)-(c) indicate that many small reserve fragments have the optimal geometry for all types of stochasticity when no migration is allowed, (d) Indicates that one large reserve is optimal when migration is allowed and the survivability is low in non-protected areas, (e) Shows that the intermediate case of a self-similar distribution is optimal when migration is allowed and the survivability in non-protected areas is high, (f) Shows that species longevity is nearly independent of geometry with the long-range dispersal of seeds appropriate for many floral species. The uncertainty was estimated to be approximately 5% for (a)-(f) from the difference in the average of the first 50 and last 50 simulations. numbers represent species longevity averaged for 100 simulations for each of the three reserve geometries in Fig. 1. The numbers in Table 1 can only be compared between different geometries for a single type of simulation only because different types of simulations are parameterized differently and are not comparable. For each model, parameters of population growth, dispersal, and disturbance have been chosen so that the species longevity was in the order of 100 time steps so that individuals were able to migrate for some distance and have their histories influenced by the particular geometry of the simulation. 5. Results The effects of different types of stochasticity are illustrated in the first three simulations [Table l(a)-(c)]. For simplicity, migration from or between reserves is not allowed for these simulations. Individuals undergo a random walk but can only move to a randomly selected neighboring cell if the cell is part of a protected area. Otherwise the individual does not move during that time step. For the case of catastrophic stochasticity, we chose to increase the population with a constant rate and wipe out all individuals at a random location on the grid within a square area governed by a power-law distribution f(A)ccA~1, where f(A) is the frequency of a natural hazard with area A. One hazard occurs per time step. A power-law distribution was chosen since many natural hazards including forest fires (Malamud et ah, 1998), earthquakes (Gutenberg & Richter, 1954), landslides (Pelletier et al, 1997) and floods (Turcotte, 1994) obey power-law frequency-size distributions. What is similar about Table 1(a)-(c) is that the reserve geometry with many small reserves [Fig. 1(b)] is clearly optimal. For the case of catastrophic stochasticity the average species longevity differs less between geometries than for environmental and demographic stochasticity. In this case the growth dynamics play no role at all. The average species longevity is determined by the average return time of the hazard which affects the entire nature reserve. The return time is inversely proportional to frequency and so will be greatest for the reserve geometry which covers the most total area (including adjacent unprotected areas): Fig. 1(b). For demographic and environmental stochasticity the results indicate that there is a significant benefit to diversifying the reserve into many fragments which individually have a high rate of local extinction but collectively enable a species to persist for a long time. When random-walk migration from or between reserve fragments is allowed, the survivability of individuals in non-protected areas becomes an important parameter. It should be emphasized that we use the phrase "from or between" reserve fragments because when migration out of reserve fragments is allowed in our model, individuals will diffuse out of protected model assessments of nature reserves 31 areas even when there is only a single large reserve and no other reserve fragments to migrate into. Thus, in the case when migration is allowed and the geometry has only a single large reserve, migration will occur from the reserve but not between reserves since there is only one reserve. When the survivability is low (50% survival probability per time step), such that migration between reserve fragments is rarely or never successful, the optimal geometry is a single large square [Table 1(d)]. Note that this case is different from the case of no migration in Table 1(a)-(c) in which individuals could not leave their original reserve fragments. Individuals in Table 1(d) may leave a reserve fragment although, because of the low survivability, few will reach other reserve fragment. In this case, minimizing the total boundary between protected and non-protected areas is the best option since the flux of individuals into unprotected areas, which is proportional to the boundary length, should be minimized to maximize species longevity. The total length of the boundary is minimized for one large reserve [Fig. 1(a)]. When the survivability is higher (80% survival probability per unit time step) the intermediate case of the self-similar distribution is optimal. In this case migration between reserves is often successful and the loss of some unsuccessful colonists in non-protected areas is offset by the benefit of diversification into several reserves. Lastly, we considered long-range dispersal. In Table l(a)-(e), all new individuals originated at the same grid point as their parent. In Table 1(f), however, individuals are dispersed some distance away according to a distribution inversely dependent on distance. In this case we found, for a wide range of survivability in protected and non-protected areas, that geometry made no difference to within the 5% uncertainty in the average longevity computed with 100 simulations. It should be noted that our model does not take into account any habitat deterioration in protected areas as a result of environmentally destructive activity in adjacent non-protected areas. If activities in non-protected areas can influence the habitat quality in protected areas (Bierregaard, 1992), small isolated reserve fragments may be the least desirable reserve geometry since they maximize the total length of the edge zone between protected and non-protected areas. In addition, when reserve fragments become too small and the migration between them is low, inbreeding depression and loss of genetic diversity (Soule & Simberloff, 1986) can make the populations within reserve fragments unviable. Both of these factors are not included in the simulation model but may be very important in realistic reserve design. 6. Conclusions These simulations highlight the importance of the presence or absence of migration between reserve fragments and the survivability in nonprotected areas as being critical in determining the theoretically optimal design of nature reserves for species longevity. In the absence of migration, many small reserves are clearly optimal in this model system. If migration is allowed and the survivability is low then a single large reserve is best. If survivability is high then the self-similar distribution similar to the actual distribution of managed areas in the conterminous United States is optimal. For long-range dispersal of the seeds of plants there is very little difference between the three geometries. REFERENCES Abele, L. G. & Connor, E. F. (1979). Application of island biogeography theory to reserve design: making the right decision for the wrong reasons. In: Proceedings of the Conference on Scientific Research in the National Parks, (Linn, R. M., ed.), 1st Edn. Vol. 1, pp. 89-94. Washington, DC: National Park Service, U.S. Dep. Interior. Andow, D. A., Kareiva, P. M., Levin, S. A. & Okubo, A. (1990). Spread of invading organisms. Landscape Ecology 4, 177-188. Bierregaard Jr R. O. (1992). The biological dynamics of tropical rainforest fragments: a perspective comparison of fragments and continuous forest. Bioscience 42, 859-866. Boulinier, T., Nichols, J. D., Hines, J. E., Salter, J. R., Flather, C. H. & Pollock, K. H. (1998). Higher temporal variability of forest breeding bird communities in fragmented landscapes. Proc. Nat. Acad. Sci. U.S.A. 95, 7497-7501. Burkey, T. V. (1989). Extinction in nature reserves: the effect of fragmentation and the importance of migration between reserve fragments. Oikos 55, 75-81. Clark, J. S., Fastie, C, Hurtt, G., Jackson, S. T., Johnson, C, King, G. A., Lewis, M., Lynch, J., Pacala, S., Prentice, C, Schupp, E. W., Webb, T. & Wyckoff, P. (1998). Reid's paradox of rapid plant migration. Bioscience 48, 13-24. 32 J. d. pelletier Clemens, M. A., ReVelle, C. S. & Williams, J. C. (1999). Reserve design for species preservation. Eur. J. Oper. Res. 112, 273-283. Diamond, J. M. & May, R. M. (1976). Island biogeography and the design of natural reserves. In: Theoretical Ecology: Principles and Applications (May, R. M., ed.), pp. 163-186. Oxford: Blackwell. Fasham, M. J. R. (1978). The application of some stochastic processes to the study of plankton patchiness. In: Spatial Patterns in Plankton Communities (Steele, J. H. ed.), pp. 131-156. New York: Plenum Press. Gardner, R. H. (1987). Neutral models for the analysis of broad-scale landscape patterns. Landscape Ecology 1, 19-28. Gilpin, M. E. & Diamond, J. M. (1980) Subdivision of nature reserves and the maintenance of species diversity. Nature 285, 567-568. Glazier, D. S. (1986). Temporal variability of abundance and the distribution of species. Oikos 47, 309-314. Gonzalez, A., Lawton, J. H., Gilbert, F. S., Blackburn, T. M. & Evans-Freke, I. (1998). Metapopulation dynamics, abundance, and distribution in a microecosystem. Science 281, 2045-2047. Goodman, D. (1987). Consideration of stochastic demography in the design and management of biological reserves. Nat. Resour. Model. 1, 205-234. Gutenberg, B. & Richter, C. F. (1954). Seismicity of the Earth and Associated Phenomena. Princeton: Princeton University Press. Hastings, H. M., Pekelney, R., Monticciolo, R., Kannon, D. V. & Delmonte, D. (1982). Time scales, persistence, and patchiness. Biosystems 15, 281-289. Hastings, H. M. & Sugihara, G. (1993). Fractals: A User's Guide for the Natural Sciences. New York: Oxford University Press. Kondev, J. & Henley, C. L. (1995). Geometrical exponents of contour loops on random Gaussian surfaces. Phys. Rev. Lett. 74, 4580-4583. Krummel, J. R., Gardner, R. H., Sugihara, G., O'Neill, R. V. & Coleman, P. R. (1987). Landscape pattern in a disturbed environment. Oikos 48, 321-324. Kunin, W. E. (1998). Extrapolating species abundance across spatial scales. Science 281, 1513-1515. Levin, S. A. & Paine, R. T. (1974). Disturbance, patch formation, and community structure. Proc. Nat. Acad. Sci. U.S.A. 71, 2744-2747. Malamud, B. D., Morein, G. & Turcotte, D. L. (1998). Forest fires: an example of self-organized critical behavior. Science 281, 1840-1841. mandelbrot, B. (1983). The Fractal Geometry of Nature. New York: W. H. Freeman & Co. McGhie, R. G. (1996). Creation of a Comprehensive Managed Areas Spatial Database for the Conterminous United States: Summary Project Technical Report (NASA-NAGW-1743). NCGIA Technical Report 96-4, NCGIA, University of California, Santa Barbara, CA. McKinney, M. L. (1997) Extinction vulnerability and selectivity: combining ecological and paleontological views. Ann. Rev. Ecol. Syst. 28, 495-516. Milne, B. T. (1991). Lessons from applying fractal models to landscape patterns. In: Quantitative Methods in Landscape Ecology (Turner, M. G. & Gardner, R. H., eds), New York: Springer. Pelletier, J. D., Malamud, B. D., Blodgett, T. & Turcotte, D. L. (1997). Scale-invariance of soil-moisture variability and its implications for the frequency-size distribution of landslides. Engineering Geology 48, 255-268. Pressey, R. L., Possingham, H. P., Logan, V. S., Day, J. R. & Williams, P. H. (1999). Effects of data characteristics on the results of reserve selection algorithms. J. Biogeography 26, 179-191. Roughgarden, J. D. (1977). Patchiness in the spatial distribution of a population caused by stochastic fluctuations in resources. Oikos 29, 52-59. Sayles, R. S. & Thomas, T. R. (1978). Surface topography as a non-stationary random process. Nature 271, 431-434. Schoener, T. W. & Spiller, D. A. (1987). High population persistence in a system with high turnover. Nature 330, 475-477. Simberloff, D. & Abele, L. G. (1982). Refuge design and island biogeographic theory: effects of fragmentation. Am. Nat. 120, 41-50. Soule, M. E. & Simberloff, D. (1986). What do genetics and ecology tell us about the design of nature reserves? Biol. Conserv. 35, 19-40. Sousa, W. P. (1984). The role of disturbance in natural communities. Ann. Rev. Ecol. Syst. 15, 353-391. Turcotte, D. L. (1994). Fractal theory and the estimation of extreme floods. J. Res. Nat. Inst. Stan. Technol. 99, 377-389. turcotte, D. L. (1997). Fractals and Chaos in Geology and Geophysics, 2nd Edn. New York: Cambridge University Press. Wilson, E. O. & Willis, E. O. (1975). Applied biogeography. In: Ecology and Evolution of Communities (Cody, M. L. & Diamond, J. M., eds), pp. 523-534. Cambridge, MA: Harvard University Press. Wright, S. J. & Hubbell, S. P. (1983). Stochastic extinction and reserve size: a focal species approach. Oikos 41, 466-476.