Schapaugh_A MESThesis2009.pdf

Media

Part of The Dynamics and Viability of the Endangered Streaked Horned Lark (Eremophila alpestris strigata)

extracted text (extracttext:extracted_text)
THE DYNAMICS AND VIABILITY OF THE
ENDANGERED STREAKED HORNED LARK
(Eremophila alpestris strigata)

by
Adam W. Schapaugh

A Thesis
Submitted in partial fulfillment
of the requirements for the degree
Master of Environmental Studies
The Evergreen State College
June 2009

© 2009 by Adam W. Schapaugh. All rights reserved.

This Thesis for the Master of Environmental Study Degree
by
Adam W. Schapaugh

Has been approved for
The Evergreen State College by

Dr. Alison Styring
Member of the Faculty

Dr. Scott Pearson
Westside Research Team Leader and Senior Scientist
Washington Department of Fish and Wildlife

Dr. Robert Knapp
Member of the Faculty

Date

ABSTRACT
THE DYNAMICS AND VIABILITY OF THE
ENDANGERED STREAKED HORNED LARK
(Eremophila alpestris strigata)
Adam W. Schapaugh
The Streaked Horned Lark (Eremophila alpestris strigata) is an endangered songbird
endemic to prairie and open coastal habitats west of the Cascade Mountains. Its
persistence is handicapped by dramatic habitat loss, human disturbance and predation.
Fewer than 1,000 individuals remain in four isolated subpopulations across Washington
and Oregon. By combining demographic data and information on size and distribution of
local populations, I developed a stage- and space-structured demographic model to
analyze E. a. strigata’s viability in Washington State. Importantly, simulations that
include variation and correlation in survival rates yield variable, yet qualitatively
consistent, forecasts of population growth. The model predicted a continuing statewide
population decline and near certain risk of extinction over the next 25 years. To
determine where conservation efforts and data collection are best focused, I performed a
series of perturbation analyses in which the effects of changing vital rates were
quantified. I found that population growth is most sensitive to the survival of adults.
Under a modest scenario, increasing the survivorship of adults by 10% was sufficient to
lengthen the median time to extinction by more than 5 years. A more optimistic
improvement (20%) yielded a viable Coastal subpopulation. In contrast, I found that
some well-meaning management activities, in particular efforts solely targeting
fecundity, are unlikely to be either cost effective or biologically sound. Although many
anthropogenic impacts threaten E. a. strigata, the subspecies’ future could be bright,
provided that research and management focus on biologically significant aspects of its
life history. However, I demonstrate that complacency is ill-advised; decisive action to
quickly improve demographic rates is needed, given the consistent qualitative output
across models and the inherent uncertainty in predictions of future population trends.

TABLE OF CONTENTS
LIST OF FIGURES................................................................................................. v
LIST OF TABLES ................................................................................................. vi
ACKNOWLEDGEMENTS ..................................................................................... vii
INTRODUCTION ................................................................................................... 1
TAXONOMY AND CONSERVATION STATUS ............................................... 2
DISTRIBUTION AND HABITAT ASSOCIATIONS .......................................... 2
POPULATION ESTIMATES AND THREATS ................................................. 4
METHODS............................................................................................................ 5
MODELING APPROACH ............................................................................ 5
MODEL STRUCTURE ................................................................................ 5
PRIMARY DATA SOURCES........................................................................ 7
POPULATION SPATIAL STRUCTURE ......................................................... 8
PARAMETER ESTIMATES ......................................................................... 9
DEMOGRAPHIC PARAMETERS............................................................ 9
CARRYING CAPACITY, INITIAL ABUNDANCES
AND STAGE DISTRIBUTIONS ............................................................. 10

DETERMINISTIC MODEL ........................................................................ 11
STOCHASTIC MODEL ............................................................................. 12
STOCHASTIC SIMULATIONS OF MANAGEMENT SCENARIOS................... 14
RESULTS ........................................................................................................... 15
DETERMINISTIC MODEL ........................................................................ 15
STOCHASTIC MODEL ............................................................................. 18
SIMULATION SCENARIOS ....................................................................... 21
DISCUSSION ....................................................................................................... 23
LITERATURE CITED .......................................................................................... 26

LIST OF FIGURES
Figure 1a-b. Approximate historical distribution and current distribution of E. a. strigata
(adapted from Gabrielson and Jewett 1940; Behle 1942; Rogers 2000; Altman 2003;
Stinson 2005). ................................................................................................................ 3
Figure 2. Simplified life cycle of E. a. strigata. Fledglings survive to become adults
according to a juvenile survival probability, P1. Birds breed at one year of age and
survive thereafter according to an adult survival probability, P2. F1 and F2 are fecundity
transition probabilities for juveniles and adults, respectively. The ‘self-loop’ on juveniles
arises from the fact that juveniles begin the year as fledglings, but by the time of the next
census, have produced fledglings of their own ................................................................. 6
Figure 3. Spatial distribution of the subpopulations making up the E. a. strigata
population in Washington. Arrows indicate infrequent migration between
subpopulations, the consequences of juvenile dispersal ................................................... 8
Figure 4. Projection in time of each of the subpopulations (Deterministic model), with
the statewide decline shown in bold .............................................................................. 15
Figure 5. Projection in time of each of the subpopulations using pooled estimates of
survival parameters. Juvenile survival in each subpopulation was set to 0.16; adult
survival was set to 0.47 ................................................................................................. 17
Figure 6. Stochastic projection in time of the statewide population, indicated by the solid
line (mean values). The dashed lines represent the upper and lower bounds of the 95%
confidence intervals around the mean estimate ............................................................. 18
Figure 7. Stochastic projection in time of each of the subpopulations, indicated by the
solid line. Again, the dashed lines represent the lower and upper bounds of the 95%
confidence intervals around each mean estimate. ........................................................... 19
Figure 8. Stochastic projection of the statewide population, indicated by the solid line,
using pooled estimates of survival from Camfield et al. 2007. The dashed lines represent
the 95% confidence intervals around the mean estimate................................................. 20
Figure 9. Cumulative extinction probability of the statewide population, according to the
initial model (solid line) and model using pooled estimates of survival taken from
Camfield et al. 2007 (dashed line). ................................................................................ 20
Figure 10. The stochastic growth rate (λs) produced by independent perturbations of
entries in the projection matrix. The baseline λ s for each subpopulation is indicated by
the solid line ................................................................................................................. 21
Figure 11. The stochastic growth rate produced by simultaneous perturbations of entries
in the projection matrix. Increases shown in the legend indicate, respectively, how pairs
of vital rates were perturbed. The baseline λ s for each subpopulation is indicated by the
solid line ....................................................................................................................... 22

LIST OF TABLES
Table 1. Vital-rates used to parameterize a stage- and space-structured model for E. a.
strigata. ....................................................................................................................... 10
Table 2. Subpopulation-specific parameters used in simulations .................................. 10
Table 3. Local deterministic growth rates (λd) and sensitivity and elasticity coefficients
of vital rates (defined in Table 1). 95% confidence limits for elasticity values are shown
in parentheses. s1(i) = annual survivorship of juveniles; s2(i) = annual survivorship of
adults; ff = annual fecundity; mba = probability of juvenile moving from site a to b. ...... 16
Table 4. Results of simulations to determine the sensitivity of local population growth to
increases in vital rates. The maximum local λd, given the uncertainty in each parameter,
is shown in boldface. .................................................................................................... 17
Table 5. Some results of the simulations (values from 10,000 replicates). ................... 18
Table 6. Sensitivity and elasticity coefficients of local stochastic growth rates (λ s) to
model parameters and variances. .................................................................................. 20

ACKNOWLEDGMENTS
Many thanks to Dr. Alison Styring (The Evergreen State College, MES program)
for serving as my primary advisor and her two years of support, excellent review, and
encouragement during my time in the MES program. I also thank Dr. Scott Pearson of
the Washington Department of Fish and Wildlife for providing the data that formed the
core of my thesis and for his excellent feedback, ideas, and support. I also thank Dr. Rob
Knapp (The Evergreen State College, MES program) for serving as a reader and for his
timely and constructive comments. Lastly, for her enduring love and encouragement, I
thank my wife, Jami.

INTRODUCTION
Globally, habitat loss and degradation has precipitated widespread declines in
many grassland birds (Askins 1995, Donovan and Flather 2002). This is true of the
Streaked Horned Lark (Eremophila alpestris strigata), a rare subspecies of Horned Lark
endemic to the Pacific coast of North America (Pearson and Altman 2005). Historically,
E. a. strigata occupied prairie and open coastal habitats west of the Cascade Mountains in
British Columbia, Washington and Oregon (Stinson 2005). However, the amount of
suitable habitat has contracted over the past century, due primarily to the intensification
of agriculture and urban development. Not coincidentally, E. a. strigata has been
extirpated from much of its historic range and is now confined to small pockets of
isolated habitat in Washington and Oregon (Stinson 2005). The global population is
likely below 1,000 individuals (Pearson and Hopey 2005); distributed among four
spatially distinct subpopulations in the south Puget Sound, lower Columbia River,
Washington coast, and Willamette Valley (Rogers 2000; Stinson 2005).
Despite documented shifts in the subspecies breeding range (e.g. Rogers 2000),
the rapid decline of E. a. strigata has gone largely unnoticed. Not discounting nearly a
decade of applied research by the Washington and Oregon Department’s of Fish and
Wildlife, long-term ecological information is limited, as have been quantitative
evaluations of the subspecies status (but see Pearson et al. 2008). Here, I present a
viability analysis of the three remaining subpopulations of E. a. strigata in Washington
State. I combined demographic data and information on size and distribution of each of
the subpopulations to construct a stage- and space-structured stochastic model, which I
used to (1) determine the likely outcome of current population processes; (2) examine
explicitly the separate influences of fecundity, adult survival, and juvenile survival on
population dynamics; and (3) predict the consequences of alternative recovery actions.
My use of stage-based models for population viability analysis differs from many
published population models in two ways. First, I pay particular attention to the
importance of temporal variability, using straightforward estimates of survival variation
from field data. Second, I use the models not only to address threats, but also to ask how
best to prioritize future data collection. Although there is more demographic data on E.
a. strigata than is often available for endangered species, significant gaps and
uncertainties do compromise some of the results. Thus, I use the models to draw
conclusions about the relative benefits of further research on different life stages. Ideally,
this work will prove to be beneficial for the conservation of E. a. strigata, but also to

