Appendix F: 

Tinker, M.T., D.F. Doak, J.A. Estes, and A.H. Kage. 2004. A spatially 

explicit simulation model to predict range expansion of the southern sea 

otter. 
This page intentionally left blank. 
 

A SPATIALLY EXPLICIT SIMULATION MODEL TO PREDICT 

RANGE EXPANSION OF THE SOUTHERN SEA OTTER.

Tinker, M.T., D.F. Doak, J.A. Estes, and A.H. Kage. 2004. 
Abstract 

1. Reliable projections of future population growth and southward range expansion of 

the southern sea otter population would be useful for a number of management purposes, 

and are of interest to several state and federal agencies, including MMS. 

2. We developed a spatially explicit simulation model to project population growth and 

southward range expansion. Our model represents a unique synthesis of a multi-state 

dispersal matrix with the integrodifference equation approach to calculating invasion 

speed. The model was parameterized with data from the current study. 

3. We used the results of repeated simulations with this model to predict future patterns 

of range expansion while accounting for uncertainty in model parameters. Although 

model output closely matched historical data on rates of southward range expansion, 

simulation results were highly variable, reflecting the uncertainty of estimates of input 

parameters as well as the uncertainty regarding ultimate causes of variation in survival 

rates. 

4. Males are more likely than females to move long distances, and most of the individuals 

that travel south of Point Conception are males, yet movement rates of females, 

particularly those of juvenile and sub-adult females, had much more impact on both 

population expansion and southward range expansion than those of males. 

5. The survival rate of juvenile/sub-adult females in the southern part of the range 

emerged as a key parameter influencing the rate of range expansion. Fieldwork to 

improve estimates of this parameter would do much to reduce uncertainty in the model’s 

predictions. 

6. Although our model is relatively simple and does not account for many important 

aspects of sea otter biology and ecology, it makes use of all existing demographic and 

movement data and provides a robust and generalizable approach to understanding and 

predicting population dynamics in southern sea otters. This new tool for conservation 

biologist and managers can be easily expanded and improved as additional data and more 

precise parameter estimates become available. 
Introduction 

Data on age- and sex-specific probabilities of survival, reproduction and 

movement provide the basic tools for understanding past and present population 

dynamics (Caswell 2001, Doak and Morris 2002), and can also be used for predicting 

future population dynamics. Such tools are often the basis for conservation and 

management decisions. In the case of the southern sea otter, developing a realistic 

projection of future population growth and range expansion at the southern end of the 

current distribution would facilitate the informed assessment of potential impacts of sea 

otters on important industries (e.g. fisheries, eco-tourism), potential negative effects of 

human activity on sea otters (e.g. risks associated with the near-shore transport and 

extraction of petroleum, entanglement in fishing gear, etc.), and the eventual recovery 

and de-listing of this threatened species (USFWS 2003). One of the most important 

outcomes of the current research efforts has been the accumulation of a detailed and 

extensive database of spatially-explicit demographic and movement information for the 

southern sea otter. Here, we use these data to develop a spatially structured simulation 

model for predicting population growth and southward range expansion. 

 Stage-based projection matrices provide a means of integrating information on 

population structure, individual survival and reproduction in an intuitive and 

mathematically useful way. Projection matrices are commonly used to predict future 

population dynamics, measure the sensitivity of these dynamics to particular vital rates, 

and elucidate the underlying processes responsible for patterns of interest (Caswell 2001). 

Multi-state projection matrices represent an extension of the basic matrix approach, and 

can be utilized when population structure or environmental conditions vary with 

geographic location, or when the effects of individual movements between subpopulations are thought to be important (Lebreton and Gonzalez-Davila 1993). Multistate projection models facilitate the quantitative interpretation of source-sink dynamics 

(e.g. Doak 1995), and can help to clarify the relative importance of survival and dispersal 

in driving population trends (e.g. Lebreton and Gonzalez-Davila 1993). 

Multi-state matrices are often used to study metapopulation dynamics and the 

colonization rate of unoccupied habitat (Caswell 2001). The latter phenomenon can also 
be modeled as a continuous variable using integrodifference equations, as described by 

Neubert and Caswell (2000). This relatively new technique utilizes stage-specific data on 

dispersal and vital rates to derive the asymptotic speed at which the population front (or 

“traveling wave”) will invade empty habitat (Neubert and Caswell 2000). Note that for 

the purpose of our analyses we will define the term “dispersal” in a purely quantitative, 

descriptive sense, referring to the average linear distance moved (or mean net 

displacement) between the location of an individual at time t=1 and the location at t=2. 

This definition makes no reference to the biological cause or behavioral significance of 

the movements, which may often differ between age and sex classes. 

 In the case of the southern sea otter, analyses of multiple data sets indicate that 

demographic rates are not constant across the sea otters range, but vary between broad 

geographic areas (Figure 1; Chapter 2, this report). Moreover, it appears that individual 

dispersal distances also vary by age and sex class (Chapter 3, this report). A multi-state 

projection matrix model would therefore be an appropriate tool for elucidating the 

relative importance of dispersal and survival in driving population dynamics at the range 

peripheries. The interaction between stage-specific dispersal and vital rates will also 

likely determine the rate and pattern of population expansion into unoccupied territory, a 

phenomenon best described using integrodifference equation models (Neubert and 

Caswell 2000). We couple these two techniques in a spatially explicit simulation model, 

parameterized using data from the current study (Chapter 2, this report), and use the 

results of repeated simulations to predict future patterns of range expansion while 

accounting for uncertainty in model parameters. We then use sensitivity analysis to 

determine the relative importance of each model parameter, in order to highlight specific 

areas where further study will be particularly useful. 
Methods 

Matrix Structure 

A stage-based, 2-sex projection matrix (Caswell 2001) was used to describe 

annual transitions between 4 age-classes: juveniles (defined as 1 year post-weaning), subadults (2 and 3-year-olds), prime-age adults (4–10-year-olds) and aged adults (11 years of 

age or older). Specifically, we constructed an 8×8 matrix of the form: 

A = 

 Sex = f Sex = m 

Sex = f 

4,3 4,4

3,2 3,3

2,1 2,2

1,2 1,3 1,4

