Willow ptarmigan Lagopus lagopus is considered a popular small game species by many hunters in Scandinavia. A simple harvest strategy would be to prohibit harvest in parts of the total area. We used a spatial model of a fluctuating population of willow ptarmigan, divided into 25 subareas to investigate the possible advantages of buffer zones in managing harvest. We let the breeding success be the source of annual environmental stochasticity but without any spatial variation. Survival was assumed to be density dependent over the total area, whereas dispersal was modelled as density independent. We then compared four major scenarios in which we let dispersal and harvest vary. About 75% of the area could be left open to hunting even if the level of harvest was close to the extinction level if executed in all grids. This harvest strategy would be particularly advantageous if the goal is to provide as many hunting opportunities as possible, rather than to harvest a maximum sustainable yield. Furthermore, it is quite simple and does not need a resource-demanding control system. We believe that a harvest strategy which sets aside a part of the area as a buffer, and places a limit to the harvest effort in the grids that are open for hunting, would be a cost-efficient system with only a small risk of overharvesting.
In 1993, more than 60,000 km2 of the state-owned Swedish mountain range were opened to the public for small game hunting, with willow ptarmigan Lagopus lagopus being the most popular game species. This change generated a discussion as to whether this would result in a low risk of overharvesting, simple access for hunters as well as minimised disturbance for reindeer husbandry and other forms of land use. Wildlife managers in the region have no scientific system for deciding when there is a risk of overharvesting and thus a need to control harvest, but at present judge this from either total hunter effort and/or game bagged in an area.
Traditionally, harvesting small game through sport hunting has been assumed to have little effect on harvested populations (Allen 1954, Ellison 1991). It was believed that compensatory mechanisms and/or low hunter efficiency guaranteed a minimal risk of over-harvesting The understanding of which mechanisms that can compensate for hunting mortality in grouse is, however, poor (Ellison 1991 and references therein). Recent findings suggest that the compensatory forces may be overestimated and that hunters can locally remove a large part (>40%) of the autumn population (Kastdalen 1992, Smith & Willebrand 1999, T. Willebrand, unpubl. data).
Constant, threshold and proportional harvest strategies have been theoretically evaluated by Lande, Engen & Sæther (1995). Using a threshold with a proportional harvest above this threshold performed best in most cases. Aanes, Engen, Sæther, Willebrand & Marcström (in press) suggested an alternative approach as the best strategy when evaluating a long-term time series on willow ptarmigan population dynamics and harvest. The best strategy was to allow a proportional harvest up to a certain bag and then close the area, although the difference from a threshold strategy in combination with proportional harvest was small.
Our model assumes a closed population without any emigration or immigration. Migration has, however, been shown to strongly affect the population dynamics of subdivided populations (Ranta, Kaitala & Lundberg 1998, Björnstad, Ims & Lambin 1999, Kendall, Björnstad, Bascompte, Keitt & Fagan 2000). Spatial processes should therefore be included in realistic harvesting models (Jonzén, Lundberg & Gärdmark 2001). Accordingly, Cary, Small & Rusch (1992) showed that a local overharvest of ruffed grouse Bonasa umbellus was compensated for by immigration from large surrounding areas where limited hunting took place. They also stressed the need to consider spatial scale and dispersal between subpopulations in harvest models. McCullough (1996) suggested that a maximum sustainable yield (MSY) could be achieved by closing a part of an area to harvesting, but also pointed out that the effects were dependent on the assumptions about the dispersal rate between hunted and closed areas (see also Joshi & Gadgil 1991).
In this paper, we aim to evaluate the effects of closing parts of an area, i.e. the use of buffer zones, on reducing the risk of overharvesting willow ptarmigan. Such a simple management strategy would require less monitoring and still provide unregulated access to the open parts.
Methods
Densities
Spring densities of willow ptarmigan can range within 2–20 pairs/km2 (Marcström & Höglund 1980, Brittas 1988). Similar spring densities have been reported in Norway (Steen 1989). The highest recording of spring numbers was made on Tranöy, averaging more than 70 pairs/km2 during 1960–1980 (Myrberget 1988). Willow ptarmigan seem to fluctuate in synchrony over large areas, but regional synchrony is less pronounced than that of snowshoe hares Lepus americanus in North America. The synchronising agent for willow ptarmigan seems to be the vole cycle which causes shifts in the level of predation on eggs and chicks (Myrberget 1988).
Breeding success
Willow ptarmigan produce one brood per year and the annual variation in production of young is large. Several long-term studies show variation ranging within 0.5–6 young/pair in late summer/early fall (Marcström & Høglund 1980, Myrberget 1988). There are no indications that the production of young is density-dependent in Scandinavia. Myrberget (1988) showed that even if the population at Tranöy was reduced by 50% there was no change in production of young per female. Data from mid-Sweden also show a lack of correlation between density and the production of young (V. Marcström, pers. comm.). Furthermore, Ellison (1991 and references therein) concluded in his review of harvesting tetraonids that although breeding has been extensively studied, compensatory natality has rarely been detected.
Survival
Predation is the major cause of non-hunting mortality in willow ptarmigan (Marcström, Kenward & Engen 1988, Brittas & Willebrand 1992, Smith & Willebrand 1999) accounting for 70–100% of natural mortality depending on age and season. Despite the relatively consistent proportion of deaths caused by predation, annual survival rates of grouse range within 0.1 -0.7 (Choate 1963, Jenkins, Watson & Miller 1963, Zwickel & Bendell 1972, Martin & Hannon 1987). Annual survival estimates for willow ptarmigan in Scandinavia are mostly based on mark/recapture estimates. On Tranöy the apparent annual adult and young survival was 50 and 27%, respectively (Myrberget 1988, Steen & Erikstad 1996). In a large study using radio-tagged individuals, the estimate of annual survival was close to 40%; no difference in survival between adults and young of the year (measured from the first autumn) was found (Smith & Willebrand 1999). Using 2.8 young/pair as a longterm production average, a stable population should have an average annual survival of 42%, assuming no environmental stochasticity.
If there is no strong density dependence in the production of young in willow ptarmigan, increased mortality caused by harvest has to be compensated for either by reduced natural mortality or, locally, by increased immigration. Increased natural survival is often assumed to be the possible responsible compensatory agent in harvested local populations (Myrberget 1985, Kastdalen 1992, Mortensen 1994), and increased immigration will eventually increase the overall survival. Myrberget (1985) concluded that hunting had not resulted in any long-term declines in Norway, and even areas with harvest rates of 30–50% did not show a decline in population density (Kastdalen 1992).
Movement and dispersal
Dispersal of yearlings has been suggested as the principle mechanism by which some tetraonid populations are controlled at high population densities (Boag & Schroeder 1987, Moss & Watson 1990, Hudson 1992). In Sweden, there is no seasonal migration of willow ptarmigan. The great majority of young females move away from their natal site, whereas young males do not (mean natal dispersal 11.4 km. vs 2.6 km; Smith 1997). In a recent study by Smith & Willebrand (1999) on radio-marked willow ptarmigan, no net effect of movement into a harvested area from the surrounding non-hunted areas was observed. Movements and dispersal of adults and young males could redistribute birds at most on just a very local scale (distances of <5 km). Only dispersing young females (and possibly some adult females) could be responsible for the exchange of grouse between populations at larger scales: young female willow ptarmigan appear to be obligate dispersers at all densities (Smith 1997).
The harvest
Willow ptarmigan in Scandinavia are usually managed in subunits of 10–100 km2 depending on the type of land ownership. Harvest regulation is mostly imposed by season length, but limits on bag or effort are sometimes used. Limits are rarely based on what is considered normal from historical records. Harvest rates of grouse are usually in the range of 5–10%, but rates up to 50% of the autumn population have been reported (Kastdalen 1992, Smith & Willebrand 1999). Smith & Willebrand (1999) suggested that hunting mortality was almost completely additive to natural mortality, and that immigration was the compensatory force in their study of willow ptarmigan harvest in state-managed areas. However, the adjacent areas that were not hunted did not appear to be the major source of immigration, and they suggested that immigrants had to come from a large area.
Model structure
The population
We placed our model population of willow ptarmigan in a rectangular world represented by 5 x 5 km grids that could be linked by dispersal (Fig. 1). The area is considered closed: edge grids were projected to the grids on the opposite edge (shaded grids in Figure 1). Each grid was assumed to have an average population of 1,000 individuals with an equal sex ratio and a ratio of 58% young. The model kept track of the sex and age ratio in each grid (a spatial Birth-Immigration-Death-Emigration-model (BIDE-model); Renshawl991). The population was evaluated once a year just before the opening of the hunting season. Hunting was assumed to be instantaneous and thereafter the density-dependent survival rate was determined. We assumed that natural mortality took place before dispersal. We could either let the harvest be at the same level in all 25 grids (no spatial stratification), or we could let a different number of grids be harvested at some predetermined level. We always started the harvest in grid 0 and continued with 1, then 2 and so forth.
Breeding success
Demographic stochasticity was introduced into the model by a random annual variation in chick production. We used a normal distribution with a mean of 2.8 chicks/pair and a standard deviation of 0.5. We also adjusted the mean by ±0.4 following a sinusoidal wave of 3.5 years to mimic a cyclic component. All values less than 0.2 were set to 0.2 chicks/adult. An even sex ratio of chicks was assumed.
Survival
Survival was the only part of the BIDE-model that was density dependent. We assumed that chick production and dispersal have negligible density dependence when compared with survival. Two earlier studies suggested that density-dependent mechanisms act on a larger scale than normally assumed (Small, Holzwart & Rusch 1991, Smith & Willebrand 1999), and that estimates of density-dependent survival are masked by movement on the local level. At the very local level, however, these are not necessarily density dependent.
The annual survival of willow ptarmigan was as high as 80% in one of our study areas with radio-marked birds after densities had declined for several years (T. Willebrand, unpubl. data). We assumed that a brood size of 2.8 would keep the population at equilibrium level, corresponding to a survival rate of 0.4167. Thus, we had two data points and could estimate a two-parameter relationship with density and survival, but we have limited information on which relationship to use.
We used Hollings Disk equation,
because this makes it possible to state the maximum survival at low densities (a), in our case 0.85; in the equation, S is survival, m and a constants, Nx the population in grid X, and Hx the bag size in grid X. All grids experienced the same survival (%) which was calculated using the total population in all 25 grids. Selecting the form and strength of the relationship is crucial for estimating the size of MS Y and its corresponding population level (Willebrand 1997). However, as our aim is to make statements of buffer zones, we use the same relationship in all the scenarios we compare.
Dispersal
Most vertebrate species show dispersal patterns that are related to density (Hansson 1991). Such density-dependent processes strongly influence the local as well as the global dynamics of spatially subdivided populations (Lande, Engen & Sæther 1999, Sæther, Engen & Lande 1999). This also seems to be true for many harvested game species (Beasom & Robertson 1985). However, we chose to explore the situation where this effect was negligible compared to density-dependent survival because we have no data to support any relationship between density on one hand and emigration/immigration on the other. There is no difference in survival between residents and dispersing grouse in our model, and natural mortality takes place before dispersal starts. We have three different levels of dispersal in our model; no dispersal, young are divided equally over all 25 grids, and only young females disperse to the four closest grids.
Scenarios
We wanted to investigate the outcome under the following scenarios:
A) Global harvest/no movement: these simulations were run to establish a classical MSY relationship (Sutherland 2001) between population level and bag size. Three levels of harvest based on these results where later used in the two simulations below. One level was just below MSY, the second and third levels were substantially higher. The second level was as high as possible without making the population going extinct, and the third level was a harvest that would always lead to extinction.
B) Local harvest/no movement: in these simulations we did not allow any movement between grid cells despite a density-dependent survival that was based on the total population. In the simulation we arbitrarily chose a harvest of 10%. We started with grid cell 0 and followed the numbering in Figure 1.
C) Local harvest/young disperse: here we used the three harvest levels selected in Scenario A. All young remaining after harvest were pooled and divided equally among the 25 grid cells. As in Scenario B, we increased the number of grid cells that were harvested following the numbering in Figure 1.
D) Local harvest/restricted dispersal: the same three harvest levels as in Scenarios A and C were used in this scenario, but here we put a strong restriction on how young grouse could move between grid cells. Only young females could move between cells, and movement made each cell receive one fourth from each of the four closest cells. All young females that were born in a grid and were not harvested, emigrated out of the grid. As in the two previous scenarios, we increased the number of grids that were harvested in each simulation.
We ran the simulation for 1,500 years and extracted the last 1,000 values of population and bag levels in each simulation. We compared the local population and bag in the grids 0, 5,10 and 20 in addition to the total population and the total bag for all 25 grids. In Scenario D we also extracted the sex ratio. The mean and standard deviation were calculated for all the extracted variables in each simulation. The model was developed in the simulation software POWERSIM version 2.5C (Powersim 1996), and the SAS statistical package (version 6.12, SAS 1989) was used for all statistical calculations.
Results
Scenario A - Global harvest/no movement
Maximum Sustainable Yield (MSY) was achieved at 24% harvest, and the bag was 2,948 (SD = 863) at a population size of 12,182 (SD = 3,920). Harvesting the population at 48% was still possible (N = 632, SD = 316, H = 291, SD = 159), but increasing the harvest rate to 50% made the population rapidly go extinct. Survival cannot increase to more than 0.85 in our model. Because we assume that compensation only occurs through reduced natural mortality, there is a limit to how much harvest the population can withstand. A deterministic model predicts extinction when the harvest is larger than 50.9%. We chose to use 22,48 and 55% harvest in the later simulations (Scenarios C & D).
Scenario B - Local harvest/no movement
Harvested grid populations went extinct when some grids were not harvested in a system where no dispersal was allowed. This was true even if the harvest was as low as 10%. Survival was determined by the total population size in all the 25 grid cells, but the local reduction due to harvest will have a small effect on the total population change. The increase in local survival was insufficient to compensate for the harvest. The harvested grid(s) entered a path to extinction whereas populations on the grids that were not hunted would increase to larger and larger values as more and more grid cells were harvested. In the extreme case, one remaining grid will host a population that previously was resident in all the grids, a 25-fold increase! The system could not possibly provide a sustainable harvest.
Scenario C - Local harvest/young disperse
Introducing dispersal between the subpopulations changed the results dramatically. We assumed that the production of young remaining after the harvest is equally shared among all grid cells. This made it possible to harvest 23 of 25 grids at 55%. This was a harvest level that would make the global population rapidly go extinct in Scenario A. At 55% the MSY was 2,987 (SD = 1,094) at a population size of 12,763 (SD = 3,705) and 16 grids harvested (Fig. 2). Reducing the harvest rate to 48% will change the corresponding values to 2,993 (SD = 1,063), 11,785 (SD = 3,543) and 18 grids. At a more modest harvest of 22% all grids were harvested without reaching a peak since the harvest level was below the MSY of the global population. Harvesting all grids at 22% gave a bag of 2,849 (SD = 902) at a population size of 12,945 (SD = 4,099). In this scenario the population consisted of two categories, either a non-harvested grid, or a harvested grid with a reduced density compared with the grids that were not harvested. As more grids were harvested both types declined. The maximum bag from a grid was obtained when it entered the subset of harvested grids and then it declined as more grids entered. The yield was the combination of the number of grids harvested and the bag in each grid.
Scenario D - Local harvest/restricted dispersal
The results were qualitatively similar when only young females were allowed to move among neighbouring grids. In this situation, it was possible to harvest 24 of 25 grids at 55%. At 55% harvest, the MSY was 2,365 (SD = 976) at a total population size of 10,882 (SD = 3,288) and 18 grids harvested. Reducing the harvest rate to 48% changed the values of MSY to 2,492 (SD = 962), 12,078 (SD = 3,536) and 18 grids. The lower harvest of 22% gave a bag of 2,936 (SD = 877) at a population size of 13,340 (SD = 3,987). A considerable difference in behaviour compared to Scenario C was that grids could decline in numbers but then start to increase again depending on the harvest in the surrounding grids. This pattern was also true for bag size.
The overall sex ratio showed a slow decline (Fig. 3) as more and more grids were harvested. The sex ratio declined to 0.68 females/male at a 48% harvest when 24 of 25 grids were harvested. However, the most pronounced effect of this scenario was the drastic increase in female to male ratio in harvested grids. The harvested grids show a high female to male ratio due to a net influx of young females from neighbouring grids with higher densities, especially if they were not harvested. The overall decline in sex ratio is explained by an increased survival of grouse in grids that are not harvested but which experience a negative balance in the exchange of young females. Thus, there were grids that were harvested which showed a high female to male ratio, whereas the grids that were not harvested had a higher density in combination with a female to male ratio that was less than one.
Figure 3.
Change in female to male ratio in the total population (o) and in grid 10 (A) as the number of harvested grids increased in a system where only young females disperse. All values are presented as the average and standard deviation of the last 1,000 measurements of 1,500 years' simulation. The harvest level was chosen as 48% to show the return to equal sex ratio when all grids are harvested.
Discussion
Our analysis suggests that the use of buffer zones may provide a simple strategy for managing the harvest of our stochastic model population. It is particularly advantageous when, rather than to harvest a maximum sustainable yield, the goal is to provide as many hunting opportunities as possible without the need for a resource-demanding control system. This management system would only need to ensure that the refuges are not utilised once the needed size and distribution of buffers are established. Additional control over the system could be obtained by monitoring harvest effort and age ratio in the bag without consuming large amounts of resources. The harvest effort could be used to estimate harvest rate since there is a linear relationship up to about 30% harvest (Kastdalen 1992, T. Willebrand & M. Hörnell, unpubl. data). This relationship is not greatly affected by variation in grouse density and reproductive success (T. Willdebrand & M. Hörnell, unpubl. data). About 75% of the area could be open to hunting even if the harvest was close to the extinction harvest level when executed in all grids. We believe that a harvest strategy which sets aside part of the area as a buffer and places a limit to the harvest effort in the grids that are open for hunting would be a cost-efficient system with only a small risk of overharvesting. However, to date, probably very few areas in Sweden open to willow ptarmigan harvest have experienced harvest rates larger than 35% (T. Willebrand & M. Hörnell, unpubl. data), so there is no immediate need to change the present system. Parts of the area are already closed for other reasons, and we are at present evaluating the size and distribution of these areas in relation to willow ptarmigan habitat.
Maximising the total bag would not give the highest bag for the grids that are harvested. The grid that is harvested intensively and surrounded by other harvested grids will have a low density and provide a small bag, especially when dispersal is restricted. This could be the case when a road system or lodging facilities create an access point, as shown by Brøseth & Pedersen (2000). However, we anticipate that hunters would adjust their effort in relation to grouse encounter rate (positively) and distance from roads (negatively). Access to, or effort in the open grids would probably not have to be regulated as to optimise harvest effort in the open area, and monitoring the effort is likely sufficient.
This harvest strategy would not work if ownership was divided so that each landowner had a single grid to manage. A grid in our system (in Scenarios C or D) could not be managed on its own. The population in a grid could change rapidly even if there is no change in hunting practices if the harvest in the surrounding areas changed. Our method would work best in settings where, for example, the government has control over large areas. A larger yield would probably result from following the strategies suggested by Sæther, Engen & Lande (1996), but would also require greater knowledge about the harvested population. Such information could be expensive and difficult to obtain. Furthermore, it would assume the control over such a large area that the surrounding areas would not affect its dynamics.
McCullough (1996) encouraged the development of spatial approach with buffer areas similar to our grid system to avoid the accelerating overharvest that is characteristic of quota harvest systems when a harvest error is not discovered quickly. However, he concluded that true metapopulations have little potential for harvest, and that it would be problematic to design a simple and robust harvest model. On the other hand, the advantage of buffer areas when developing harvest strategies for wildlife management has not been widely adopted in research. In fisheries, reserves or refuges have been proposed to reduce the risks of overharvesting (Roberts 1997), which has also received some support under certain conditions from analyses of general harvest models (see Jonzén et al. 2001). The use of buffer areas in game harvesting management has a close relation to the source-sink concept (Pulliam 1988). A local harvest creates an artificial sink that makes the buffer area develop into a source. Because we assume that the density-dependent processes operate on a larger scale, this would be a true sink and not a pseudo-sink (Watkinson & Sutherland 1995, Dias 1996).
Our conclusions are only valid as long as our model sufficiently correctly describes the processes in the harvested population. We did not include any density dependence in dispersal movements and chose a relationship between density and survival that had a large-scale effect despite a local reduction. This assumption is supported by the results of the analysis of long-term data on the population dynamics of willow ptarmigan (T. Willebrand, V. Marcström, R. Brittas & M. Hörnell, unpubl. data). We have only observed high survival (>80%) in an area with a large-scale decline in the population due to poor breeding (T. Willebrand, V. Marcström, R. Brittas & M. Hörnell, unpubl. data), where as we do not see any significant increase in survival in response to a local reduction by harvest. However, it is not realistic to have the same survival (%) in neighbouring grids with extreme differences in density. Including a mechanism that would allow predators to remove a larger proportion of grouse in grids with higher density would strengthen the resilience to harvest.
A higher capacity to withstand harvesting would also be the case if we included a density-dependent dispersal between grids. Our models give a net influx of grouse into harvested grids since dispersal was the mediator of density-dependent survival in the population. This effect was achieved by letting the grids share their post-harvest production either by dividing the young equally among grids or exchanging it with the closest neighbouring grids. Thus, grids with few young would benefit from areas with more young. In several grouse studies, resident grouse (adults) rarely move to new areas to breed. Some adult females show a seasonal migration that probably retrace their natal dispersal from their first wintering site since other females will remain in the area they leave, but only young females seem to move into new areas. The dispersal in Scenario D is what we believe to be the most realistic in our models, although the large bias in the sex ratio of harvested grids would require a high harvest for at least 30 years. We also suggest a field experiment where an area is depleted of grouse to test if the attraction to dispersing females is lost, i.e. if there is a lower threshold where immigration would be inversely density dependent.
Many of the processes in our model act on such a large scale that it will be difficult to study them in traditional experimental research. Evaluating competing models by comparing with available data as suggested by Hilborn & Mangel (1997) is probably a more fruitful approach, but it would require the assistance of managers willing to adopt harvest strategies that provide data to compare different models (Walters 1986). This work is an early attempt to find appropriate these competing models.
Acknowledgements
we are grateful for help provided by the county officials in Jämtland, Västerbotten and Norrbotten. Dr. J. Ball, prof. B-E. Sæther and an anonymous referee provided helpful comments on the manuscript. The MISTRA-fund, the Swedish Environmental Board and the Swedish Hunters' Association provided funding for the study.