evaluate the applicability of theoretical concepts (e.g. nonequilibrium dynamics) applied
to the management of endangered songbirds in fragmented landscapes.

TAXONOMY AND CONSERVATION STATUS
The Horned Lark (Eremophila alpestris) is a member of the family Alaudidae
(larks) in the order Passeriformes (Sibley 2000). It is a small ground-nesting passerine
native to both the New and Old Worlds (Stinson 2005). Across North America, there are
21 described subspecies based primarily on differences in size and plumage color
(American Ornithologists’ Union 1957). E. a. strigata was first described by Henshaw
(1884), based on specimens collected in Pierce County, Washington (Stinson 2005).
Each individual is marked by a dark breast-band, lores, and malar stripe that contrast with
the yellow to white supercilium. The sexes can be distinguished by the male’s brighter
plumage and larger body (Beason 1995). In size and appearance, E. a. strigata most
closely resembles E. a. insularis of the California Channel Islands, but is characterized by
a less heavily streaked breast and bright yellow on the abdomen and flanks (Behle 1942).
Today, E. a. strigata is a federal candidate for listing under the Endangered
Species Act by the U.S. Fish and Wildlife Service. In Canada, it is recognized by the
Committee on the Status of Endangered Species and is Red-listed in British Columbia.
In Oregon, it is considered a state-sensitive species by the Department of Fish and
Wildlife. Following the completion of a range-wide assessment and preliminary
conservation strategy by Pearson and Altman (2005), the subspecies was listed as
endangered by the Washington Department of Fish and Wildlife on March 2, 2006.

DISTRIBUTION AND HABITAT ASSOCIATIONS
Historically, E. a. strigata bred in prairie and open coastal habitats from the
southwestern corner of British Columbia, through the Puget Trough, and as far south as
Eugene, Oregon (Fig. 1a)( Stinson 2005). Early authors describe the breeding range as
including the Sierra Nevada and Channel Islands of California (Gabrielson and Jewett
1940), but Behle (1942) later assigned these groups to E. a. insularis and E. a. sierra. In
Washington, E. a. strigata was largely confined to the glacial outwash prairies of the
south Puget Sound region. However, the subspecies could be found in areas of the
northern Puget Trough, along the Pacific coast extending northward from Willapa Bay,
and on islands of the Columbia River (Jewett et al 1953). In Oregon, E. a. strigata was a
common breeder throughout the Willamette and Rogue River valleys in the western half
of the state (Stinson 2005). As recently as the 1940’s, individuals were commonly seen

in the rangelands of Linn and Benton Counties and the grasslands east of Medford
(Gabrielson and Jewett 1940).
Today, more than 90% of the original grassland habitat has been lost in the south
Puget Sound (Crawford and Hall 1997). Across Washington State, distribution surveys
conducted by Rogers (1999) and Maclaren and Cummins (2000) found only 11 nesting
sites occupied. Subsequent survey work by Pearson and Hopey (2004, 2005) indicate a
small likelihood of finding additional nesting locations, as most sites with potentially
suitable habitat have already been surveyed. In Oregon, nearly all (>99%) of the presettlement grasslands are gone (Johannessen et al. 1971; Towle 1982), thus confining the
subspecies almost entirely to the Willamette Valley (Fig 1b) (Rogers 2000). Most current
sightings come from Marion and Polk Counties, especially on and around Basket Slough
National Wildlife Refuge (Altman 2000). Other small, scattered groups may be found in
the lightly populated areas between Peoria and Harrisburg, and southeast of Portland near
Estacada (Altman 1999).

Vancouver
Island

Vancouver
Island

Puget
Sound

Puget
Sound

Pacific
Ocean

Washington

Pacific
Ocean

15 km

15 km

Columbia
River

Columbia
River

a.

Washington

Oregon

b.

Oregon

Figure 1a-b. Approximate historical distribution (a.) and current distribution (b.) of E. a.
strigata (adapted from Gabrielson and Jewett 1940; Behle 1942; Rogers 2000; Altman 2003;
Stinson 2005).

POPULATION ESTIMATES AND THREATS
It is likely impossible to establish a reliable estimate of historic population size.
As an example, it is unknown what proportion of the 170,000 acres of grassland soils in
the Puget Sound supported suitable habitat (Chappell et al. 2001), or if the proportion
changed or shifted due to dynamic abiotic factors (Stinson 2005). Today, the work of
Altman (1999), Rogers (2000), Pearson (2003) and Pearson and Hopey (2004, 2005),
indicate that there are approximately 800 E. a. strigata remaining, nearly half of which
breed in the Willamette Valley (398 birds). One-third occupy habitats in the south Puget
Sound (222 birds) while the remaining individuals are split nearly equally between the
lower Columbia River and outer coast (68 and 86 birds, respectively).
As is true for many songbirds, habitat loss has, and will continue, to compromise
the viability of E. a. strigata. Presently, grasslands west of the Cascade Mountains are
among the most endangered ecosystem types in the region (Dunn and Ewing 1997;
Rogers 1999, 2000). Across Washington and Oregon, grasslands composed
predominantly of native species occupy just 3% of their pre-settlement range (Crawford
and Hall 1997). In addition to the usual sources of habitat loss, such as urban
development, fire suppression now severely threatens open habitat by allowing the
establishment of both alien and native flora (Rogers 1999). In grasslands that do persist,
encroachment by alien species such as Agrostis tenuis increase the density of vegetation
beyond what is typically utilized for nesting (Beauschesne and Cooper 2003). This is
especially problematic in the south Puget Sound, where most of the remaining habitat is
within Fort Lewis Military Base, and has been managed under a policy of fire
suppression for decades (Altman 1999).
The coastal nesting areas have also experienced a similar, albeit less extensive,
contraction of suitable habitat. Even though there has been no systematic attempt to
estimate the degree of landscape alteration, the loss of habitat to the invasion of nonnative beachgrasses (Ammophila spp.) and erosion is likely considerable (Pearson and
Altman 2005; Stinson 2005). European beachgrass was introduced as part of dune
reclamation programs on the west coast of the United States beginning in the 19th century
(Wiedemann 1987), and has increased as much as 574% over the past fifty years along
portions of the Washington coast (Buell et al. 1995).

METHODS
MODELING APPROACH
Population Viability Analyses (PVA) are a collection of quantitative methods
used to predict the likely future status of a population or collection of populations of
conservation concern (Morris and Doak 2002). The earliest applications of PVA were
the stochastic models of Shaffer (1981, 1983), who explored the consequences of grizzly
bear (Ursus arctos) management in Yellowstone National Park. By that time,
deterministic demographic analyses had been used for nearly a decade in the management
of endangered species (Miller and Botkin 1974). Shaffer offered a new approach in
population modeling when he developed a stochastic simulation that incorporated chance
events (demographic and environmental stochasticity), producing probabilistic estimates
of extinction risk (Shaffer and Samson 1985). PVA became a heuristic concept when
Gilpin and Soule (1986) broadened it to include an examination of the synergism of
interrelated risk factors through multivariate modeling. Since that time, the idea of PVA
as a process of risk analysis has remained, where deterministic and stochastic factors are
identified, risks associated with each are considered, and a model describing both factors
and risks is developed (Beissinger and Westphal 1998).
Perhaps the most useful application of PVA, given the urgency and limited
resources for managing declining species, is an evaluation of the likely outcomes of
different management options on specific populations. Doing so provides a way to
prioritize efforts in a manner that maximizes the likelihood of population persistence.
(Pearson et al. 2008). Although some critics of PVA doubt its utility in conservation
planning, alternative methods of making conservation decisions are often poorly adapted
to account for uncertainty and less transparent about their reliability (Brook et al. 2000).

MODEL STRUCTURE
I developed a stage- and space-structured demographic model to quantitatively
investigate E. a. strigata’s viability in Washington. Only females were included in its
formulation because I found no evidence of male limitation on female reproductive rates.
I examined deterministic and stochastic versions; the latter of which allowed for year-toyear changes in survival rates due to fluctuations in the environment over time. Each
matrix was parameterized according to a post-reproductive census as described in
(McDonald and Caswell 1993): the population is censused in late summer after fledging
but before any juvenile mortality or dispersal has occurred, and individuals must survive

to the next breeding season before reproducing. I distinguished two life-history stages:
(a) juveniles (fledged individuals < 1-yr-old) and (b) adults (≥ 1-yr. old), based on the age
of sexual maturity (Fig. 2).
I used Lefkovitch stage-structured matrices, which provide a way of
incorporating demographic data into a structured population growth model (Caswell
2001). Specifically, demographic rates determine a transition matrix A containing the
annual probabilities of surviving and reproducing. Multiplying the matrix by a vector nt,
containing the distribution of individuals among stage classes in year t, gives the resulting
distribution the following year. Successive years can then be modeled by iterating this
equation. In their unmodified form, Lefkovitch matrices assume a single, well-mixed
population lacking spatial structure and density-dependence, assumptions that I relax
below. They also assume homogeneous probabilities of survival and reproductive
success within each stage (Caswell 2001). All simulations were carried out using
MATLAB (Mathworks 2008).
ff