0 0

0 0

0 0

0

G P

G P

G P

R R R

0 0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

Sex = m

0 0 0 0

0 0 0 0

0 0 0 0

0 R5,2 R5,3 R5,4

8,7 8,8

7,6 7,7

6,5 6,6

0 0

0 0

0 0

0 0 0 0

G P

G P

G P

1 

where each element of A, aj,i, represents the transition from age/sex class i to age/sex 

class j. Note that i ≤ 4 correspond to female age-classes, while i ≥ 5 correspond to male 

age classes. Three types of transition are identified in A: G, P and R. The first type of 

transition, G, represents survival and growth: in other words, the probability that 

individuals survive for 1 year and advance to the next age-class. The second type, P, 

represents “persistence”, or survival without transition to the next age class. The final 

type of transition, R, represents survival and successful reproduction: for our purposes, an 

individual female is considered to have successfully reproduced if she gives birth and 

successfully weans a pup (i.e. she contributes a single viable juvenile to the population). 

To estimate P, G and R we used standard equations for deriving fixed-stageduration transition probabilities from underlying vital rates (Caswell 2001): 

















−

− = ⋅ −

−

( ) 1

( ) ( ) 1

1

, i i

i i i i

s T

s T s T

i i i P s

λ

λ λ 2 

















−

− = ⋅

−

( ) 1

( ) ( ) 1

, i i

i i i i

s T

s T s T

j i i G s

λ

λ λ 3 
j i i i wi R = s ⋅ b ⋅ 2

1 , 4 

where Ti is the stage duration (in years) for age/sex class i, λ is the annual rate of 

population growth, si is the annual survival rate for an individual of stage i, bi represents 

the birth rate for a female of stage i and wi represents the weaning success rate for a 

female of stage i. We assumed a 1:1 sex ratio at birth. 

 Matrix A was used to model the basic demographic processes for sea otters at a 

particular time and place. The next step was to introduce spatial structure to the model 

by constructing a multi-state matrix. Consider a structured population consisting of three 

sub-populations, each of which exhibits unique demographic properties, and between 

which there is no potential for dispersal. A multi-state matrix for such a population 

would be: 

B = 

3,3

2,2

1,1

0 0

0 0

0 0

Α

Α

Α

5 

where each cell of matrix B represents a transition from sub-population x to subpopulation y. The diagonal of B consists of sub-matrices Ay,x (where y = x), each of 

which describes demographic processes for a single sub-population according to equation 

1. All other elements of B are set to 0 because we specified no dispersal between subpopulations. 

Unfortunately, although equation 5 is conceptually very simple, it is also 

unrealistic: for more realistic dynamics we must allow dispersal between subpopulations. Accordingly, let mi,x→y represent the probability of moving to subpopulation y for an individual of stage i that starts the year in sub-population x. These 

probabilities can be incorporated into a new matrix, M, that has a general structure 

identical to A but whose elements correspond to stage-specific dispersal probabilities 

(mi,x→y) rather than survival or reproduction. Dispersal probabilities between subpopulations can be calculated from empirically-derived dispersal distance distributions. 

We assumed that the annual dispersal distance for an individual sea otter of stage i at 

location x′ (i.e. somewhere within sup-population x) could be described by a Laplace 

probability distribution with shape parameter σi,x. The Laplace distribution essentially 
consists of two back-to-back exponential distributions, and was considered appropriate 

for modeling sea otter movements in California because animals are restricted to 

dispersal in one of two directions (north or south along the coast). 

We defined each sub-population spatially by its northern and southern boundaries 

along the ATOS line (the “as the otter swims” line, corresponding to a series of points at 

500m increments, north to south, along the 10m bathymetric contour); thus subpopulation y was spatially defined by boundary points yN and yS, and spanned 0.5(yS – yN) 

km of coastline. The probability that an individual located at point x′ would disperse to 

sub-population y was estimated as the absolute difference between the Laplace 

cumulative distribution function evaluated at values yN–x′ and yS – x′ : 































 − ′ − − 





























 − ′ ′ → = −

i x

S

i x

N

i x y

y x y x

m

, ,

( ) exp

2

1

exp 1

2

1 1

σ σ

6 

and the net probability of dispersal from sub-population x to sub-population y was 

estimated as: 

∑( ) ′=

→ ′ → = ′ S

N

x

x x

i x y i x y m m p(x ) , ( ) 7 

where p(x′ ) sums to 1 for xN ≤ x′ ≤ xS and represents the probability that an individual 

animal from sub-population x would be located at x′. For simplicity we used discreet 

summation, rather than continuous integration, thereby assuming that all points within a 

500m interval (or ATOS unit) would be adequately represented by a single integer value 

of x′ . A uniform spatial distribution of individuals within each sub-population would 

have allowed us to define p(x′ ) = 1/(xS–xN+1); however, examination of annual rangewide survey counts suggested that sea otters were not uniformly distributed but were 

often skewed towards one boundary or clumped within one region. To account for this 

non-uniformity, we used beta probability functions to represent the spatial distribution of 

otters within each sub-population. Specifically, we fit beta distributions to the raw 

survey data for each sub-population (using data from the 5 most recent spring counts; 

Figure 2), having first standardized the location data to values between 0 and 1 (where 0 

corresponded to the northern boundary and 1 corresponded to the southern boundary of 
the sub-population). Using the resulting beta functions we could then calculate the 

probability that an individual from sub-population x would be located at x′ : 

















− +

′ − −















− +

′ − + ′ = 1 , , 1

1 ( ) , ,

N S

N

N S

N

x x

x x

a b

x x

x x

p x β a b β 8 

where β(a,b,Z) represents the beta cumulative distribution function with parameters a and 

b, evaluated at value = Z. 

 Equations 6 – 8 were solved for each pair of sub-populations, including the 

special case when y = x (which corresponds to the probability of remaining within the 

same sub-population). To ensure that ∑ mi,x→y = 1 for each combination of i and x, we 

made the simplifying assumption that individuals whose dispersal distance brought them 

to the range end would “bounce off” this boundary. This was accomplished by a minor 

adjustment to equation 6 when the target sub-population was a range end: when the target 