Puget
Lowlands

s2(b)

Coastal

ff

mba

Juveniles

ff

Adults

s1(b)

s1(a)
mca
s2(a)

Adults

Juveniles

ff

mab
mbc
mcb

mac
Adults

Juveniles

s1(c)

ff

ff

s2(c)

Columbia
River
Figure 2. Simplified life cycle diagram of E. a. strigata used to guide the structure of the models
developed in this paper. Fledglings survive to become adults according to a local juvenile
survival probability,s1(i). Birds breed at one year of age according to a fecundity transition
probability, ff, and survive thereafter according to a local adult survival probability, s 2(i). The
probability of juvenile dispersal is indicated by mij and represented by the dotted lines.
Parameter abbreviations used in the models are presented in italics.

As was you will see below, the location of active nesting areas and the observed
variation in local demography (see Parameters below) suggests the presence of at least
three subpopulations, a violation of the assumption of homogeneity in unmodified
Lefkovitch matrix models. Specifically, both juvenile and adult survival depends on
which subpopulation an individual occupies. Thus, I developed a spatially-structured
matrix accounting for differences and spatial correlations in survivorship and movement
between subpopulations. The matrix had the following form:
𝒇𝒇𝒔𝟏 (𝟏 − 𝒎𝒃𝒂 − 𝒎𝒄𝒂) 𝒇𝒇𝒔𝟐 (𝟏 − 𝒎𝒃𝒂 − 𝒎𝒄𝒂)

𝐹1 𝑚𝑎𝑏

𝐹2 𝑚𝑎𝑏

0

0

𝒔𝟏(𝒂)

𝒔𝟐(𝒂)

𝐹1 𝑚𝑏𝑎

𝐹2 𝑚𝑏𝑎

0

0

𝒔𝟏(𝒃)

𝒔𝟐(𝒃)

𝐹1 𝑚𝑐𝑎

𝐹2 𝑚𝑐𝑎

𝐹1 𝑚𝑐𝑏

𝐹2 𝑚𝑐𝑏

0

0

0

0

𝐹1 𝑚𝑎𝑐

𝑭𝟏 (𝟏 − 𝒎𝒂𝒃 − 𝒎𝒄𝒃 ) 𝑭𝟐 (𝟏 − 𝒎𝒂𝒃 − 𝒎𝒄𝒃 )

𝐹2 𝑚𝑎𝑐

0

0

𝐹1 𝑚𝑏𝑐

𝐹2 𝑚𝑏𝑐

0

0

𝑭𝟏 (𝟏 − 𝒎𝒂𝒄 − 𝒎𝒃𝒄 ) 𝑭𝟐 (𝟏 − 𝒎𝒂𝒄 − 𝒎𝒃𝒄 )
𝒔𝟏(𝒄)

𝒔𝟐(𝒄)

where si is the survival rate of individuals in class i; ff is the number of female fledglings
per reproductive female; and m is the probability of juvenile dispersal, mij, from each site
j to site i. The product {ffsi} is the reproduction element (Fi) of each stage-class. Note
that the reproduction element for juveniles arises from the fact that birds in their first year
of life begin the year as fledglings, but by the end of that year have produced fledglings
of their own (Figure 1). The 2 x 2 sub-matrices (in bold) on the diagonal from the upper
left to lower right corners represent the Puget Lowlands, Coastal, and Columbia River
subpopulations (denoted by subscripts a, b, and c, respectively).
I incorporated density-dependence in each of the two models using a population
ceiling that limited the number of females. This approach was feasible for E. a. strigata
due to its territorial behavior, which may limit the number of breeders as determined by
the availability of habitat. The ceiling was calculated as two times the initial number of
females in each subpopulation and was representative of the carrying capacity (K) of each
subpopulation. The ceiling was enforced only when the number of females exceeded K.
If so, the model proportionally truncated individuals of both stages to the number allowed
before proceeding with matrix-vector multiplication for the next time step.

PRIMARY DATA SOURCES
At a minimum, realistic population modeling requires accurate estimates of vital
rates for the population in question (Doak et al. 1994). To estimate survival rates for E.
a. strigata, I relied heavily on mark-resight data gathered for the Washington Department

;

of Fish and Wildlife under the supervision of Dr. Scott Pearson. During the breeding
seasons of 2002-2006, E. a. strigata were studied at seven sites in southwestern
Washington, yielding data I used to estimate site- and stage-specific survival and
dispersal rates (see Pearson et al. 2008 for description of sites). Estimates of fecundity
were taken directly from Camfield et al. (2007). I relied on unpublished government
reports (e.g. Pearson and Hopey 2004, 2005) for estimates of additional parameters.

POPULATION SPATIAL STRUCTURE
Many of the seven sites studied by Pearson et al. (2008) exist in the form of
clusters, where discrete nesting areas are situated close to each other. Observations of
banded individuals suggest that among-area movements are common and that physically
isolated habitats remain functionally connected. Thus, each study site could not be
regarded as a distinct subpopulation. To delineate subpopulations, I created a 15-km
buffer around each recorded nesting site and considered those sites whose buffers
overlapped to constitute a single subpopulation. Consequently, the model contained 3
subpopulations (Fig. 3), where the minimum distance separating any two after applying
the buffer was 70-km. I selected a 15-km buffer because E. a. strigata have a relatively
high degree of site fidelity to previous nesting sites (Pearson et al. 2008). If breeders do
change colonies within the same year, they are unlikely to move substantial distances.
Inter-annual movements, which may involve longer distances, were modeled as dispersal.
Puget
Sound

Coastal
Pacific
Ocean

Puget
Lowlands

Washington

15 km

Columbia
River
Oregon

Figure 3. Spatial distribution of the subpopulations (shaded
areas) making up the E. a. strigata population in
Washington. Arrows indicate infrequent migration between
subpopulations, the consequences of juvenile dispersal.

PARAMETER ESTIMATES
Demographic Parameters
I utilized the program MARK (White and Burnham 1999) to estimate annual
rates of apparent survival and movement probabilities for juveniles and adults using the
first-order Markovian multi-strata model (Brownie et al. 1993 and Hestbeck et al. 1991).
This model is an extension of the Cormack-Jolly-Seber model (Cormack 1964, Jolly
1965, Seber 1965) that permits inclusion of categorical data (e.g. discrete subpopulations)
in the encounter histories that can change during the life of an individual (Lebreton and
Pradel 2002). The primary advantage of multi-strata models is that they yield estimates
of apparent survival that are specific to each categorical state, while also estimating the
probabilities of changing states (Sandercock 2006). I assumed that apparent survival in
the interval time t to t + 1 did not depend on the site at time t – 1 (Brownie et al. 1993). I
also assumed that individuals made site transitions near the end of each time interval and
that individuals did not immigrate to sites outside of the study area (Sandercock 2006).
Apparent survival rates and movement probabilities were calculated using a
reduced model where the probability of encounter was held constant (Table 1).
Information-theoretic procedures using AICc’s provided by Program MARK were
utilized for model selection (Burnham and Anderson 2002). I used variance
discounting to separate sampling variance from the process variance originating from
demographic and environmental stochasticity (Burnham et al. 1987). In Program
MARK, this procedure is fairly straightforward using random-effects models, which
estimate and minimize sampling variation (White et al. 2001). Ideally, the result is an
uncontaminated estimate of the environmental variance around each mean survival rate.
Dispersal in both models referred to the proportion of juveniles that transitioned
from one subpopulation to another by the start of the following year (Bailey and Mock).
This approach was especially feasible because the resighting probability in the markresight study was high (Pearson et al. 2008). For site transitions that had no observed
dispersal (b→a; c→a, b→c; c→b; a→c), I arbitrarily assigned an annual movement
probability of 0.0001 (Table 1). Fecundity rates were taken directly from Camfield et al.
(2008), who estimated the number of female fledglings per female per year. Because this
analysis did not provide site-specific estimates of fecundity, I assigned a uniform value
across subpopulations (Table 1). Thus, the reproductive rate in each subpopulation was
similar, with differences due only to the values of stage specific survival rates (see matrix
formulation).

Table 1. Vital-rates used to parameterize a stage- and space-structured model for E. a. strigata.
Parameter
Stochastic Form
Mean
s1(a) = annual survivorship of juveniles, Puget Lowlands*
beta random variable (rv) 0.137
s1(b) = annual survivorship of juveniles, Coast*
beta rv
0.210
s1(c) = annual survivorship of juveniles, Columbia River*
beta rv
0.100
s2(a) = annual survivorship of adults, Puget Lowlands*
beta rv
0.444
s2(b) = annual survivorship of adults, Coast*
beta rv
0.714
s2(c) = annual survivorship of adults, Columbia River*
beta rv
0.272
ff = number of fledglings produced per female, per year†
constant
0.910
mab = probability of juvenile moving from site b to a*
constant
0.0001
mac = probability of juvenile moving from site c to a*
constant
0.0001
mba = probability of juvenile moving from site a to b*
constant
0.1100
mbc = probability of juvenile moving from site c to b*
constant
0.0001
mca = probability of juvenile moving from site a to c*
constant
0.0001
mcb = probability of juvenile moving from site b to c*
constant
0.0001

Variance
0.110
0.130
0.071
0.163
0.116
0.179

* Source: Estimated from field data provided by S.F. Pearson.
† Source: Camfield et al. 2007.

Carrying Capacity, Initial Abundances, and Stage Distribution
As described earlier, the carrying capacity of each subpopulation was defined as
the maximum number of territories for breeding females, and was assumed to be two
times the initial number of females in each subpopulation (Table 2). This assumption
was largely arbitrary. There is little consensus on the average territory size of a nesting
female, mostly because of the tendency of territory sizes to vary with habitat type and
quality (Stinson 2005). Thus, I chose to set the carrying capacities at such levels where
they would not qualitatively alter the conclusions of each model. However, I did test for
the significance of this assumption, as described below in the sensitivity analysis. I
estimated initial abundances of females in each subpopulation by dividing the numbers
reported in Pearson and Hopey (2004, 2005) by two, assuming a one-to-one sex ratio.
Because both the deterministic and stochastic models eventually converge to a stablestage distribution, I distributed this total number to stage-classes by solving analytically
for the dominant right eigenvector of each sub-matrix (Caswell 2001).
Table 2. Subpopulation-specific parameters used in simulations.
Initial Abundance
a
Subpopulation
Total
Juveniles Adults
Carrying Capacity
Puget Lowlands
111
53
58
222
Washington Coast
43
20
23
86
Columbia River
34
16
18
68
a
Sources: Pearson and Hopey 2004, 2005

DETERMINISTIC MODEL
The mean parameter values shown in Table 1 were used to construct the
deterministic model. Because the dynamics of a population in a constant environment are
captured by the iteration of a single transition matrix (Caswell 2001), modeling was a
basic exercise in matrix multiplication. The deterministic growth rate λd (annual rate of
population change) for each subpopulation was calculated numerically by projecting the
number of individuals over 25 years (Caswell 2001). I quantified the dependence of the
deterministic model’s behavior on each parameter value by performing a two-step
prospective perturbation analysis (Horvitz et al. 1997). First, I calculated sensitivity
coefficients of λd to 5% changes (Mills et al. 1999) in: (1) the reproductive rate of
females; (2) the survival rates of juveniles and adults, in each subpopulation; and (3)
dispersal rates. This form of sensitivity analysis differs from analytical sensitivities
based on eigenvalues (Caswell 2001) in that the parameters analyzed are not limited to
the matrix elements and thus may include underlying vital rates and other parameters. I
altered one parameter at a time, resetting each to its mean value before changing to the
next and running the model again. The sensitivity coefficient of λd to each parameter was
estimated as:

𝑆𝑥𝑑 =

λ 𝑑 ,𝑝𝑒𝑟𝑡𝑢𝑟𝑏𝑒𝑑 −𝜆 𝑑 ,𝑚𝑒𝑎𝑛
𝑥 𝑝𝑒𝑟𝑡𝑢𝑟𝑏𝑒𝑑 −𝑥 𝑚𝑒𝑎𝑛

=

𝛥𝜆 𝑑
Δ𝑥

;

where x is the parameter being perturbed. After calculating the 𝑆𝑥𝑑 values for each

parameter, deterministic elasticities were found using:

𝐸𝑥𝑑 =

𝑥
𝜆𝑑

𝑆𝑥𝑑 .

Next, because the elasticity and sensitivity results are dependent on the vital-rate
estimates used to make the matrix, inaccuracy in these estimates may have a substantial
influence on the sensitivity results. The simplest way to explore the importance of
uncertainty was to modify the method just described by simulating a set of random
matrices where the parameters were drawn from a range of values for each vital-rate,
assuming a uniform distribution. This is essentially a parametric bootstrap method to
place confidence limits on the elasticities by asking how much of the variation in each
growth rate is explained by the uncertainty in each vital-rate (Morris and Doak 2002).
Minimum and maximum survival values were set to the 95% confidence limits calculated
directly from the field data. The range of fecundity values was taken directly from
Camfield et al. (2007).

STOCHASTIC MODEL
One option for incorporating temporal variability is to simulate variation by
randomly drawing from two or more projection matrices, each containing transition
probabilities estimated over a separate time step (Morris and Doak 2002). There are two
problems with this approach. First, data sets such as the one available for E. a. strigata
do not allow the estimation of every entry in the matrix for all time steps or locations.
This can result from poor experimental design or simply small sample sizes. If only the
demographic data corresponding to a complete projection matrix is used, a great deal of
information is lost (i.e. all data from sites or years that did not yield complete projection
matrices). Second, the used of fixed matrices, each corresponding to an actual year of
field work, enforces a very specific set of correlations in the variation of vital rates.
Although these correlations exactly reflect the observed data, the use of fixed matrices
prevents an exploration of the importance of correlation structure itself, apart from the
variances represented among matrices (Doak et al. 1994).
To more flexibly simulate environmental stochasticity in survival rates, variation
in the si(t)’s was represented as realizations of the beta distribution. The beta distribution
is a family of continuous probability distributions confined to the interval [0, 1] and is
appropriate for simulating survival rates, which are probabilities of binary events. I
simulated random beta values with the means and variances specified in Table 1 using a
method outlined in Morris and Doak (2002). This approach relies on the formula for the
cumulative distribution function of the beta:
𝑝 = 𝐹 𝑥 𝑎, 𝑏 =

𝑥 𝑎 −1
1
𝑡 (1
𝛽 (𝑎,𝑏) 0

− 𝑡)𝑏−1 𝑑𝑡;

where 𝛽 𝑎, 𝑏 is the beta function built-in to MATLAB and a and b are transformations
of the mean and variance, respectively, given by:
𝑎 = 𝑠𝑖

𝑠𝑖 (1−𝑠𝑖 )
𝑣𝑎𝑟 (𝑠𝑖 )

− 1 ; and 𝑏 = (1 − 𝑠𝑖 )

𝑠𝑖 (1−𝑠𝑖 )
𝑣𝑎𝑟 (𝑠𝑖 )

−1 .

For a beta distribution with parameters a and b, and a value p between 0 and 1, this
function gives the probability, 𝐹 𝑝 𝑎, 𝑏 , that a random value chosen from the
distribution will be less than or equal to p. In five steps: (1) I began with a uniform
random number between 0 and 1 of F(p│a ,b); (2) I chose a second random number x1
between 0 and 1 and calculated F( x1│a ,b); (3) if F( x1│a ,b) was greater than F(p│a ,b),
I chose a random number x2 between 0 and x1; if F( x1│a ,b) was less than F(p│a ,b), I
chose a random number x2 between x1 and 1; (4) I then set x1 = x2 and recalculated

F( x1│a ,b); (5) I repeated steps 3 and 4 until F( x1│a ,b) - F(p│a ,b) was less than 0.01
which provided a good approximation of p. Rather than searching for a new value for
each beta-distributed vital rate in each year, I began by storing a set of 99 beta values for
each survival rate at the beginning of a simulation. These values were spaced uniformly
by their cumulative distribution function values [F(p│a ,b) = 0.01, 0.02,…,0.99]. I then
drew random numbers for the F(p│a ,b) value of each rate in each year and chose the
stored rate with the closest F value. Morris and Doak (2002) found with a large number
of stored values, this simplification does not influence the outcome of a stochastic
simulation in any appreciable way and also provides the correct association to simulate
correlated beta-distributed values.
To generate values of correlated beta random variables, I followed an approach
developed in Gross et al. (1998) and expanded by Morris and Doak (2002). Because I
lacked time specific estimates of fledgling production, I included only correlations
between juvenile and adult survival within and across subpopulations, while treating
fecundity as an independent parameter. Program MARK provided Pearson correlation
coefficients for every pairwise combination of juvenile and adult survival, which were
used to create a correlation matrix (C). This is a two-step process that uses correlated
normal random variables to generate correlated survival rates with the beta distribution. I
began with the correlation matrix estimated from the set of survival data. Taking the
eigenvalues and right eigenvectors, I decomposed C into two matrices: W was a matrix
with columns as the right eigenvectors of C, and D was a matrix with diagonal elements
equal to the eigenvalues of C. Then, C1/2 = W * D1/2 *W’ where D1/2 was a matrix with
the square-roots of each eigenvalue on the diagonal and W’ was the transpose of matrix
W. Using C1/2, it was possible to take a set of uncorrelated standard normal values and
convert them into a set of correlated normal values. With m as a column vector of
uncorrelated random values from a standard normal distribution, y = C1/2 * m was a
column vector containing now correlated standard normal variables, according to C.
Next, by generating multiple m vectors with different random values and then
applying y = C1/2 * m to each, I simulated a new vector of correlated values, yt for each
year of a stochastic simulation. Then, I used these correlated standard normal values in
each year yt to make a vector of corresponding values from the cumulative distribution
function Fx for each standard normal. Abramowitz and Stegun (1964) provide an
approximation for Fx that essentially computes the relative position of each random value
along the range of possible values that a normal variable can exhibit. Finally, I used

these Fx values to find corresponding values for each of the survival rates by searching
through the 99 beta values stored at the beginning of each simulation. This provided a
good approximation of each pairwise correlation while preserving the cumulative
distribution function of each beta random variable. This approach is solely
phenomenological in that it describes the pattern of correlation evident in the data, but
does not address any specific mechanism underlying that correlation (Doak et al. 1994).
I used stochastic simulations to estimate sensitivity and elasticity coefficients
using the same method described for the deterministic model. In addition to the stageand subpopulation-specific vital rates, I calculated sensitivity coefficients of λs to 5%
changes in: (1) the variance of each simulated survival rate; and (2) the carrying capacity
of each subpopulation. For stochastic simulations, I began with the subpopulation sizes
shown in Table 2. Each simulation consisted of 10,000 replications for 25 years. I
avoided projecting further into the future for two reasons: (1) the very low growth rates
seen in the deterministic model made projections beyond this time frame irrelevant; and
(2) given the results of Fieberg and Ellner (2000), the confidence intervals around the
probability of extinction can easily encompass the entire range from zero to one as the
time horizon approaches a century. Therefore, seeking to ensure the persistence of E. a.
strigata by predicting and managing over an initial short time interval seemed more
appropriate given the limitations of modeling.