was the northern-most end of the range, yN–x′ was set to ∞, and when the target was the 

southern-most end of the range, yS–x′ was set to ∞. 

Combining the two types of matrix, M and A, we constructed a multi-state matrix 

which allowed for movements between sub-populations (illustrated again with just three 

sub-populations for simplicity): 

B = 

3,1 3,1 2,3 2,3 3,3 3,3

2,1 2,1 2,2 2,2 2,3 2,3

1,1 1,1 1,2 1,2 1,3 1,3

Α Μ Α Μ Α Μ

Α Μ Α Μ Α Μ

Α Μ Α Μ Α Μ

o o o

o o o

o o o

9 

Each cell of matrix B consists of the Hadamard product of a demographic matrix and a 

movement matrix (Ay,x ° My,x, where ° represents element-by-element multiplication). 

The resulting elements of each sub-matrix therefore represent joint probabilities of 

moving from x to y and then successfully making the transition from stage i to stage j. 

Note that the diagonal of B is mathematically identical to all other cells, though it 

actually represents the special case where individuals do not disperse. For computational 

simplicity we required that individuals disperse at the start of each new year, after which 

survival, growth and reproduction occur at the new location: this results in a demographic 

sub-matrix Ay,x that is identical for all cells in a given row of B. 
Simulating Range Expansion 

The multi-state matrix model in equation 9 accounts for dispersal and 

demographic processes within the existing range of the southern sea otter at time t. It 

does not, however, account for the expansion of the existing range boundaries to the 

north and south. In order to predict the expansion of the population into un-occupied 

territory, we developed an integrodifference equation model (Neubert and Caswell 2000). 

This approach utilizes a stage-based demographic matrix, in conjunction with the 

moment generating functions of stage-specific dispersal kernels (in this case, the Laplace 

distribution functions described above) to solve for the asymptotic speed of the “traveling 

wave” formed by the population front as it moves into empty habitat. The asymptotic 

wave speed has been found to correspond well to the rate of population range expansion 

in both numerical simulations and empirical data sets (Neubert and Caswell 2000); we 

therefore used it to estimate the rate at which the range boundaries of the sea otter 

population will move to the south and north over time. 

To predict southward range expansion, we populated demographic matrix A with 

vital rates corresponding to the southern-most sub-population, and used this in 

conjunction with the appropriate Laplace distribution parameters to solve the 

integrodifference equation for the asymptotic wave speed, following the methods 

outlined by Neubert and Caswell (2000). We used the same approach to predict 

northward range expansion, using vital rates and dispersal kernels corresponding to the 

northern-most sub-population. The resulting estimates of range expansion speed were 

used to re-set the northern-most and southern-most sub-population boundaries on an 

annual basis; this of course had the effect of altering the predicted rates of dispersal to 

and from these sub-populations, and so equations 6 through 9 were re-solved after each 

year of population projection. 

Model Parameterization 

The results of maximum likelihood analyses of mark-recapture data and 10 years 

of carcass age structure data (see Chapter 2) suggest that the southern sea otter population 

consists of 4 sub-populations, identified based on consistent differences in vital rates: 

these correspond to the northern, north-central, south-central and southern portions of the 

current sea otter range (Figure 1). We identify 2 additional areas for the purpose of this 
simulation: these are the population “frontal” areas, defined as recently-occupied areas at 

the northern and southern ends of the range that are currently utilized seasonally and 

exclusively by males (each frontal area spans 30 km of coastline; Figure 1). The resulting 

6 sub-populations were defined spatially as follows: 

1. Northern front, North of Pigeon Pt. (ATOS 90–149) 

2. Northern periphery, Pigeon Pt. – Santa Cruz (ATOS 150–230) 

3. North-central, Santa Cruz – Pt. Sur (ATOS 231–500) 

4. South-central, Pt. Sur – Pt. Buchon (ATOS 501–844) 

5. Southern periphery, Pt. Buchon – Pt. Conception (ATOS 845–1110) 

6. Southern front, South of Pt. Conception (ATOS 1111–1170) 

Note that the ATOS outer boundaries for 1 and 6 correspond to the current range limits 

(2003–2004), and were used to initiate forward simulations. We define the southern and 

northern range boundaries as the two points on the ATOS line spanning 99.5% of the 

spring survey count (allowing for up to 4 outlying individual animals at each end), 

recognizing that this is somewhat arbitrary and that some animals will occasionally be 

observed well beyond these boundaries. Also, we assumed that sea otters in the frontal 

areas would exhibit identical vital rates to those in the neighboring sub-populations 

(although the age/sex structure would be different). Thus vital rates from sub-population 

2 were used to parameterize 1, and vital rates from 5 were used to parameterize 6. 

 To account for uncertainty associated with future population dynamics we used a 

re-sampling approach, utilizing the range of available vital rate estimates to parameterize 

the demographic matrices uniquely for each new simulation. Analyses of carcass age 

structure data provided 10 years (1992-2001) of estimates for each sub-population 

(Chapter 2). Mark-recapture analyses of telemetry provided two more sets of estimates, 

one for the 1980’s (Siniff and Ralls 1988, Siniff and Ralls 1991) and one for 2001–2003 

(Chapter 2). Accumulating evidence suggests that there has been very little variation in 

reproduction parameters over the past 20 years, so we used a single set of age-specific 

rates for all simulations: these were set according to the birth rates and weaning success 

rates calculated from radio-tagged study animals at San Simeon (Chapter 2) and were 

also consistent with values reported in the literature (Siniff and Ralls 1991, Jameson and 
Johnson 1993, Riedman et al. 1994). All of the estimates for demographic parameters 

that we used for simulations are summarized in Table 1. 

 Movement probabilities were parameterized by fitting Laplace probability 

distributions to annual dispersal distances that had been recorded from radio-tagged study 

animals (Figure 3). Raw data were available for 72 study animals from the current study, 

with movements restricted primarily to the south half of the range; these data were 

augmented by data from a concurrent study in the north half of the range (Bodkin and 

Staedler, unpublished data) and from a similar telemetry study in the 1980’s (Siniff and 

Ralls 1988, Ralls et al. 1996). Further details about the collection and analyses of 