STOCHASTIC SIMULATIONS OF MANAGEMENT SCENARIOS
Using the same methods of the prospective perturbation analysis, I explored what
levels of adult and juvenile survival and fecundity are required to yield estimates of λs ≥ 1
in each subpopulation. This procedure involved running a series of simulations with
perturbed estimates of each of the three vital-rates (s1, s2, ff) based on biologically
realistic estimates producing a viable population of E. a. articola, a closely related
subspecies of Horned Lark found in British Columbia (Camfield et al. 2007). First, I ran
two simulations for each vital-rate increased by 10 and 20%. Next, I ran a set of
simulations where I incrementally perturbed two vital rates in concert by 5 and 15%. A
total of twelve simulations represented each possible combination of two vital rate
manipulations (s1 - s2; s1 - ff; s2 – ff) increased to both 5 and 15%. I avoided running
simulations where the variance or spatial correlations between vital rates was perturbed.
While both have serious implications for the persistence of small populations, it is rare
that either can be altered through management (Morris and Doak 2002).

RESULTS
DETERMINISTIC MODEL
Projections of the deterministic model suggest that the three subpopulations will
not sustain themselves for much longer than two decades (Fig. 4). Each declined to
extinction by the twenty-fifth year regardless of the initial subpopulation size or stagedistribution. While the statewide population was projected to decline 39% per year
(λd = 0.61), the rate of decline of each subpopulation varied considerably (Table 4).
Because the annual fecundity rate was set equal in each sub-matrix, disparities among
local growth rates were attributable entirely to differences in survival and movement
probabilities. The most striking result was how rapidly the Columbia River and Puget
Lowlands subpopulations declined to extinction: neither persisted more than 10 years.
Thus, for most of the projection, an extant statewide population was made possible only
by the persistence of a small number of individuals (< 20) in the Coastal subpopulation.
I determined the sensitivity and elasticity of the deterministic model to small
changes in vital rates as described previously. Sensitivity is defined as the absolute
change in population growth rate per unit change in a parameter while all others are held
constant (Wootton and Bell 1992). Elasticity measures the proportional response of
population growth given a proportional change in a parameter, again, all other rates are
held constant (Doak et al. 1994). Note that I calculated the sensitivities and elasticities of
the measurable vital rates (Table 1), not the combinations of these rates that comprise the
actual elements of the projection matrix.
Puget Lowlands

Coastal

Columbia River

Statewide

200

No. of Juveniles and Adults

180
160
140
120
100
80
60
40
20
0
1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

No. Years

Figure 4. Projection in time of each of the subpopulations (deterministic
model), with the statewide decline shown in bold.

The prospective perturbation analysis indicates that population growth is far
more sensitive to changes in survivorship than to fecundity (Table 3). The maximum
sensitivity always corresponded to the survival rates of adults in each subpopulation. The
model was moderately sensitive to juvenile survival, and less sensitive to fecundity and
juvenile dispersal. Also, the relatively small response of each local deterministic growth
rate to changes in fecundity and the large sensitivity to changes in survival rates
emphasizes the need to better quantify both juvenile and adult survival, given the low
sample sizes currently available to calculate these rates.
Although elasticities aptly summarize the relative effects of small modifications
in different demographic rates, they do not directly reflect how uncertainty in each rate
may alter the reliability of the perturbation analysis. As described previously, I assessed
the effects of this uncertainty by simulating selected vital rates from uniform distributions
bounded by the 95% confidence interval around each mean estimate. The maximum
local λd’s, given the uncertainty in survival and fecundity rates are shown in Table 4. The
variation that did occur in each deterministic growth rate was largely generated by
uncertainty in adult survival values, indicating that error in this rate was the most
important in influencing the results of the sensitivity analysis.

Table 3. Local deterministic growth rates (λd) and sensitivity and elasticity coefficients of vital
rates (defined in Table 1). 95% confidence limits for elasticity values are shown in
parentheses. s1(i) = annual survivorship of juveniles; s2(i) = annual survivorship of adults; ff =
annual fecundity; mba = probability of juvenile moving from site a to b.
Parameter
Sensitivity
Elasticity (95% C.I.)
Puget Lowlands (λd =0.55)
s1(a)
s2(a)
ff
mba
Coastal (λd =0.89)
s1(b)
s2(b)
ff
mba
Columbia River (λd =0.41)
s1(c)
s2(c)
ff

0.809
1.000
0.122
-0.122

0.201
0.799
0.201
-0.025

(0.075, 0.437)
(0.562, 0.924)
(0.075, 0.437)
(-0.153, -0.008)

0.910
1.000
0.214
0.091

0.2143
0.7857
0.2143
0.031

(0.110, 0.493)
(0.508, 0.889)
(0.110, 0.492)
(0.005, 0.168)

0.910
1.000
0.111

0.316 (0.053, 0.911)
0.684 (0.088, 0.946)
0.316 (0.053, 0.911)

Note: movement probabilities whose values were arbitrarily assigned are not included. Elasticity
coefficients do not sum to 1 because, as described above, I calculated the sensitivities and elasticities of
the underlying vital rates, as opposed to the elements of the projection matrix.

Table 4. Results of simulations to determine the sensitivity of local population growth to
increases in vital rates. The maximum local λd, given the uncertainty in each parameter, is
shown in boldface.
Vital Rate
Puget Lowlands
s1(a)
s2(a)
ff
Coastal
s1(b)
s2(b)
ff
Columbia River
s1(c)
s2(c)

Maximum vital
rate value

Maximum potential value of
local λd

Percentage variation in
local λd explained

0.265
0.740
1.080

0.681
0.865
0.588

6
90
1

0.455
0.927
1.080

1.128
1.122
0.946

21
76
3

0.327
1.000

0.516
1.091

9
90

1.080

0.440

0.1

ff

More importantly, in simulations of the Coastal and Columbia River
subpopulations, uncertainty in either juvenile or adult survival yielded growth rates in
excess of one, thus encompassing increasing population dynamics. Consequently, I
decided to rerun the deterministic model using estimates of survival pooled over all sites
(projection shown in Fig. 5, survival values taken directly from Camfield et al. 2007).
The time to extinction projected by the two models differed considerably, suggesting that,
particularly for the Coastal subpopulation, the initial model is probably optimistic.
Though it fails to account for spatial differences in demography, due to the imprecision
generated by small sample sizes, the second projection may provide a more reliable
summary of the probable statewide population decline.
Puget Lowlands

Coastal

Columbia River

Statewide

No. of Juveniels and Adults

250

200

150

100

50

0
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

No. Years

Figure 5. Projection in time of each of the subpopulations using pooled estimates of
survival parameters. Juvenile survival in each subpopulation was set to 0.16; adult
survival was set to 0.47 (taken directly from Camfield et al. 2007).

STOCHASTIC MODEL
The stochastic model predicted population declines similar to those suggested by
the deterministic model. While the average statewide and local growth rates were
relatively unaffected by the addition of environmental stochasticity (compare Tables 4
and 5), it is clear that the addition of temporal variation tends to broaden the distribution
of population sizes (Fig. 6). Although this distribution ultimately shrank as a result of the
convergence on zero produced by extinctions of the Puget Lowlands and Columbia River
subpopulations, over the first decade, much uncertainty was evident.
Table 5. Some results of the simulations (values from 10,000 replicates).
Stochastic
Time of first
Median time to
growth rate extinction (yrs.) extinction (yrs.)*
Statewide
0.58
13
17 (13-22)
Puget Lowlands
0.51
7
9 (7-12)
Coastal
0.83
13
17 (13-22)
Columbia River
0.39
3
4 (3-7)
* Range shown in parentheses

Extinction
prob. at 25 yrs.
100
100
100
100

Here, I show the full distribution of population trajectories (Fig.’s 6 & 7), rather
than simply expressing the mean ending population size, because my concern was with
the possibility of extinction or near extinction in any one realization of the future (i.e.,
any one simulated replication). Stochastic projections for each subpopulation are shown
in Fig. 7. The Columbia River subpopulation was the first to decline to extinction,
persisting on average less than five years. The Puget Lowlands also declined to
extinction rapidly: by year 10 approximately 90% of replicates had resulted in extinction.
The consequence of both early extinctions was startling: by the fifth year the statewide
population had declined from 191 individuals to, on average, 35 individuals. As was the
case with the deterministic model, the Coastal subpopulation was the most persistent
(Table 5) and accounted for all individuals statewide beyond year 12.
200

Figure 6. Stochastic
projection in time of the
statewide population,
indicated by the solid
line (mean values).
The dashed lines
represent the upper and
lower bounds of the
95% confidence
intervals around the
mean estimate.

No. of Juveniles and Adults

180
160
140
120
100
80

60
40
20
0
0 1

3

5

7

9

11

13

No. Years

15

17

19

21

23

25

Figure 7. Stochastic projection in time of each of the subpopulations, indicated by
the solid line. Again, the dashed lines represent the lower and upper bounds of the
95% confidence intervals around each mean estimate.

No. of Juveniles and Adults

50

100
80
60
40
20

45
40
35
30
25
20

15
10
5

0

0
0

1

2

3

4

PUGET LOWLANDS

5

6

7

8

9 10 11 12 13 14 15

0 1 2 3 4 5 6 7 8 9 10111213141516171819202122232425

COASTAL

No. Years

No. Years

40

No. of Juveniles and Adults

No. of Juveniles and Adults

120

35
30
25
20
15
10
5
0

0

1

2

COLUMBIA RIVER

3

4

5

6

7

8

9

10

No. Years

Sensitivity analysis of the stochastic model was qualitatively similar to the
deterministic model, in that a given proportional change in adult survival caused a much
larger effect on population growth than did changes in other parameters (Table 6). The
model was moderately sensitive to fecundity, juvenile dispersal, and the variances of each
survival rate. More striking, however, was the complete lack of sensitivity of each
subpopulation to the local carrying capacity, strongly suggesting that each is
demographically, rather than habitat, limited. In other words, the number of females
needed to achieve population stability is not limited by the availability of nesting habitat.
As I did for the deterministic model, I reran the stochastic model using estimates
of survival pooled over all sites. The statewide projection is shown in Fig. 8. Again,
though failing to account for spatial differences in demography, this simulation suggests
that the initial model is likely optimistic. By switching to the pooled estimates of
survival, by the tenth year more than three-quarters of the 10,000 replicates had declined
from 191 individuals to less than 5.

Table 6. Sensitivity and elasticity coefficients of local stochastic growth rates (λs) to
model parameters and variances.
Sensitivity

Elasticity

0.809
0.961
0.130
-0.254
-0.330
-0.160
0.000

0.201
0.740
0.208
-0.051
-0.035
-0.011
0.000

0.841
1.020
0.195
0.236
-0.230
-0.172
0.000

0.190
0.812
0.201
0.029
-0.015
-0.012
0.000

0.881
1.011
0.103
-0.162
-0.353
0.000

0.204
0.797
0.217
-0.021
-0.048
0.000

200

1

180

0.9

Cumulative prob. of extinction

No. of Juveniles and Adults

Parameter
Puget Lowlands
s1(a)
s2(a)
ff
mba
variance of s1(a)
variance of s2(a)
K
Coastal
s1(b)
s2(b)
ff
mba
variance of s1(b)
variance of s2(b)
K
Columbia River
s1(c)
s2(c)
ff
variance of s1(c)
variance of s2(c)
K

160
140
120
100
80
60
40
20
0

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

0

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15

No. Years

Figure 8. Stochastic projection of the statewide
population, indicated by the solid line, using
pooled estimates of survival from Camfield et al.
2007. The dashed lines represent the 95%
confidence intervals around the mean estimate.

0 1

3

5

7

9

11

13

15

17

19

21

23

25

No. Years

Figure 9. Cumulative extinction probability of
the statewide population, according to the
initial model (solid line) and model using
pooled estimates of survival taken from
Camfield et al. 2007 (dashed line).

SIMULATION SCENARIOS
One advantage of modeling is that it allows simulation of hypothetical scenarios
as a way to quantify the consequences of changes in the environment or in management
strategies. As described previously, I considered different scenarios related to the
effectiveness of management activities by simulating a priori positive changes in the
demographic rates of each subpopulation. Local stochastic growth rates for two levels of
independent vital-rate improvements are shown in Figure 10a-c. Clearly, the Coastal
subpopulation is the most responsive to management. Modest increases (10%) of each
demographic rate yielded growth rates of λs ≥ 0.90. Though the subpopulation was still
projected to decline to extinction, the increase in λs was sufficient to lengthen the median
time to extinction by more than 5 years. A more optimistic improvement (20%) in adult
survival was sufficient to yield an average estimate of λ s >1. Thus, it appears that
increasing adult survival by 15-20% might, alone, be enough to halt declines of the
Coastal subpopulation. The Puget Lowlands and Columbia River subpopulations were
less responsive to independent manipulations. Even under the optimistic scenario (20%
increases), the benefits realized were not encouraging. Still, the results of these large
scale perturbations agree well with the sensitivity analysis; adult survival is of primary
importance to E. a. strigata’s viability.
(B.) COASTAL

0.68
0.64
0.6
10% Increase

0.56

20% Increase

0.52

Stochastic Growth Rate

Stochastic Growth Rate

(A.) PUGET LOWLANDS

0.48

1.05
1
0.95
10% Increase

0.9

20% Increase

0.85
0.8

Juvenile
Survival

Adult Fecundity
Survival

Juvenile
Survival

Adult Fecundity
Survival

Stochastic Growth Rate

(C.) COLUMBIA RIVER
0.53
0.51
0.49
0.47
0.45
0.43
0.41
0.39
0.37

10% Increase
20% Increase

Juvenile
Survival

Adult Fecundity
Survival

Figure 10. The stochastic growth rate (λs) produced by independent perturbations of entries in
the projection matrix. The baseline λs for each subpopulation is indicated by the solid line.

Local stochastic growth rates produced by four levels of simultaneous vital-rate
improvements are shown in Figure 11a-c. Again, management activities targeting any
demographic parameter were most effective when concentrated on the Coastal
subpopulation. Any one of four scenarios (Fig. 11b) was sufficient to produce estimates
of λs ≥ 1. In contrast, the difficulty of achieving stable growth (λs = 1) of the other two
subpopulations reiterates the need for dramatic and timely management action.
Comparisons of the different scenarios suggest that the most effective strategy
might be to target adult survival and fecundity in each subpopulation. An increase of
these two vital rates by 15% would raise the statewide stochastic growth rate
approximately 25%, and increase the median time to extinction by almost 10 years.
Generally, smaller yet simultaneous increases yielded greater responses by each
stochastic growth rate than independent perturbations, though both approaches had
difficulty reversing the declines suggested by the baseline simulation.

0.63
5- & 5%

0.58

5- & 15%

0.53

15- & 5%
15- &15%

0.48
s1: s2

s1: ff

Stochastic Growth Rate

(B.) COASTAL

0.68

1.1
1.05
1
0.95
0.9
0.85
0.8
0.75

s2: ff

5- & 5%
5- & 15%
15- & 5%
15- &15%

s1: s2

s1: ff

s2: ff

(C.) COLUMBIA RIVER
Stochastic Growth Rate

Stochastic Growth Rate

(A.) PUGET LOWLANDS

0.6
0.5
0.4

5- & 5%

0.3

5- & 15%

0.2

15- & 5%

0.1

15- &15%

0
s1: s2

s1: ff

s2: ff

Figure 11. The stochastic growth rate produced by simultaneous perturbations of
entries in the projection matrix. Increases shown in the legend indicate, respectively,
how pairs of vital rates were manipulated. The baseline λs for each subpopulation is
indicated by the solid line.

DISCUSSION
Clearly, the subpopulations of E. a. strigata in Washington are in grave danger.
My demographic analyses concur with the work of Pearson et al. (2008), who found that
statewide, the population is declining by approximately 40% per year. Perhaps more
interesting than this troubling conclusion is the way in which environmental stochasticity
affects E. a. strigata’s prognosis. When variation in survival rates at levels estimated
directly from field data was simulated, I obtained broad confidence limits for population
sizes over the next decade. More often than not, environmental stochasticity had a
cumulative tendency to hasten declines and shorten the median time to extinction. Thus,
for E. a. strigata, and perhaps for many other endangered species (Menges 1992; Doak et
al. 1994), models that fail to account for variability may mislead managers into thinking
there is ample time to intervene before a population reaches zero, when in actuality
extinction could be imminent, due simply to the inconsistencies of the environment.
In spite of the increasing number of empirical studies of population structure and
dynamics, the need of adding realism to theoretical models is clear (Gaona et al. 1998).
Intensive field work by Dr. Scott Pearson and others enabled me to avoid some usual
limitations of population models, such as those resulting from the use of generic
computer packages (e.g. VORTEX-Lacy 1993; RAMAS-Akcakaya and Ferson 1992;
ALEX-Possingham and Davies 1995), which out of necessity contain a large number of
assumptions and model behavior rather simplistically (Lindenmayer et al. 1995). For
example, in many models of spatially-structured populations, similar values of
demographic parameters are assigned across subpopulations (e.g., Akcakaya et al. 1995;
Lindenmayer and Possingham 1996), making it difficult to identify disparities among
local growth rates that are attributable to differential survival and fecundity.
Nonetheless, these models are not free of assumptions and limitations; and each
should be considered carefully when evaluating the robustness of the results. First,
assuming that yearly environmental states are independent of one another (i.e. no
environmental autocorrelation) reduces the validity of the stochastic model if fluctuations
are caused by systematic trends in environmental conditions, such as anthropogenic
habitat deterioration or climate change. I also opted against simulating catastrophic
events, which are a form of variation distinguished from environmental stochasticity by
the magnitude of their effects on demography (Beissinger and Westphal 1998). They
typically result in large population declines and greatly increase the chance of extinction
(Mangel and Tier 1994). In the case of E. a. strigata, I saw little value in simulating such

events since the entire population would likely immediately go extinct, and would thus
complicate the analyses of the simulated management scenarios.
Finally, some demographic parameters, notably the survival rates of the Coastal
and Columbia River subpopulations, are derived from small sample sizes, and the
sensitivity analysis indicates these rates could have a strong influence on population
dynamics. The ability to develop accurate estimates of both the means and variance
about vital rates is dependent on the number of years that a study has been conducted
(Beissinger and Westphal 1998). Variance in growth rates does not begin to stabilize, if
at all, until 8 to 20 years of data have been collected (Pimm 1991). Thus, the four annual
transitions captured by the field study almost certainly do not represent the full range of
environmental variation experienced by this subspecies.
The reality that such limitations and scarce data may compromise the validity of
any viability analysis is unquestionable (Morris and Doak 2002). Indeed, a number of
authors have raised important criticisms about applying quantitative population models to
rare species (e.g. Boyce 1992, Beissinger and Westphal 1998). Thus, a sensible way to
proceed is to consider relative risks (Beissinger and Westphal 1998), rather than basing
decisions on absolute measures of viability. Knowing the exact, quantitative value of E.
a. strigata’s growth rate is of lesser importance than considering a simple, more
qualitative assessment of whether the population will tend to grow or decline. In this
way, my confidence is strengthened by the consistent qualitative output across models.
Despite the fact that quantitative predictions often disagreed, rarely were differences in
the model structure or parameter values sufficient to yield estimates of λ > 1.
While the ultimate purpose of this research is to contribute to the conservation of
E. a. strigata, providing managers with definite proposals of action is beyond the scope
of this paper. Nevertheless, some reflections concerning conservation follow naturally
from the results above. Although the current demographic rates for E. a. strigata warn of
continued decline, my results suggest that preventable anthropogenic impacts may play a
large role in that decline. Indeed, management to bolster adult survival and perhaps
fecundity may foster some level of population stability. However, results of the
sensitivity analysis also emphasize that some well-meaning management strategies, in
particular efforts solely targeting fecundity, are unlikely to be either cost-effective or
biologically sound. The simulation scenarios indicate that the effect of management
activities will differ depending on location. For instance, increases of survivorship and
fecundity only resulted in estimates of λ > 1 in the Coastal subpopulation. By identifying