movement data can be found in Chapter 3. For this particular analysis we used maximum 

likelihood methods to fit probability distributions for 4 age/sex classes: juvenile/sub-adult 

females, adult females, juvenile/sub-adult males and adult males. The juvenile/sub-adult 

age classes and adult/aged-adult age classes were pooled because there were insufficient 

sample sizes (particularly for juveniles) to allow calculation of separate distributions. 

Dispersal distance kernels were calculated separately for the north half and south half of 

the range, and for the southern range front in the case of males. Preliminary analysis 

indicated that only the data for adult males differed significantly between subpopulations, and so data were pooled across areas for the other age/sex classes. For each 

probability distribution we calculated the 95% confidence intervals around σi,x (Table 1), 

and used this range of values to parameterize movement matrices and dispersal kernel 

functions. 

Running Simulations 

Projecting population dynamics can be accomplished simply by matrix 

multiplication with a population vector. The population vector consists of the number of 

animals in each stage-class, thus the vector length must equal the number of rows in the 

projection matrix: in this case, 8 values (4 stages for each sex) for each sub-population, 

giving a total length of 48. One common way to initialize such a population vector is to 

multiply an estimate of population size (in this case the survey count for 2003) by the 

stable stage distribution (SSD) calculated from the matrix using standard algebraic 

techniques (Caswell 2001). However, this approach requires the assumption that 

demographic rates have been approximately stable for long enough that the age 
distribution has converged on the SSD: in the case of the southern sea otter, there is 

considerable evidence that this has not been the case (Estes et al. 2003a). Consequently, 

prior to running forward simulations (i.e. to project future population growth and range 

expansion), it was necessary to run a historical or “hind-cast” simulation to initialize the 

age-structure for each sub-population. To accomplish this, we utilized the historical 

demographic rates presented in Table 1 to simulate population dynamics from 1989 to 

2003. 

We initialized the 1989 population vectors for each sub-population by multiplying 

the 1989 spring census count by the SSD associated with the 1980’s demographic rates 

(Table 1, estimate 12). Movement matrices were parameterized using the best-fit 

dispersal kernels for each age/sex class. We then projected 14 years of population 

dynamics (Figure 4), calculating all demographic transitions, dispersal and range 

expansion rates as explained above. We adjusted the demographic rates for the 4th–14th

years of the projection (1992 to 2003), setting them to equal the appropriate maximum 

likelihood estimates (Chapter 2; Table 1, estimates 1–11). The outcome of this historical 

projection was an expected population vector for 2003, which was used to initialize stage 

distributions for all forward simulations. An additional result was a comparison of 

expected vs. observed population counts and expected vs. observed range expansion, 

which we used as a way of graphically evaluating the efficacy of our model structure, 

assumptions and parameter values. 

We conducted forward simulations in a similar way, projecting 15 years of 

population dynamics and range expansion using matrix multiplication. We first created 

500 unique dispersal kernels by randomly selecting stage- and location-specific Laplace 

distribution parameters (σi,x) from within the ranges listed in Table 1. For each of the 

resulting 500 movement matrices, we ran 20 simulations using different demographic 

matrices: the first 10 iterations were parameterized using the best-fit maximum likelihood 

values from 1992-2001 (Chapter 2; Table 1, estimates 1–10), while for iterations 11–20 

we randomly selected vital rates from within the 80% confidence intervals associated 

with mark-recapture parameter estimates (Chapter 2; Table 1, estimates 10–11). 

Confidence intervals for each estimate were calculated from standard errors using a logitbased “back transform” method (Burnham and Anderson 1998). The random 
combinations of dispersal and demographic estimates resulted in 10,000 unique iterations 

of the simulation model. 

We summarized simulation results graphically by plotting frequency histograms 

of three key results: the predicted number of independent otters south of Pt. Conception 

(i.e. in sub-population 6) after 10 years, the predicted number after 15 years, and the 

predicted rate of range expansion to the south (in units of km/year). We estimated the 

mean, median, mode and variance for these three variables, as well as their 95% 

confidence limits. To calculate confidence intervals we assumed a negative binomial 

probability distribution in the case of the number of otters south of Pt. Conception, and a 

Weibull probability distribution in the case of the range expansion speed. We also 

estimated summary statistics for the net population growth rate (λ, calculated as the 

geometric mean rate of growth for each simulation), which was normally distributed. 

Finally, we calculated the sensitivity of simulation results to all model parameters 

using multiple regression analysis: specifically, we calculated the proportion of variance 

in three response variables (the predicted number of otters south of Pt. Conception after 

15 years, the rate of southward range expansion, and λ) explained by each of the 

demographic and dispersal parameters, after accounting for variance due to all other 

parameters. Individual variance components were estimated by their partial coefficients 

of determination (r2

p), following (Neter et al. 1990). 
Results 

The historical projection simulation resulted in population dynamics that were 

consistent with observed survey counts over the same period (Figure 4a). While this was 

not especially surprising (the survey counts were one of the data sets used to fit the 

demographic rates, along with carcass data), it nonetheless suggested that the resulting 

stage distribution vector for 2003 was reasonably accurate, and also clearly demonstrates 

the range of different growth rates possible under the simulation parameters. Also 

encouraging was the close agreement between actual southward range expansion over the 

past 15 years and the predicted population wave speed. Although the position of the 

southern range boundary from year-to-year was highly variable, the long-term trend was 

fit by a linear expansion rate of approximately 4.9 km/year (R2

 = 0.59). The mean 

predicted rate of expansion over the same period, as calculated from stage-specific 

dispersal and demographic rates, was 3.95, a value not significantly different from the 

observed trend (Figure 4b). 

The net annual rate of population increase (λ) for all forward simulations was 

1.01, and 95% of the simulations resulted in λ of 0.971–1.052. The rate of population 

growth to the south of Pt. Conception surpassed that of the rest of the population in 

almost all instances, with 95% of the simulations showing a rate of increase south of Pt. 

Conception of 4–20% per year. The median number of independent otters south of Pt. 

Conception after 10 years was 117, and after 15 years this value had increased to 131 

(Figure 5). The rapid growth to the south was partly attributable to dispersal from other 

portions of the population, but also reflected a high intrinsic rate of growth. The 

interaction between dispersal and intrinsic population increase resulted in continued 

range expansion to the south in virtually all simulations: the median predicted wave 

speed was 4.9 km/year over the 15 year projection (Figure 6). Interestingly, this wave 

speed is precisely the same as the average rate of expansion over the past 15 years (Figure 

4). Continued range expansion at this median rate would mean that after 10 years the 

southern range boundary will have moved to a location near Santa Barbara harbor mouth 

(ATOS = 1267), and after 15 years to Carpinteria (ATOS = 1316). There is a great deal 
of uncertainty around these estimates: the 95% confidence interval for the 15 year 

estimate was ATOS = 1183–1584. Table 2 summarizes all simulation statistics. 

The simulation results were sensitive to both variation in dispersal and variation 

in survival parameters, but the relative magnitude of sensitivities was quite different for 

different response variables. Variation in dispersal parameters had the most substantial 

impact on the predicted number of individuals south of Pt. Conception, but had a 

negligible effect on net population growth (Figure 7). Not surprisingly, variation in 

survival rate parameters at the south end of the range had a strong effect on all three 

response variables; however, while variation in survival rates at the center of the range 

had minimal effect on future range expansion and population growth south of Pt. 

Conception, their impact on net population growth was three times greater than survival 

rates at the south end of the range (Figure 7). 

A closer inspection of stage-specific sensitivities showed that, in terms of 

dispersal, juvenile/sub-adult female movement rates had the greatest effect on population 

growth and range expansion to the south (Figure 8). Dispersal of Juvenile/Sub-adult 

males had a significant effect on the expected number of individuals south of Pt. 

Conception, but virtually no effect on the rate of southward range expansion. Adult male 

dispersal had almost no effect on the simulation results, despite the long-distance 

movements frequently conducted by this class of animals. Stage-specific survival rates 

showed a similar pattern of sensitivities: variation in sub-adult female survival had the 

greatest impact on simulation results, while juvenile and adult female survival had less of 

an effect (Figure 9). The only result showing any sensitivity to male survival rates was 

the number of otters south of Pt. Conception, and variation in male survival had almost 

no effect on the rate of range expansion or on net population growth. 

 
Discussion 

The predictions of our hind-cast model closely matched the historical data on 

rates of southward range expansion, suggesting that estimation of asymptotic wave speed 

(Neubert and Caswell 2000) is an appropriate technique for simulating range expansion 

of southern sea otters. This method is particularly appropriate for a population that is 

expanding along a 1-dimensional axis, as is the case with the southern sea otter. Because 

it incorporates information on stage-specific dispersal probabilities, demographic rates 

and population structure, the integrodifference approach is also likely to provide a better 

approximation to range expansion dynamics than the 1-dimensional diffusion model used 

previously to model invasion speed in sea otters (Lubina and Levin 1988). 

Explicit analysis of uncertainty can provide useful insights to managers (Doak 

and Mills 1994, Pascual and Adkison 1994, Ralls and Taylor 2000). The best way to 

incorporate uncertainty into management decisions is to consider, as in our analysis, the 

full range of expected outcomes (Gerber et al. in press). Projections of both the number 

of independent otters at Point Conception and the rate of southward range expansion 

were highly variable, reflecting the uncertainty in input parameter estimates and 

uncertainty about the ultimate causes of fluctuations in survival rates, such as density 

dependence, disease, and fishing interactions. 

Despite this variability, sensitivity analysis of the model’s predictions gave us a 

greatly improved understanding of the processes underlying population growth and range 

expansion. Sensitivity analysis in this case serves two main purposes. First, it identifies 

the parameters to which the model is most sensitive: better estimates of these parameters 

will therefore do most to improve the precision of the model predictions. Our analysis 

identified the survival rate of juvenile and sub-adult females at the end of range as a key 

parameter influencing both population growth and range expansion to the south of Pt. 

Conception. Hence, fieldwork designed to improve estimates of survival rates of young 

females in southern areas would do most to reduce uncertainty in these particular 

predictions. 

Second, sensitivity analysis highlights the particular components of the population 

that are driving range expansion and/or population growth. These results are sometimes 
not intuitively obvious. For example, although males are more likely to move long 

distances than females and most of the individuals that travel south of Point Conception 

are males, male movements proved much less important than female movements (Figure 

8). Movement rates of juvenile and sub-adult females had the greatest effect on both 

population growth and range expansion to south, whereas dispersal of their male 

counterparts had no impact on the rate of southward range expansion, and variation in 

adult male survival had almost no effect on either population growth or range expansion. 

This last result is not so surprising considering that range expansion by males alone 

would provide no intrinsic population growth (i.e. reproduction) at the ends of the range: 

because reproduction is ultimately what dives population growth and subsequent range 

expansion, it is the movement and survival of females that is the limiting factor for both 

processes. This method of sensitivity analysis also allows for evaluation of spatial 

patterns: for instance, female survival at the center of the range probably has little effect 

on the rate of range expansion, but is the most important demographic parameter for 

predicting growth of the population as a whole (Figure 7). 

Our simple multi-state matrix model does not explicitly account for a variety of 

important aspects of sea otter biology and ecology: these include density dependence 

(Laidre et al. 2001), spatial and temporal variation in habitat quality (Doak 1995, Thomas 

and Kunin 1999, Virgl and Messier 2000), responses to ephemeral phenomenon such as 

episodic prey recruitment events (Watt et al. 2000), seasonal reproductive peaks and 

movement patterns (Jameson 1989), and important behavioral characteristics such as 

dietary specializations (Estes et al. 2003b), territoriality (Jameson 1989), contagious 

distribution, and male/female (or age class) segregation at smaller spatial scales. It is 

worth noting that the model actually does implicitly account for some of these factors 

(such as density dependence and habitat quality), in so far as these factors have affected 

past and present vital rates and movement probabilities within the existing range. 

Despite the above-mentioned limitations, our model provides a robust and 

generalizable approach to understanding and predicting population dynamics in southern 

sea otters by making use of all existing demographic and dispersal data. It represents a 

unique synthesis of a multi-state dispersal matrix and the integrodifference equation 