divergent local dynamics, management actions can be focused where they make the
greatest contribution to the viability of the statewide population (Wootton and Bell 1992;
Doak and Mills 1994; Gaona et al. 1998).
My work also provides information needed to prioritize future research efforts.
For instance, the method I used to simulate the effect of management was quite crude; the
results of such a model will be much more useful if conservation measures are evaluated
empirically in terms of their impact on demographic parameters. Further information on
the spatial and temporal variation in juvenile and adult survival are also especially
important, given the low sample sizes currently available to calculate these rates and their
importance in determining population growth. These data would decrease the uncertainty
about vital rates, which contributed most of the imprecision to model results.
Investigating the causes for relatively low fecundity, a comparatively small yet
significant inhibitor of population growth, should also be a research priority, as should
specific methods to improve each of the depressed vital rates.
In summary, I believe that the case of E. a. strigata is a hopeful one.
Fortunately, the data needed to evaluate the viability of E. a. strigata is vastly better than
that available for most endangered vertebrates. The subspecies has received state
protection and may soon receive attention at the federal level. The Washington
Department of Fish and Wildlife is eagerly investigating and implementing strategies that
may improve demographic rates and my models demonstrate that management
interventions may have positive effects for E. a. strigata. However, present legal
protection aimed at averting the extinction of the subspecies is almost certainly
inadequate. Both the deterministic and stochastic models demonstrate that complacency
is ill-advised; decisive action to quickly improve demographic rates is needed, given the
consistent qualitative conclusions across models and the uncertainty inherent in
predictions of future population trends.

LITERATURE CITED
Abramowitz, M. and I. Stegun. 1964. Handbook of Mathematical Functions, with
Formulas, Graphs and Mathematical Tables. National Bureau of Standards,
Washington, D.C. (Reprinted by Dover Publications, New York, New York, USA).
Ackakaya, H.R. and S. Ferson. 1992. RAMAS/space: spatially structured population
models for conservation biology. Version 1.3. Applied Biomathematics, New York,
New York, USA.
Ackakaya, H.R., M.A. McCarthy, and J.L. Pearce. 1995. Linking landscape data with
Population Viability Analysis: management options for the helmeted honeyeater
Lichenostomus melanops cassidix. Biological Conservation 73: 169-176.
American Ornithologists' Union. 1957. Check-list of North American birds. 5th edition.
American Ornithologists' Union, Washington, D.C.
Askins, R.A. 1995. Hostile landscapes and the decline of migratory songbirds. Science
267: 1956-1957.
Altman, B. 1999. Status and conservation of state sensitive grassland bird species in the
Willamette Valley. Unpublished report to: Oregon Department of Fish and Wildlife,
Corvallis, OR.
Altman, B. 2000 Conservation Strategy for landbirds in the lowlands and valleys of
western Oregon and Washington. Version 1.0. Oregon-Washington Partners in
Flight.
Bailey, E.A. and P.J. Mock. 1998. Dispersal capability of the California gnatcatcher: A
landscape analysis of distribution data. Western Birds 29:351-360(10).
Beason, R. C. 1995. Horned Lark (Eremophila alpestris). in The Birds of North America,
No. 195. A. Poole, and F. Gill, editors. The Academy of Natural Sciences,
Philadelphia, and The American Ornithologists' Union, Washington, D.C.
Beauchesne, S., and J. Cooper. 2003. COSEWIC status report on the horned lark strigata
subspecies Eremophila alpestris strigata. Status report prepared for the Committee
on the Status of Endangered Wildlife in Canada. COSEWIC Secretariat c/o Canadian
Wildlife Service, Environment Canada.
Behle, W.H. 1942. Distribution and variation of the Horned Larks (Otocoris alpestris) of
Western North America. University of California Publications in Zoology
46:205-316.
Beissinger, S. R., and M. I. Westphal. 1998. On the use of demographic models of
population viability in endangered species management. Journal of Wildlife
Management 62:821-841 (12).
Beissinger, S. R., J. R. Walters, D. G. Catanzaro, K. G. Smith, J. B. Dunning, S. M. Haig,
B. R. Noon, and B. M. Stith. 2006. Modeling approaches in avian conservation and
the role of field biologists. The Auk 123:1-56.
Bowles, J.H. 1898. Notes on the Streaked Horned Lark. Osprey 3:53-54.

Boyce, M.S. 1992. Population Viability Analysis. Annual Review of Ecology and
Systematics 23:481-506.
Brook, B.W., J.J. O’Grady, A.P. Chapman, M.A. Burgman, H.R. Ackakaya and R.
Frankham. 2000. Predictive accuracy of population viability analysis in conservation
biology. Nature 404: 385-387 (12).
Brownie, C., J. E. Hines, J. D. Nichols, K. H. Pollock, and J. B. Hestbeck. 1993.
Capture-recapture studies for multiple strata including non-Markovian transitions.
Biometrics 49:1173-1187.
Buell, A.C., A.J. Pickart, and J.D. Stuart. 1995. Introduction history and invasion patterns
of Ammophila arenaria on the north coast of California. Conservation Biology
9:1587-1593 (6).
Burnham, K. P., and D. R. Anderson. 1998. Model Selection and Inference: A Practical
Information-Theoretic Approach. Springer Associates, New York, New York, USA.
Burnham, K. P., and D. R. Anderson. 2002. Model selection and multimodel inference:
a practical information-theoretic approach. 2nd Ed. Springer-Verlag, New York, New
York, USA.
Burnham, K. P., D. R. Anderson, G. C. White, C. Brownie, and K. H. Pollock. 1987.
Design and analysis methods for fish survival experiments based on releaserecapture. American Fisheries Society Monograph No. 5. Bethesda, Maryland,
USA.
Camfield, A.F., S.F. Pearson, and K. Martin. 2007. Comparative demography of horned
larks: Conservation implications and consequences for life history studies.
Ecology: Submitted.
Caro, T.M. and M.K. Laurenson. 1994. Ecological and genetic factors in conservation: a
cautionary tale. Science 263: 485-486.
Caswell, H. 2001. Matrix Population Models: Construction, Analysis and Interpretation.
2nd Ed. Sinauer Associates, Sunderland, Massachusetts, USA.
Caughley, G. 1994. Directions in conservation biology. Journal of Animal Ecology 63:
215-244.
Chappell, C.B., M.S. Mohn Gee, B. Stephens, R. Crawford, and S. Farone. 2001.
Distribution and decline of native grasslands and oak woodlands in the Puget lowland
and Willamette Valley ecoregions, Washington. Pages 124-139 in Reichard, S.H.,
P.W. Dunwiddie, J.G. Gamon, A.R. Kruckeberg, and D.L. Salstrom eds.
Conservation of Washington’s native plants and ecosystems. Washington Native
Plant Society, Seattle.
Cormack, R. M. 1964. Estimates of survival from the sighting of marked animals.
Biometrika 51:429-438.
Crawford, R.C., and H. Hall. 1997. Changes in the south Puget prairie landscape. Pages
2-15 in P. Dunn and K. Ewing, (eds.). Ecology and conservation of the south Puget
sound prairie landscape. The Nature Conservancy, Seattle, WA.

Davis, G.J. and R.W. Howe. 1992. Juvenile dispersal, limited breeding sites, and the
dynamics of metapopulations. Theoretical Population Biology 41: 184-207.
Dawson, L.W., and J.H. Bowles. 1909. The birds of Washington. Occidental Press,
Seattle, Washington.
Doak, D.F. and L.S. Mills. 1994. A useful role for theory in conservation. Ecology
75:615-626.
Doak, D.F., P. Kareiva and B. Klepetka. 1994. Modeling population viability for the
desert tortoise in the Western Mojave Desert. Ecological Applications 4:446-460 (3).
Donovan, T.H., and C.H. Flather. 2002. Relationship among North American songbirds
trends, habitat fragmentation, and landscape occupancy. Ecological Applications 12:
364-374.
Drovetski, S. V., S. F. Pearson, and S. Rohwer. 2005. Streaked horned lark Eremophila
alpestris strigata has distinct mitochondrial DNA. Conservation Genetics 6:875-883.
Falconer, D.S. 1989. Introduction to quantitative genetics. Third edition. Longman
Scientific and Technical, co-published with John Wiley, New York, New York,
USA.
Fox, G.A. and B.E. Kendall. 2002. Demographic stochasticity and the variance reduction
effect. Ecology 83: 1928-1934 (7).
Ferson, S. and M.A. Burgman. 1995. Correlations, dependency bounds and extinction
risks. Biological Conservation 73: 101-105(2).
Franklin, I.R. 1980. Evolutionary changes in small populations. Pages 135-139 in M.E.
Soule and B.A. Wilcox (eds.). Conservation biology: an evolutionary-ecological
perspective. Sinauer Associates. Sunderland, Massachusetts, USA.
Gaona, R., F. Ferreras and M. Delibes. 1998. Dynamics and viability of a metapopulation
of the endangered Iberian lynx (Lynx pardinus). Ecological Monographs 68:349-370.
Gilpin, M. 1996. Metapopulations and wildlife conservation: approaches to modeling
spatial structure. Pages 11-27 in D.R. McCullough (ed.). Metapopulations and
wildlife conservation. Island Press, Washington D.C., USA.
Gilpin, M. and M.E. Soule. 1986. Minimum viable populations: processes of species
extinction. Pages 19-34 in M.E. Soule (ed.). Conservation biology: the science of
scarcity and diversity. Sinauer Associates, Sunderland Massachusetts, USA.
Ginzburg, L.R., S. Ferson and H.R. Ackakaya. 1990. Reconstructibility of density
dependence and the conservative assessment of extinction risks. Conservation
Biology 4:63-70.
Gross, K.G., J.R. Lockwood, C. Frost and W.F. Morris. 1998. Modeling controlled
burning and trampling reduction for conservation of Hudsonia montana.
Conservation Biology 12: 1291-1301.

Hanksi, I. 1997. Metapopulation dynamics: From concepts and observations to predictive
models. Pages 69-91 in I. Hanski and M.E. Gilpin (eds.). Metapopulation biology:
ecology, genetics and evolution. Academic Press. San Diego, California, USA.
`Henshaw, H.W. 1884. The shore larks of the United States and adjacent territory. Auk
1:254-268.
Hestbeck, J. B., J. D. Nichols, and R. A. Malecki. 1991. Estimates of movement and site
fidelity using mark-resight data of wintering Canada geese. Ecology 72:523-533.
Jewett, S.G., W. P. Taylor, W.T. Shaw, and J.W. Aldrich. 1953. Birds of Washington
State. University of Washington Press, Seattle, Washington.
Johannessen, C.L., Davenport, W.A., Millet, A., McWilliams, S. 1971. The vegetation of
the Willamette Valley. Annals of the Association of American Geographers 61:286302.
Jolly, G. M. 1965. Explicit estimates from capture-recapture data with both death and
immigration in a stochastic model. Biometrika 52:225-247.
Lacy, R.C. 1993. VORTEX: a computer simulation for population viability analysis.
Wildlife Research 20:45-65.
Lande, R. 2002. Incorporating stochasticity in population viability analysis. Pages 18-40
in S.R.Beissinger and D.R. McCullough (eds.), Population viability analysis. Univ.
of Chicago Press.
Lebreton, J.D. and R. Pradel. 2002. Multistate recapture models: modelling incomplete
individual histories. Journal of Applied Statistics 29: 353-369.
Lindenmayer, D.B., H.R. Ackakaya, R.C. Lacy and H.P. Possingham. 1995, A review of
three models for metapopulation viability analysis-ALEX, RAMAS/space, and
VORTEX. Ecological Modelling 82:161-174.
Lindenmayer, D.B. and H.P. Possingham. 1996. Ranking conservation and timber
management options for Leadbeater’s possum in southeastern Australia using
population viability analysis. Conservation Biology 10:235-251.
Ludgwig, D. and C.J. Walters. 1985. Are age-structured models appropriate for catcheffort data? Canadian Journal of Fisheries and Aquatic Sciences 42: 1066-1072.
MacLaren, P.A. and E.B. Cummins. 2000. Streaked Horned Lark Surveys in Western
Washington. Report funded by USFWS.
Mangel, M. and C. Tier. 1994. Four facts every conservation biologist should know about
persistence. Ecology75: 607-614 (2).
McDonald, D.B. and H. Caswell. 1993. Matrix methods for avian demography. Pages
139-185 in D. Power, editor. Current Ornithology. Plenum Press, New York, New
York, USA.
Menges, E.S. 1992. Stochastic modeling of extinction in plant populations. Pages 253276 in P.L. Fiedler and S.K. Jain (eds.), Conservation Biology: the theory and

practice of nature conservation, preservation and management. Chapman and Hall,
New York, New York, USA.
Miller, R.S. and D.B. Botkin. 1974. Endangered species models and predictions.
American Scientist 62: 172-181.
Mills, L.S., D.F. Doak and M. Wisdom. 1999. The reliability of conservation actions
based upon elasticities of matrix models. Conservation Biology 13:815-829(9).
Morris, W. F., and D. F. Doak. 2002. Quantitative Conservation Biology: Theory and
practice of population viability analysis. Sinauer Associates, Inc., Sunderland,
Massachusetts, USA.
Pearson, S.F. 2003. Breeding phenology, nesting success, habitat selection, and census
methods for the streaked horned lark in the Puget lowlands of Washington. Natural
Areas Program Report 2003-2. Washington Department of Natural Resources.
Olympia, WA.
Pearson, S.F. and B. Altman. 2005. Range-wide streaked horned lark (Eremophila
alpestris strigata) assessment and preliminary conservation strategy. Washington
Department of Fish and Wildlife, Olympia, WA.
Pearson, S.F. and M. Hopey. 2004. Streaked horned lark inventory, nesting success and
habitat selection in the Puget lowlands of Washington. Washington State Department
of Natural Resources. Natural Areas Report. 2004-01.
Pearson, S.F. and M. Hopey. 2005. Streaked horned lark nest success, habitat selection
and habitat enhancement experiments for the Puget lowlands, coastal Washington
and Columbia River islands. Washington State Department of Natural Resources.
Natural Areas Report. 2005-01.
Pearson, S. F., A. F. Camfield, and K. Martin. 2008. Streaked horned lark (Eremophila
alpestris strigata) fecundity, survival, population growth and site fidelity: Research
progress report. Washington Department of Fish and Wildlife, Wildlife Science
Division.
Pearson, S. F., M. Hopey, W. D. Robinson, and R. Moore. 2005. Range, abundance and
movement patterns of wintering streaked horned larks in Oregon and Washington.
Natural Areas Program Report 2005-2. Washington Department of Natural
Resources.
Pimm, S.L. 1991. The balance of nature? University of Chicago Press, Chicago, Illinois,
USA.
Possingham, H.P. and I. Davies. 1995. ALEX: a model for viability analysis of spatially
structured populaltions. Biological Conservation 73:143-150.
Ralls, K.J. and B.L. Taylor. 1997. How viable is population viability analysis. Pages 228235 in S.T.A. Pickett, R.S. Ostfeld, M. Shachal and G.E. Likens (eds.). The
Ecological Basis of Conservation. Chapman and Hall, New York, New York, USA.
Ralls, K.J., S.R. Beissinger and J.F. Cochrane. 2002. Guidelines for Using Population
Viability Analysis in Endangered-Species Management. Pages 521-550 in S.R.

Beissinger and D.R. McCullough (eds.). Population Viability Analysis. University of
Chicago Press, Chicago, Illinois, USA.
Ricklefs, R. E. and G. Bloom. 1977. Components of avian breeding productivity. Auk
94:86–96.
Rogers, R.E. 1999. The Streaked Horned Lark in Western Washington. Unpublished
Report. Washington Department of Fish and Wildlife, Olympia, WA.
Rogers, R. E. 2000. The status and microhabitat selection of streaked horned lark,
western bluebird, Oregon vesper sparrow, and western meadow lark in western
Washington. MSc. Thesis. The Evergreen State College, Olympia, Washington.
Sandercock, B.K. 2006. Estimation of demographic parameters from live encounter
data: a summary review. Journal of Wildlife Management 70:1504-1520
Seber, G. A. F. 1965. A note on the multiple recapture census. Biometrika 52:249-259.
Schrott, G. R., K. A. With and A. W. King. 2005. Demographic limitations of the ability
of habitat restoration to rescue declining populations. Conservation Biology 19:
1181-1193.
Shaffer, M.L. 1981. Minimum population sizes for species conservation. Bioscience
31:131-134.
Shaffer, M.L. 1983. Determining minimum viable population sizes for the grizzly bear.
International Conference on Bear Research and Management 5: 133-139.
Shaffer, M.L. and F.B. Samson. 1985. Population size and extinction: a note on
determining critical population size. American Naturalist 125: 144-152.
Sibley, D.A. 2000. National Audubon Society: The Sibley Guide to Birds. Chanticleer
Press. New York.
Stinson, D. W. 2005. Draft Washington State status report for the mazama pocket gopher,
streaked horned lark, and Taylor’s checkerspot. Washington Department of Fish and
Wildlife.
Towle J.C. 1982. Changing geography of the Willamette Valley woodlands. Oregon
Historical Quarterly, 83:66-87.
Tuljapurkar, S.D. 1990. Population Dynamics in Variable Environments. SpringerVerlag. New York, New York, USA.
USFWS (U. S. Fish and Wildlife Service). 2002. Biological and Conference Opinions for
the Columbia River Channel Improvements Project. U.S. Fish and Wildlife Service,
Portland, Oregon.
White, G. C. and K. P. Burnham. 1999. Program MARK: survival estimation from
populations of marked animals. Bird Study 46 Supplement:120-138.
White, G. C., K. P. Burnham, and D. R. Anderson. 2001. Advanced features of Program
Mark. Pages 368-377 in R. Field, R. J. Warren, H. Okarma, and P. R. Sievert,
editors. Wildlife, land, and people: priorities for the 21st century. Proceedings of the

Second International Wildlife Management Congress. The Wildlife Society,
Bethesda, Maryland, USA.

Wright, S. 1969. Evolution and genetics of populations. Volume 2. The theory of gene
frequencies. University of Chicago Press, Chicago, Illinois, USA.
Wootton, J.T. and D.A. Bell. 1992. A metapopulation model for the Peregrine Falcon in
California: viability and management strategies. Ecological Applications 2:307-321.