approach to calculating invasion speed. Our model should provide a useful and flexible 
tool for conservation biologists and managers, and can be easily expanded upon or 

improved as additional data and more precise parameter estimates for southern sea otters 

become available. 
References 

Burnham, K. P., and D. R. Anderson. 1998. Model selection and inference: a practical 

information-theoretic approach, 2nd edition. Springer-Verlag, New York, NY. 

Caswell, H. 2001. Matrix population models: construction, analysis, and interpretation, 

2nd ed edition. Sinauer Associates, Sunderland, MA. 

Doak, D. F. 1995. Source-sink models and the problem of habitat degradation: general 

models and applications to the Yellowstone grizzly. Conservation Biology 

9:1370-1379. 

Doak, D. F., and L. S. Mills. 1994. A useful role for theory in conservation. Ecology 

75:615-626. 

Doak, D. F., and W. F. Morris. 2002. Quantitative conservation biology: theory and 

practice of population viability analysis. Sinauer Associates, Sunderland, MA. 

Estes, J. A., B. B. Hatfield, K. Ralls, and J. Ames. 2003a. Causes of mortality in 

California sea otters during periods of population growth and decline. Marine 

Mammal Science 19:198-216. 

Estes, J. A., M. L. Riedman, M. M. Staedler, M. T. Tinker, and B. E. Lyon. 2003b. 

Individual variation in prey selection by sea otters: Patterns, causes and 

implications. Journal of Animal Ecology 72:144-155. 

Gerber, L. R., T. Tinker, D. Doak, and J. Estes. in press. Mortality sensitivity in life-stage 

simulation analysis: A case study of southern sea otters. Ecological Applications. 

Jameson, R. J. 1989. Movements, home range, and territories of male sea otters off 

central California. Marine Mammal Science 5:159-172. 

Jameson, R. J., and A. M. Johnson. 1993. Reproductive characteristics of female sea 

otters. Marine Mammal Science 9:156-167. 

Laidre, K. L., R. J. Jameson, and D. P. DeMaster. 2001. An estimation of carrying 

capacity for sea otters along the California coast. Marine Mammal Science 

17:294-309. 

Lebreton, J. D., and G. Gonzalez-Davila. 1993. An introduction to models of subdivided 

populations. Journal of Biological Systems 1:389-423. 

Lubina, J. A., and S. A. Levin. 1988. The spread of a reinvading species: Range 

expansion in the California sea otter. American Naturalist 131:526-543. 
Neter, J., W. Wasserman, and M. H. Kutner. 1990. Applied linear statistical models: 

Regression, analysis of variance, and experimental designs, 3rd edition. Irwin, 

Chicago, IL. 

Neubert, M. G., and H. Caswell. 2000. Demography and dispersal: Calculation and 

sensitivity analysis of invasion speed for structured populations. Ecology 

81:1613-1628. 

Pascual, M. A., and M. D. Adkison. 1994. The decline of the Steller sea lion in the 

northeast Pacific: demography, harvest or environment? Ecological Applications 

4:393-403. 

Ralls, K., T. C. Eagle, and D. B. Siniff. 1996. Movement and spatial use patterns of 

California sea otters. Canadian Journal of Zoology 74:1841-1849. 

Ralls, K., and B. L. Taylor. 2000. Better policy and management decisions through 

explicit analysis of uncertainty: new approaches from marine conservation 

biology. Conservation Biology 14:1240-1242. 

Riedman, M. L., J. A. Estes, M. M. Staedler, A. A. Giles, and D. R. Carlson. 1994. 

Breeding patterns and reproductive success of California sea otters. Journal of 

Wildlife Management 58:391-399. 

Siniff, D. B., and K. Ralls. 1988. Population status of California sea otters: MMS 14-12-

001-30033. U.S. Department of Interior, Minerals Management Service, 

California. 

Siniff, D. B., and K. Ralls. 1991. Reproduction, survival and tag loss in California sea 

otters. Marine Mammal Science 7:211-229. 

Thomas, C. D., and W. E. Kunin. 1999. The spatial structure of populations. Journal of 

Animal Ecology 68:647-657. 

USFWS. 2003. Final revised recovery plan for the southern sea otter (Enhydra lutris 

nereis). U.S. Fish and Wildlife Service. Portland, OR. 

Virgl, J. A., and F. Messier. 2000. Assessment of source-sink theory for predicting 

demographic rates among habitats that exhibit temporal changes in quality. 

Canadian Journal of Zoology 78:1483-1493. 

Watt, J., D. B. Siniff, and J. A. Estes. 2000. Inter-decadal patterns of population and 

dietary change in sea otters at Amchitka Island, Alaska. Oecologia 124:289-298. 
Tables 

Table 1. Parameter estimates used for the simulation model. Numbers in parentheses (following 

 mean values) are standard errors, while two numbers seperated by a hyphen indicate the range of

 values used in simulation runs.

Model Parameter Juveniles Sub-adults Adults Aged Adults

Annual birth rates 0 0.4 0.98 0.9

Wean success rates 0 0.4 (0.10) 0.061 (0.07) 0.8 (0.07)

Female annual survival rates

 Estimates 1-10: see Ch. 2, Appendix B 

 sub-populations 1-2 0.838 - 0.858 0.847 - 0.869 0.843 - 0.870 0.509 - 0.556

 sub-population 3 0.833 - 0.853 0.842 - 0.864 0.836 - 0.866 0.504 - 0.550

 sub-population 4 0.847 - 0.867 0.854 - 0.876 0.847 - 0.874 0.505 - 0.554

 sub-populations 5-6 0.848 - 0.869 0.856 - 0.878 0.849 - 0.876 0.508 - 0.557

 Estimate 11: 2001-2003 1

 sub-populations 1-2, 5-6 0.85 (0.145) 0.88 (0.145) 0.91 (0.088) 0.55 (n/a)

 sub-populations 3-4 0.84 (0.060) 0.84 (0.060) 0.84 (0.060) 0.55 (n/a)

 Estimate 12: 1980's, all sub populations 0.85 (0.145) 0.88 (0.145) 0.91 (0.088) 0.55 (n/a)

Male annual survival rates

 Estimates 1-10: see Ch. 2, Appendix B 

 sub-populations 1-2 0.782 - 0.809 0.782 - 0.811 0.746 - 0.784 0.328 - 0.371

 sub-population 3 0.776 - 0.802 0.776 - 0.805 0.739 - 0.778 0.322 - 0.365

 sub-population 4 0.793 - 0.820 0.791 - 0.821 0.751 - 0.791 0.324 - 0.370

 sub-populations 5-6 0.795 - 0.822 0.794 - 0.823 0.754 - 0.793 0.327 - 0.373

 Estimate 11: 2001-2003 1

 sub-populations 1-2, 5-6 0.88 (0.179) 0.88 (0.179) 0.87 (0.095) 0.35 (n/a)

 sub-populations 3-4 0.88 (0.179) 0.88 (0.179) 0.84 (0.060) 0.35 (n/a)

 Estimate 12: 1980's, all sub populations 0.88 (0.179) 0.88 (0.179) 0.70 (0.167) 0.35 (n/a)

Laplace Dispersal Parameters (σ) 2

 Females, all sub-populations 32.6 - 83.1 32.6 - 83.1 7.5 - 11.7 7.5 - 11.7

 Males, sub-population 1–3 63.2 - 171.0 63.2 - 171.0 7.6 - 20.6 7.6 - 20.6

 Males, sub-population 4 63.2 - 171.0 63.2 - 171.0 39.5 - 95.9 39.5 - 95.9

 Males, sub-populations 5–6 63.2 - 171.0 63.2 - 171.0 79.7 - 150.0 79.7 - 150.0

1 Estimates correspond to 1980-'s values for locations or stages not measured in 2001-2003

2 Units = 500m increments (ATOS values)

Table 2. Summary of results from simulations

Variable Mean Std. dev. Median Mode L95 U95

Net rate of increase (λ) 1.011 0.021 1.012 1.01 0.971 1.052

S. of Pt. Conception, 10yrs 120.69 38.477 117 107 57 207

S. of Pt. Conception, 15yrs 136.63 53.212 131 121 53 259

Southern Exp. Speed(km/yr) 5.0679 3.3469 4.86 1.43 0.422 13.82

Range End, 10 yrs 1271 1267 1199 1178 1446

Range End, 15 yrs 1322 1316 1213 1183 1584
Figure Captions 

Figure 1. Map of central California showing current range of the southern sea otter (excluding 

San Nicolas Island), and identifying the spatial arrangement of the six sub-populations identified 

for the simulation model. 

Figure 2. The spatial distribution of sea otters along the coast (on a north-to-south axis) is 

plotted as a histogram for each of the six sub-populations (see Figure 1), based on the results of 

the five most recent spring surveys. Beta probability density functions were fit to each data set 

and are overlain on the histograms: the parameter values for each function are displayed above 

each curve. The Beta functions were used in the calculation of movement rates between subpopulations (see text). 

Figure 3. Annual dispersal distance frequency histograms are shown for adult females (pink 

bars) and juvenile females (blue bars). Laplace probability functions were fit to each of these 

distributions (dashed and solid lines), and used in the calculation of stage-specific movement 

rates. Note that the distribution for juvenile females shows greater dispersion, which is reflected 

by a higher value of the scale parameter σ. Similar functions were calculated for males (not 

shown here). 

Figure 4. Results of a historical simulation of population dynamics for the southern sea otter 

population over the years 1989–2003. Predicted population counts, based on the simulation, are 

shown at top, with observed counts for comparison. Predicted range expansion to the south 

(increasing ATOS values over time) is shown at bottom, with observed range-end boundaries 

shown for comparison. The range end was defined as the point along the coast at which 99.75% 

of the sea otter population was to the north, based on the annual spring survey. A linear leastsquares curve was fit to the observed range-end dataset, and is plotted (along with the 95% 

prediction interval) to illustrate the correspondence between the predicted and observed mean rate 

of expansion. 

Figure 5. A frequency distribution of predicted outcomes is shown for two of the key simulation 

results: the expected number of independent sea otters south of Pt. Conception after 10 years (top) 

and after 15 years (bottom). The distributions were well described by negative-binomial 

probability distributions. 

Figure 6. A frequency distribution of the predicted rate of southward range expansion is shown, 

based on the results of 10,000 replicate simulations. A weibull probability distribution was fit to 

the raw data. 

Figure 7. Results of a sensitivity analysis, showing the relative proportion of the variance in 

simulation dynamics explained by three groups of model parameters: dispersal rates, survival 

rates at the south of the range (sub-populations 5 and 6) and survival rates at the center of the 

range (sub-populations 3 and 4). Sensitivities are shown for three response variables: A) the 

number of individual otters south of pt. conception after 15 years, B) average southward wave 

speed, or rate of range expansion to the south, and C) the overall rate of population increase over 

the simulation period. 

Figure 8. Results of a sensitivity analysis of dispersal parameters, showing the relative 

proportion of the variance in simulation dynamics explained by stage- and sex-specific dispersal 

rates. Sensitivities are shown for two response variables: A) the number of individual otters south 

of pt. conception after 15 years, B) average southward wave speed, or rate of range expansion to 

the south. 
Figure 9. Results of a sensitivity analysis of survival parameters, showing the relative proportion 

of the variance in simulation dynamics explained by stage- and sex-specific survival rates. 

Sensitivities are shown for three response variables: A) the number of individual otters south of 

pt. conception after 15 years, B) average southward wave speed, or rate of range expansion to the 

south, and C) the overall rate of population increase over the simulation period. 
 Figure 1 

120o W

30o N

Study 

Area Santa Cruz

Monterey Bay

Pt. Sur

Pt. Buchon

Pt. Conception

Santa Barbara

San Francisco Bay

0 40 80 120 160 20 Kilometers

Sea Otter Range, 2003

Sub-populations

1. Northern Front

2. Northern Periphery

3. North-Central

4. South-Central

5. Southern Periphery

6. Southern Front

120o W

30o N

Study 

Area Santa Cruz

Monterey Bay

Pt. Sur

Pt. Buchon

Pt. Conception

Santa Barbara

San Francisco Bay

0 40 80 120 160 20 Kilometers

Sea Otter Range, 2003

Sub-populations

1. Northern Front

2. Northern Periphery

3. North-Central

4. South-Central

5. Southern Periphery

6. Southern Front
 Figure 2 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Sub-population 1, β param. a = 5.45, b = 3.78

Sub-population 2, β param. a = 1.15, b = 0.529

Sub-population 3, β param. a = 1.2, b = 1.16

Sub-population 4, β param. a = 1.32, b = 1.17

Sub-population 5, β param. a = 0.739, b = 1.38

Sub-population 6, β param. a = 2.15, b = 5.76

Relative position along coastline

North boundary South boundary

Relative density of sea otters (1999-2003 spring surveys) 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Sub-population 1, β param. a = 5.45, b = 3.78

Sub-population 2, β param. a = 1.15, b = 0.529

Sub-population 3, β param. a = 1.2, b = 1.16

Sub-population 4, β param. a = 1.32, b = 1.17

Sub-population 5, β param. a = 0.739, b = 1.38

Sub-population 6, β param. a = 2.15, b = 5.76

Relative position along coastline

North boundary South boundary

Relative density of sea otters (1999-2003 spring surveys) 
 Figure 3 

-30 -20 -10 0 10 20 30 40

0

0.02

0.04

0.06

0.08

0.1

Annual female dispersal distance (km)

Relative frequency

Laplace distribution, juveniles

Laplace distribution, adults

-30 -20 -10 0 10 20 30 40

0

0.02

0.04

0.06

0.08

0.1

Annual female dispersal distance (km)

Relative frequency

Laplace distribution, juveniles

Laplace distribution, adults
 Figure 4 

1600

1800

2000

2200

2400

Counts of Independent otters 

Observed Counts

Expected Counts

950

1000

1050

1100

1150

1200

1988 1990 1992 1994 1996 1998 2000 2002 2004

Year

Southern range boundary (ATOS value) 

Observed Range End (3yr Avg.)

Predicted Range End (Hind-cast Simulation)

Observed Linear Trend (R2 = 0.592)

95% Regression Prediction Interval

1600

1800

2000

2200

2400

Counts of Independent otters 

Observed Counts

Expected Counts

1600

1800

2000

2200

2400

Counts of Independent otters 

Observed Counts

Expected Counts

950

1000

1050

1100

1150

1200

1988 1990 1992 1994 1996 1998 2000 2002 2004

Year

Southern range boundary (ATOS value) 

Observed Range End (3yr Avg.)

Predicted Range End (Hind-cast Simulation)

Observed Linear Trend (R2 = 0.592)

95% Regression Prediction Interval

950

1000

1050

1100

1150

1200

1988 1990 1992 1994 1996 1998 2000 2002 2004

Year

Southern range boundary (ATOS value) 

Observed Range End (3yr Avg.)

Predicted Range End (Hind-cast Simulation)

Observed Linear Trend (R2 = 0.592)

95% Regression Prediction Interval

Observed Range End (3yr Avg.)

Predicted Range End (Hind-cast Simulation)

Observed Linear Trend (R2 = 0.592)

95% Regression Prediction Interval
 Figure 5 

200

400

600

800

1000

1200

1400

Ten Year Projection

0 50 100 150 200 250 300 350 400 0 

200

400

600

800

1000

1200

1400

Number of independent otters south of Pt. Conception

Number of simulations

Fifteen Year Projection

200

400

600

800

1000

1200

1400

Ten Year Projection

0 50 100 150 200 250 300 350 400 0 

200

400

600

800

1000

1200

1400

Number of independent otters south of Pt. Conception

Number of simulations

Fifteen Year Projection
 Figure 6 

0 2 4 6 8 10 12 14 16 18 20

0

200

400

600

800

1000

1200

Mean population wave speed (km/year)

Number of simulations

Rate of Southward Range Expansion

0 2 4 6 8 10 12 14 16 18 20

0

200

400

600

800

1000

1200

Mean population wave speed (km/year)

Number of simulations

Rate of Southward Range Expansion
 Figure 7 

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

Dispersal distance Survival, south end

of range

Relative sensitivity (r2

p coefficient)

A. Number of individual animals south of Pt. Conception

B. Rate of southward range expansion

C. Range-wide rate of population growth (λ) 

Survival, center

of range

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

Dispersal distance Survival, south end

of range

Relative sensitivity (r2

p coefficient)

A. Number of individual animals south of Pt. Conception

B. Rate of southward range expansion

C. Range-wide rate of population growth (λ) 

Survival, center

of range
 Figure 8 

B. Rate of Southward Range Expansion

0

0.05

0.1

0.15

0.2

Juv/Sub-A. Fm. Adult Fm.

A. Number of individual animals south of Pt. Conception

0

0.04

0.08

0.12

0.16

0.2

Relative sensitivity to dispersal (r2

p coefficient)

Juv/Sub-A. Ma. Adult Ma.

B. Rate of Southward Range Expansion

0

0.05

0.1

0.15

0.2

Juv/Sub-A. Fm. Adult Fm.

A. Number of individual animals south of Pt. Conception

0

0.04

0.08

0.12

0.16

0.2

A. Number of individual animals south of Pt. Conception

0

0.04

0.08

0.12

0.16

0.2

Relative sensitivity to dispersal (r2

p coefficient)

Juv/Sub-A. Ma. Adult Ma.
 Figure 9 

A. Number of individual animals south of Pt. Conception

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

B. Rate of Southward Range Expansion

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

Relative sensitivity to survival (r2

p coefficient)

C. Range-wide Rate of Population Growth (λ) 

0

0.02

0.04

0.06

0.08

0.1

0.12

Juv. Fm. Sub-A. Fm. Adl. Fm. Juv. Ma. Sub-A. Ma. Adl. Ma.

A. Number of individual animals south of Pt. Conception

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

B. Rate of Southward Range Expansion

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

Relative sensitivity to survival (r2

p coefficient)

C. Range-wide Rate of Population Growth (λ) 

0

0.02

0.04

0.06

0.08

0.1

0.12

Juv. Fm. Sub-A. Fm. Adl. Fm. Juv. Ma. Sub-A. Ma. Adl. Ma.