Population dynamics of walruses in Greenland

The historical and current dynamics of the three Atlantic walrus (Odobenus rosmarus rosmarus) populations that occur in Greenland are estimated using ageand sex-structured population models with exponential growth, density-regulated growth and selection-delayed dynamics. These models are integrated with data in a Bayesian framework, where the likelihood of the simulated population trajectories are evaluated from recent abundance estimates and age-structure information from a selective hunt. The overall decline in the Baffin Bay population caused by historical catches is unclear due to incomplete catch reporting prior to 1950s. However, it is estimated that the population declined by 40% from the 1960s to 2005; decreased catches (≈140 to ≈70) have subsequently allowed this population to increase. The 2012 abundance estimate is 1,400 (95% CI: 1,000–2,000) individuals, and the annual natural growth rate in this population is now 7.7% (95% CI: 6.7–8.9%). Averaging across models, it is estimated that West Greenland/Baffin Island walruses declined by 80% from 7,000 (95% CI: 5,400–10,000) in 1900 to 1,350 (CI: 950–1,950) in 1960. Hereafter they increased to 3,100 (95% CI: 2,500–4,400) in 1993, and owing to increased catches they have experienced a minor decline between 1994 and the early 2000s. Annual catches were then cut from 190 to the current quota of 61, and the population is again increasing with a 2012 estimate of 3,900 (95% CI: 2,500–5,300) individuals. A 2012 estimate of 1,400 (95% CI: 700–3,100) walruses in East Greenland is recovered relative to 1888; the year prior to our first historical catches by European sealers. The historical trajectory, however, is uncertain. Density regulation estimates a relatively flat trajectory, with a maximal depletion in 1890 to 80% of the initial abundance, and a slow continuous increase to almost no current growth. A recovered population is also estimated by selection-delayed dynamics. However, this model estimates a continued increase, and a historical depletion to 2% in 1957. These results are only partially comparable with an earlier assessment for East Greenland. Updated abundance estimates for West Greenland, and modelling with age-structured data from Baffin Bay, have greatly improved our status estimates for Baffin Bay and West Greenland/Baffin Island. Herein, we perform population dynamic modelling based on the historical harvest to enhance our understanding of walrus population dynamics in Greenland. We aim to provide scientific harvest recommendations to the North Atlantic Marine Mammal Commission (NAMMCO) and, thus, it is essential that the uncertainties of our models are quantified (Harwood and Stokes 2003; Clark and Bjørnstad 2004). This is achieved by integrating abundance and age-structured harvest data with population dynamic models embedded in a Bayesian framework. Our analyses incorporate parameter uncertainty, observation errors and model uncertainty, but we do not consider demographic and environmental stochasticity, which is often referred to as process variation. Among the four sources of uncertainty, model uncertainty has generally received the least attention (Buckland et al. 2007). This is unfortunate because, as long as demography is not seriously affected by random variation, it is the underlying population dynamic mechanisms that define the population trajectories. If population models are wrong, the estimated trajectories and their associated uncertainty cannot be trusted, even when the statistical analyses are conducted in an ideal manner. Single species models of animal population dynamics have usually assumed either exponential or density-regulated growth. However, recent findings show that population dynamics and evolutionary changes often occur on similar time-scales (Thompson 1998; Hairston et al. 2005; Saccheri and Hanski 2006; Coulson et al. 2011; Schoener 2011), so it has become increasingly essential to explore the population dynamic consequences of evolutionary processes. So far, there has been a general interest in the population dynamic consequences of evolution caused by fishery-induced selection (e.g., Law 2000; Browman 2000; Carlson et al. 2007; Edeline et al. 2007; Bromaghin et al. 2011; Scott and Sampson 2011). However, ordinary evolution occurs by natural selection and not by anthropogenic selection, and Sinervo et al. (2000) have demonstrated that natural selection is a sufficient mechanism for population cycles in side-blotched lizards (Sinervo et al. 2000), as suggested originally by Voipio (1950) and Chitty (1960). Removal of individuals by harvest does not induce direct anthropogenic selection unless the harvest is selective on some phenotypic trait. Nevertheless, even a non-selective harvest will induce natural selection because it changes the density-dependent selection pressure within the population. This has been investigated across seven species of large baleen whales (suborder Mysticeti) that were heavily depleted by past commercial whaling (Witting 2013). By applying age-structured population models of density-regulated growth, density-regulated growth with depensation, predation by killer whales and selection-delayed dynamics, it was found that only the latter mechanism sufficiently reconciled current rates of increase with past levels of exploitation. Selection-delayed dynamics have a density-regulated exponential growth rate that is accelerated and decelerated via intra-specific selection by density-dependent competitive interactions (Witting 2000a,b; 2002; 2013). Depending upon the selection response of the population, it results in a damped to stable population cycle that is associated with a cycle in one or several life history parameters. The body mass, e.g., is expected to be largest during the late peak and early declining phases of a cycle and smallest during the late low and early increasing phases (see prediction in Fig. 4 in Witting 2013, and Daphnia data in Fig. 6 in Witting 2000b). Such life history cycles are widespread in species with population dynamic cycles (e.g., Chitty 1952; Hansson 1969; Krebs and Myers 1974; Boonstra and Krebs 1979; Mihok et al. 1985; Myers 1990; Lidicker and Ostfeld 1991; Stenseth and Ims 1993; Watson et al. 1994; Hodges et al. 1999; Simchuk et al. 1999; Sinervo et al. 2000; Ergon et al. 2001; Norrdahl and Korpimäki 2002; Lambin et al. 2006) and, as they generally are unexplained in the absence of density-frequency-dependent selection, they indicate that several natural populations may be driven by selection-delayed dynamics. So, to explore the range of natural single species dynamics we apply short-term modWalrus of the North Atlantic 192 els of exponential growth, as well as longer-term models of density-regulated growth and selection-delayed dynamics. The models developed in this study for walruses are age and sex-structured, and the Bayesian integration between our models and data resemble hidden process models (Newman et al. 2006; Buckland et al. 2007). Hence, our data do not directly describe the different states of each process, but only overall estimates at certain points in time (i.e., total abundance through surveys and the age-structure of a selective hunt). Our models produce trajectories given their assumptions and prior distributions of the parameters, and the likelihood of these time-series of theoretical truths are subsequently evaluated by comparison with our observed time-series of abundance and age-structure data. A population dynamic assessment based on comparable methods in 2005 indicated that two of the three walrus populations in Greenland were likely over-exploited (Witting and Born 2005). While the assessment was based on the best available data at that time, our knowledge of these populations and our population modelling abilities have increased since that time. Based on satellite tracking, it is now clear that walruses in West Greenland are part of a larger population that summers primarily around Baffin Island (Dietz et al. 2009). And where the assessment in 2005 was based primarily on best-guess-estimates of abundance, we now have more solid abundance estimates for all these populations (Heide-Jørgensen et al. 2014; Born et al. 2009; Heide-Jørgensen et al. 2009b; Stewart et al. 2009a). Our current population models can also be fitted to age-structured data from a selective hunt, and this provides better estimates of the population dynamic growth rate. This study can thus be seen as an update of the 2005 assessment and an improvement in the estimates on the status of the three populations. The major data uncertainties in this paper are associated with abundance estimates and catch histories, including unknown loss and reporting rates. Many of the historical catches for Baffin Bay are model constructed and, thus, associated with a larger degree of uncertainty than the catch histories from the two other populations. Our description of the dynamics rest on the available data and our population dynamic models, and they do not provide a complete description of the dynamic history.


INTRODUCTION
The three populations of Atlantic walrus (Odobenus rosmarus rosmarus) that occur in Greenland have been subject to exploitation for centuries. First at a limited scale by Greenlanders, and then by European whalers and sealers who heavily impacted these populations. From the beginning of the 20 th Century Greenlanders have hunted walruses with increasing efforts; since the introduction of fire-arms and motorized vessels. The populations are still exploited for subsistence purposes in Greenland and Canada, with an additional small scale sport hunt in Canada.
Herein, we perform population dynamic modelling based on the historical harvest to enhance our understanding of walrus population dynamics in Greenland. We aim to provide scientific harvest recommendations to the North Atlantic Marine Mammal Commission (NAMMCO) and, thus, it is essential that the uncertainties of our models are quantified (Harwood and Stokes 2003;Clark and Bjørnstad 2004). This is achieved by integrating abundance and age-structured harvest data with population dynamic models embedded in a Bayesian framework. Our analyses incorporate parameter uncertainty, observation errors and model uncertainty, but we do not consider demographic and environmental stochasticity, which is often referred to as process variation.
Among the four sources of uncertainty, model uncertainty has generally received the least attention (Buckland et al. 2007). This is unfortunate because, as long as demography is not seriously affected by random variation, it is the underlying population dynamic mechanisms that define the population trajectories. If population models are wrong, the estimated trajectories and their associated uncertainty cannot be trusted, even when the statistical analyses are conducted in an ideal manner.
Single species models of animal population dynamics have usually assumed either exponential or density-regulated growth. However, recent findings show that population dynamics and evolutionary changes often occur on similar time-scales (Thompson 1998;Hairston et al. 2005;Saccheri and Hanski 2006;Coulson et al. 2011;Schoener 2011), so it has become increasingly essential to explore the population dynamic consequences of evolutionary processes. So far, there has been a general interest in the population dynamic consequences of evolution caused by fishery-induced selection (e.g., Law 2000;Browman 2000;Carlson et al. 2007;Edeline et al. 2007;Bromaghin et al. 2011;Scott and Sampson 2011). However, ordinary evolution occurs by natural selection and not by anthropogenic selection, and Sinervo et al. (2000) have demonstrated that natural selection is a sufficient mechanism for population cycles in side-blotched lizards (Sinervo et al. 2000), as suggested originally by Voipio (1950) and Chitty (1960).
Removal of individuals by harvest does not induce direct anthropogenic selection unless the harvest is selective on some phenotypic trait. Nevertheless, even a non-selective harvest will induce natural selection because it changes the density-dependent selection pressure within the population. This has been investigated across seven species of large baleen whales (suborder Mysticeti) that were heavily depleted by past commercial whaling (Witting 2013). By applying age-structured population models of density-regulated growth, density-regulated growth with depensation, predation by killer whales and selection-delayed dynamics, it was found that only the latter mechanism sufficiently reconciled current rates of increase with past levels of exploitation.
Selection-delayed dynamics have a density-regulated exponential growth rate that is accelerated and decelerated via intra-specific selection by density-dependent competitive interactions (Witting 2000a,b;2013). Depending upon the selection response of the population, it results in a damped to stable population cycle that is associated with a cycle in one or several life history parameters. The body mass, e.g., is expected to be largest during the late peak and early declining phases of a cycle and smallest during the late low and early increasing phases (see prediction in Fig. 4 in Witting 2013, andDaphnia data in Fig. 6 in Witting 2000b). Such life history cycles are widespread in species with population dynamic cycles (e.g., Chitty 1952;Hansson 1969;Krebs and Myers 1974;Boonstra and Krebs 1979;Mihok et al. 1985;Myers 1990;Lidicker and Ostfeld 1991;Stenseth and Ims 1993;Watson et al. 1994;Hodges et al. 1999;Simchuk et al. 1999;Sinervo et al. 2000;Ergon et al. 2001;Norrdahl and Korpimäki 2002;Lambin et al. 2006) and, as they generally are unexplained in the absence of density-frequency-dependent selection, they indicate that several natural populations may be driven by selection-delayed dynamics. So, to explore the range of natural single species dynamics we apply short-term mod-Walrus of the North Atlantic els of exponential growth, as well as longer-term models of density-regulated growth and selection-delayed dynamics.
The models developed in this study for walruses are age and sex-structured, and the Bayesian integration between our models and data resemble hidden process models (Newman et al. 2006;Buckland et al. 2007). Hence, our data do not directly describe the different states of each process, but only overall estimates at certain points in time (i.e., total abundance through surveys and the age-structure of a selective hunt). Our models produce trajectories given their assumptions and prior distributions of the parameters, and the likelihood of these time-series of theoretical truths are subsequently evaluated by comparison with our observed time-series of abundance and age-structure data.
A population dynamic assessment based on comparable methods in 2005 indicated that two of the three walrus populations in Greenland were likely over-exploited (Witting and Born 2005). While the assessment was based on the best available data at that time, our knowledge of these populations and our population modelling abilities have increased since that time. Based on satellite tracking, it is now clear that walruses in West Greenland are part of a larger population that summers primarily around Baffin Island ). And where the assessment in 2005 was based primarily on best-guess-estimates of abundance, we now have more solid abundance estimates for all these populations (Heide-Jørgensen et al. 2014;Born et al. 2009;Heide-Jørgensen et al. 2009b;Stewart et al. 2009a). Our current population models can also be fitted to age-structured data from a selective hunt, and this provides better estimates of the population dynamic growth rate. This study can thus be seen as an update of the 2005 assessment and an improvement in the estimates on the status of the three populations.
The major data uncertainties in this paper are associated with abundance estimates and catch histories, including unknown loss and reporting rates. Many of the historical catches for Baffin Bay are model constructed and, thus, associated with a larger degree of uncertainty than the catch histories from the two other populations. Our description of the dynamics rest on the available data and our population dynamic models, and they do not provide a complete description of the dynamic history.

Stock structure
There are three populations of Atlantic walrus in Greenland. West Greenland/Baffin Island walruses occur from fall to spring in the central part of West Greenland at the edge of the pack ice from c. 66°30 N to 70°30 N (Born et al. 1994;Born et al. 1995). During the open water season they are found mainly around the Southeastern Baffin Island. Further north in Baffin Bay and Smith Sound, walruses occur almost year-round in the North Water polynya and adjacent areas. They are, however, absent from the coastal areas of Northwest Greenland in August and September when they summer along the eastern and southern coast of Ellesmere Island (Canada) and in the Canadian High Arctic Archipelago (Born et al. 1995). We refer to walruses in these areas as the Baffin Bay population. Walruses occur year-round along the eastern coast of Greenland where they mainly are distributed inside the National Park of North and Northeast Greenland, north of the entrance to Scoresby Sound (c. 71°N; Born et al. 1995Born et al. , 1997. There is only limited exchange between East Greenland walruses and neighbouring populations, i.e., West Greenland, Baffin Bay and the Svalbard-Franz Joseph Land populations (Born et al. 1995;Andersen et al. 1998;Born et al. 2001). Information on distribution and migration (Born et al. 1994(Born et al. , 1995(Born et al. , 1997 and genetics (Andersen et al. 1998;Andersen and Born 2000;Born et al. 2001) indicates that the three aggregations of walruses in Greenland represent separate population units and therefore should be managed separately.

NAMMCO Scientific Publications, Volume 9
Abundance data Baffin Bay Based on aerial counts at terrestrial haul-outs, a multiplication factor for walruses in the water, and a multiplication factor for areas not covered, an estimate of 1,500 walruses was applied in the 2005 assessment for the population of walruses in Baffin Bay. Two fully corrected strip-census surveys (corrected for perception bias and availability bias) of the spring abundance in the North Water estimated 1,240 (cv:0.19) walruses in 2009and 1,760 (cv:0.29) in 2010, respectively (Heide-Jørgensen et al. 2009bHeide-Jørgensen, unpublished). We use both estimates in our model to represent the number of walruses in the Baffin Bay population.

West Greenland/Baffin Island
The assessment in 2005 was based on a fully corrected abundance estimate of 938 (cv:0.48) wintering walruses off Central West Greenland in 1990 (Witting and Born 2005). There have been two additional surveys in this areas since that time, in 2006 and 2008 with fully corrected strip census estimates of 2, 790 (cv:0.54) and 3, 240 (cv:0.76;Heide-Jørgensen et al. 2014). On the Canadian side, haul-out counts from Baffin Island give a minimum estimate of 1,000 walruses in 2007 (Stewart et al. 2009b), which increases to 3, 030 (cv:0.20) when corrected for walruses not hauled out.
Apart from the corrected abundance estimates there is a time-series of the relative abundance of wintering walruses off Central West Greenland based on sighting rates during 11 surveys conducted from 1981 to 2008 (Heide-Jørgensen et al. 2009a). This times-series is given both for a northern (west of Disko Island) and a southern (south of Disko Bay) wintering ground in Central West Greenland, with the relative abundance on the northern ground being correlated with sea ice and the relative abundance on the southern ground being independent of ice. A NAMMCO scientific working group concluded in 2009 that these time-series do not provide reliable information on trend (NAMMCO 2010), and we do not use these data in our modelling.

East Greenland
Based on opportunistic and systematic observations, the East Greenland population of walruses was estimated to number 1,000 individuals in the 2005 assessment. Given an aerial survey of the coast from Clavering Island to the northern border of the Northeast Water, there is now a fully corrected summer estimate of 1,430 (cv:0.45) walruses in East Greenland in 2009 ).
All abundance estimates that are used in the modelling are listed in Table 1.

Catch data
We estimated two catch histories for each population. A low catch history based on reported or landed catch (and in some years estimates of the catch), and a high catch history that also included estimates of non-reported losses, i.e., animals that were struck and lost, and animals that were landed and not reported. The catch histories are shown in Figure 1, and details on their estimation are given separately for the three areas in the sub-sections below.

Baffin Bay
Generally, the catch data from the Thule area of Northwest Greenland are insufficient, in particular prior to 1950 (Teilmann and Kapel 1998). To obtain estimates for the early period, catches were inferred from the trend in growth of the human population in the area. The Inuit population in the Thule area increased gradually with a rate of ≈0.8% per year from about 200 at 1900 to about 300 around 1950 (see Fig. 19 in Gilberg 1976). If assuming a proportional relationship between (a) the size of the human population, (b) the fraction of hunters [29-34% were men aged 15-64 years; Gilberg (1976:    od where walruses hauling out on ice flakes were shot at a distance before they are harpooned. After 1970, an overall loss rate of 30% was used based on observations in the late 1970s (Born et al. 1995).
The sex ratio (F: n=179; M: n=197) in an age-structured sample of the North Water catch did not differ from unity (2=0.43, p=0.512, df=1). Hence, we apply an even sex ratio to all the catches from this population. Post 2010, catches were set to the current quota.
West Greenland/Baffin Island Catch data were extracted from official catch statistics and written sources (Born et al. 1994;COSEWIC 2006). The catch history of our model includes all West Greenland catches south of Qaanaaq, plus Canadian catches after 1977 from Qikiqtarjuaq, Iqaluit and Pangnirtung.
Losses were assumed to average 5% (one lost for every five killed) from 1900 until 1930. During this period walruses were mainly taken by traditional means close to the coast and primarily at the terrestrial haul-outs. From 1930 onwards when walruses increasingly were hunted by use of motorized vessels operating in the offshore pack ice (Born et al. 1994) losses were assumed to average 30%. Loss rates of 20-30% are not uncommon during walrus hunts based on small vessels (cf. Born et al. 1995 and references therein;Gjertz et al. 1998).
An even sex ratio was applied to all historical catches. Post 2010, catches were set to the current quota.

East Greenland
Catches for our East Greenland model were extracted from Born et al. (1997). An average loss rate of 30% was applied to the catch by the European sealing vessels, while the loss rate of Greenlanders was set to 5% until 1949, after which a loss rate of 23% was applied (Born et al. 1997). Information is not available on the sex ratio in the catches taken by European sealers prior to 1956 when walruses were completely protected in Northeast Greenland north of ≈ 72°N (Born et al. 1997)-effectively prohibiting the catch of walrus in East Greenland by foreigners. Hence for the period 1889-1955 an even sex ratio in the catch was assumed. South of 72°N the recent catch consists of approximately 90% males (Born et al. 1997), and we therefore applied a 0.9M:0.1F ratio for the catch in East Greenland after 1956. Post 2010, catches were set to the current quota with 90% of the catches being males.

Age-structure
An age-structure ( Fig. 2) is available from a sample of 376 walruses (0-29 years of age) that were caught in the North Water between 1987 and 1991. This catch shows a clear underrepresentation of animals younger than ten years of age, while the distribution for animals older than ten years approximates an exponential decline. Hence, we use the first half of the distribution to estimate the age-structured selectivity of the hunt, and the second half of the distribution to enhance our estimates of the population dynamic growth rate, as well as adult survival.
All the models of the Baffin Bay population were fitted to the age-structure sampled from the 1987 to 1991 harvest in the North Water. Some of the models for the West Greenland/Baffin Island and East Greenland populations were constructed without the use of the age-structured data, while other models for these populations were fitted to the age-structure sampled from the North Water hunt in order to include some data signal on the growth rate. From Figure 1, it should be noted that the period from 1960 to 1980, (which determines the age-structure of animals older than ten years of age in 1988) do not show similar harvest patterns for the three populations. This implies that there is increased uncertainty associated with the age-structure based growth rate estimates for the West Greenland/Baffin Island and East Greenland populations.

Population dynamics
All of our population dynamic models were age and sex-structured, but only some of them were fitted to the age-structured harvest data. Age-structure was maintained across all models in order to allow for a better comparison between models.
The changes in the exponential growth rate of the density-regulated and selection delayed models operated through changes in the birth rate. While the birth rate is likely affected, it is unlikely that it is the only life history parameter that is affected by density regulation and densitydependent selection. Alternatives include regulation of juvenile, or adult, survival as well as the age of maturity. Regulation of only the birth rate was nevertheless maintained, partly because we wanted to keep our models simple, but also because our data would not allow for a distinction between regulation via the different parameters. While a distinction between regulation of the birth rate versus regulation of the age of maturity would affect the posterior estimates of these parameters, it should generally not affect the age-structure of the population. The relative number of age-class zero individuals relative to the rest of the population, however, would be affected by a distinction between regulation of fecundity versus regulation of juvenile survival. But this component of the age-structure is also subject to harvest selection, which implies that our data are unable to distinguish regulation of fecundity from regulation of juvenile survival.
For a mathematical description of our model, let x be the maximum lumped age-class. Let the number N m/f a,t+1 of males (m) and females (f) in age-classes 0 < a < x in year t+1 be (1) (1) where p m/f a is the age specific survival rate of males/females, and c m/f a,t is the age specific catch of males/females in year t. The age and gender dependent survival rates p m/f a = p p m/f a are given as a product between a survival scalar p and a relative (0 < p m/f a ≤ 1) survival rate, with the relative survival rates being one for all age-classes, except age-class zero. The age and gender specific catches c m/f a,t = c m/f t c m/f a in year t is given as a product between the total catch of males/females (c m/f t ), as specified by the catch history, and an age-specific catch selectivity against the take of younger individuals (3) where a is age, s i ≠ 0 the selectivity parameter for data age-structure distribution i (convex selection when s i < 0, concave when s i > 0, and approximately linear when s i ≈ 0), and a s,i −1 is the last age-class with selection, i.e., c m/f a,i =1 for a > a s,i .
The number of females and males in age-class zero is N f 0,t = N 0,t and N m 0,t = (1− )N 0,t where is the fraction of females at birth, and (4) where a m is the age of the first reproductive event and B a,t the number of births from females in age class a, is (5) where b a,t is the birth rate in year t for age-class a females should they be at their age specific reproductive peak, ̃a is a relative age-specific birth rate (set to one for all mature age-classes), and M f a,t is the number of mature females in age-class a in year t, defined as (6) Let b a,t be for exponential growth for density-regulated growth for selection-delayed dynamics (7) where b is a constant birth rate, b* is the birth rate at population dynamic equilibrium (assuming zero catch and equilibrium denoted by *), b max is the maximal birth rate, b a,t is the average intrinsic birth rate for females in age-class a in year t (initially set to b* for all age-classes), γ is the density dependence parameter, with the abundance component that imposes density dependence being the one-plus component (8) For the selection-delayed models, the average intrinsic birth rate of offspring (age-class zero) born by age-class a females (9) is the average intrinsic birth rate of the parents b a,t plus a proportional response to density dependent natural selection, with the birth rate truncated here at b max and ι defining the response to selection. Hence, the average intrinsic birth rate of age-class zero is the weighted average across all offspring (10) Assuming that there is no change in the intrinsic fecundity rate of a cohort over time, b a+1,t+1 = b a,t and (11) Given a stable age-structure and no catch, let, for a traditional model of exponential or densityregulated growth, λ be a constant defined by N ︿ t+1 = N ︿ t . The sustainable yield is then sy = N ︿ (λ − 1), and for the density-regulated model there is an optimum ∂sy/∂N ︿ = 0; the maximum sustainable yield (msy) at N ︿ msy , also known as the maximum sustainable yield rate (msyr = msy/N ︿ msy ) and the maximum sustainable yield level (msyl = N ︿ msy /N ︿ * ). However, for selection-delayed dynamics the intrinsic growth rate is an initial condition, so that it is the acceleration of the growth rate that is determined by the density-dependent ecology, and not the growth rate itself, as is assumed in density regulated growth. For a given abundance, this implies that there is no constant λ to define a constant sustainable yield and, thus, there is no single abundance curve of sustainable yields, and no easily defined maximum of sustainable yield (Witting 2002).

Assessment models
The population dynamic analyses in this study are based on eight models. For the Baffin Bay population we use a model of exponential growth (eB) to provide an estimate of the growth rate, and we use models of density-regulated growth (dB) and selection-delayed dynamics (sB) to estimate the historical trajectory. Because the age-structured catch data are from the North Water, we only apply models with age-structured data for the Baffin Bay population.
For the West Greenland/Baffin Island population we use a density-regulated model with no agestructured data (nW), a density-regulated model with age data (dW), and a selection-delayed model with age data (sW). These models allow us to compare the estimates of historical trajectories and current growth given the different assumptions.
For the East Greenland population we apply only a density-regulated model with no age-structured data (nE), and a selection-delayed model with age-structured data (sE). A density-regulated model with age data was not applied, because initial runs showed that density-regulated growth was unable to reconcile the historical catches with the current abundance and the growth rate suggested by the age-structured data. This inability may not necessarily be because the age data from the North Water are inapplicable to the East Greenland case, it may instead reflect a common inconsistency between density-regulated growth and the true historical trajectories of exploited populations of marine mammals (Witting 2003(Witting , 2013.

Statistical methods
We applied a Bayesian statistical analysis as developed for baleen whales in the Scientific Committee of the International Whaling Commission (IWC) (Punt and Butterworth 1999;Wade 2002;Witting 2003). The assessment models were fitted to data by projecting the population under the influence of the historical catches, with the initial abundance of the density-regulated and selection delayed models reflecting a pre-harvested population in dynamic equilibrium. A Bayesian statistical method (e.g., Berger 1985;Press 1989) was used, and posterior estimates of model parameters and other management-related outputs were calculated. This implied an integration of the product between a prior distribution for each parameter and a likelihood function that links the probability of the data to the different parameterisations of the model.

Priors
All priors were uniform on ordinary or logarithmic scales. The range of values for the priors, together with the value for the fraction of females at birth, were based on various studies of walrus populations and previous modelling of walrus population dynamics (Fay 1982;DeMaster 1984;Fay et al. 1989Fay et al. , 1997Chivers 1999). The annual natural survival rate of adults (p) was set in this study to range between 0.97 and 0.99, with the relative survival of all adult age classes being the same. The survival rate for age-class zero individuals (p 0 ) covered the range from 0.7 to 0.95.
The birth rate (b max , or b for the exponential model) was set by a uniform prior from 0.35 to 0.65 (Mansfield 1958;Fay 1982;Born 2001), reflecting a point estimate where each female produces an offspring every second year. The age of reproductive maturity (a m , first reproductive event) was set from five to nine years of age, and a uniform sex ratio of offspring was assumed.
The strength of density regulation (γ) in the density-regulated model was given by a uniform prior from 1.5 to 5, in order to approximate a maximum sustainable yield level (msyl) from 0.50 to 0.70. The prior on the density dependence (γ) and the selection response (ι) of the selectiondelayed model were set by trial and error. For reasons of simulation efficiency, the two parameters were given identical and fully correlated priors for the Baffin Bay and West Greenland/Baffin Island population (i.e., reducing the two parameters to one), while they were kept as independent parameters for East Greenland walruses because initial simulations indicated that this was required in order to get the model to fit the data. A uniform prior from zero to one was set on the catch history parameter c h , to reflect a linear transition from the low to the high catch history.
The priors on the initial (N 0 ), and the initial equilibrium (N * ), abundance were set by trial and error to be wider than the posteriors; so was the prior of the age-structured catch selectivity (s d ), in order to provide best estimates assuming that the last age class with selection was nine years.

Bayesian integration
The method of de la Mare (1986) was used to calculate the log likelihood for the abundance component under the assumption that observation errors are log-normally distributed (Buckland 1992) Walrus of the North Atlantic 200 (12) where N ︿ i,t is the point estimate of the ith set of abundance data in year t, cv i,t is the coefficient of variation of the estimate, and N t is the simulated abundance. For models with age-structured data, the log likelihood (ln L = ln L n + ln La) was given as a sum between the abundance (ln L n ) and an age-structured (ln L a ) component. Following Punt (2006), the age-structure was assumed to be multinomially distributed, with (13) where i denote the ith set of age-structured data (which were sampled in year t), n i is the number of individuals in the sample, ĉ i,t,a the observed fraction of individuals in age class a in year t, and c t,a the fraction simulated by the model.
The Bayesian integration was obtained by the sampling-importance-resampling routine (Jeffreys 1961;Berger 1985;Rubin 1988), where n s random parameterisations θ i (1≤ i ≤ n i ) are sampled from an importance function h(θ). This function is a probability distribution function from which a large number, n s , of independent and identically distributed draws of θ can be taken. h(θ) shall generally be as close as possible to the posterior, however, the tails of h(θ) must be no thinner (less dense) than the tails of the posterior (Oh and Berger 1992). For each drawn parameter set θ i the population was projected from the first year with a harvest estimate to the present. For each draw an importance weight, or ratio, was then calculated (14) where L(θi) is the likelihood given the data, and h(θi) and p(θi) are the importance and prior functions evaluated at θi. In the present study the importance function is set to the joint prior, so that the importance weight is given simply by the likelihood. The n s parameter sets were then re-sampled n s times with replacement, with the sampling probability of the ith parameter set being (15) This generates a random sample of the posterior distribution of size n r . Details regarding the sampling statistics are given in Table 2.
All models except the selection-delayed model for East Greenland (sE), show well sampled posterior distributions. The maximum importance weight of a parameter set relative to the sum of importance weights for all the sampled sets was between 0.001% and 4.9% across all models except sE. These models also have a proportion of unique parameter sets in the resample between 20.8% and 99.5%, and a maximum occurrence of a unique parameter set in the resample between 0.2% and 4.7%. The resample of the sE model have only 7.4% unique parameters, even though it had the largest number of initial draws from the prior, with 6.1 million initial parameterisations. The primary reason for the few unique parameter sets in the sE resample is that this is the only selection-delayed model that attempted to provide independent estimates of the two parameters ι and γ.

Posterior distributions
The realised prior and posterior distributions are shown in Figures 3 to 6 for selected models, and the posterior parameter estimates are listed in Table 3. The abundance parameters N 0 and N * are strongly updated in all cases. This may not be immediately clear from the figures; they often show a posterior that fade out at the upper range of the prior, but maintains its weight at the lower bound. The reason for this is that all parameterizations that resulted in zero abundance during their projection were removed from the realised prior. This implies that the density of the realised prior may approach zero at the lower limit in some models. The true priors, from which N 0 and N * were drawn, would have maintained their density at the lower range, showing a clearly updated posterior.
The problem with parameterizations that go extinct at the lower limit of the posterior abundance distribution arises primarily because of the lack of precision in the abundance estimates relative to the absolute estimate and the catch histories. This was a particular problem for an initial model of the West Greenland/Baffin Island population that did not include the more precise abundance estimate from Baffin Island. The median projection of this model was positively biased relative to the abundance estimates. This problem was generally avoided by including the more precise estimate, although the median of the density-regulated model with no age data remained slightly positively biased (Fig. 7).
Other general features across models include that there is no clear updating of the catch history parameter c h , which shows that our models cannot generally inform us regarding the absolute magnitude of the catches within the range considered. However, all models with age-structured data show a very clear updating of the age-structured selectivity (s b ), with all models agreeing that s b > 0. Figure 3 gives the prior and posterior distributions for the exponential model for Baffin Bay. These distributions show a strong updating of the growth rate (r) parameter, an updating that is provided primarily by the age-structured data from animals older than ten years of age. The figure also illustrates a weak, or almost no, updating of the life history parameters themselves (survival and birth related parameters). This shows not only that we have no information to determine which constellation of life history parameters that determine the growth rate within the limits considered, but it also illustrates that our priors on the life history parameters are not severely biased relative to the growth rate of the age-structured data. A biased life history prior would have an associated posterior distribution with most weight at the upper, or lower, end of the prior.
Walrus of the North Atlantic 202 Table 2. Sampling statistics for the different models (M). The number of parameter sets in the sample (nS) and the resample (nR), the maximum importance weight of a draw relative to the total importance weight of all draws, the number of unique parameter sets in the resample, and the maximum number of occurrences of a unique parameter set in the resample. nS and nR are given in thousands. eB, dB and sB are the exponential, density-regulated, and selection-delayed models for Baffin Bay walruses. nW, dW and sW are the density-regulated (no age data), density regulated, and selection-delayed models for West Greenland/Baffin Island walruses. nE and sE are the density-regulated (no age data) and selectiondelayed models for East Greenland walruses. The contrast to an updated growth rate is given by the density-regulated model for East Greenland (Fig. 4). This model is fitted only to a single abundance estimate, which implies that the data contains no information on the underlying growth rate of the population. As a result we see a posterior growth rate distribution that is practically identical with the prior. This model is extremely uninformative because it updates only the absolute abundance estimate.
When provided with age-structured data, the density-regulated models for Baffin Bay (Fig. 5) and West Greenland/Baffin Island show a clear updating of the growth rate. However, these models will generally not update the density regulation parameter (γ) that determines the strength and shape of the density dependence curve. The density dependence (γ) and selection response (ι) parameters are instead updated for the selection delayed model (Fig. 6 for Baffin Bay). This is largely because the current growth rate in the selection-delayed model to a large degree is determined by the density-dependent acceleration of the growth rate (which is defined by ι).

Age-structured selectivity
The fit of the exponential Baffin Bay model to the age-structured data from the North Water is shown in Figure 2. There is an excellent fit to the age-structure and the fit of the other models (not shown) are basically the same. The assumption that there is selection only against animals younger than ten years of age provide similar overall shapes between the data and the predicted age-structures. This assumption is in agreement with a hunt that prefers adult animals.  Table 3.
The estimated selectivity s d is practically the same in all cases (0.18 [95% CI: 0.11-0.24] for model eB; point estimates range from 0.12 to 0.18 across all models). It shows a steep concave selection profile, which is characteristic of a strong selection for full-grown adult animals, where there is a relatively strong selection even against the take of animals that are almost but not yet fully grown. Had the estimated selectivity function been convex instead, there would have been a softer selection against younger animals, where almost but not yet fully grown animals were preferred almost at the same level as full-grown animal.

Population dynamics
The estimated population trajectories of the different models are shown in Figure 7.

Baffin Bay
The exponential model estimates that the Baffin Bay population would increase with an annual growth rate of 7.7% (95% CI:6.  Both the density-regulated and the selection-delayed model suggest a relatively large decline in the number of walruses in this population during the 20 th Century. The density regulated model estimates a 2012 depletion ratio of 0.27 (95% CI:0.18-0.39), while the selection-delayed model suggests an even larger decline to 0.16 (95% CI:0.10-0.39). However, there is uncertainty in these estimates because they depend on the constructed catch history.

West Greenland/Baffin Island
The two density-regulated models and the selection-delayed model all agree that the population of West  Table 3.
catches peaked with more than 600 animals taken in 1938 and more than 9,400 animals taken in total over the three decades. The density-regulated models, with and without age data, estimate a maximal depletion ratio of 0.16 and 0.23 around 1960, and the selection-delayed model estimates a maximal depletion of 0.10.
Again, all models agree, that the population increased from the early 1960s to 1993 when catches were relatively low, with an average 1993 estimate of 3,100 (95% CI: 2,500-4,400) walruses. Then, owing to increased catches there was a minor decline until the early 2000s. Annual catches where then cut from an average of 190 per year from 1995 to 2005 (low catch history), to a current quota of 61 animals, and the population is again increasing with a 2012 estimate of 3,900 (95% CI: 2,500-5,300) individuals.
Both of the density-regulated models estimate a smaller current growth rate than the exponential model for the Baffin Bay population. With no age data the density-regulated model estimates 5.8% (95% CI:1.0-8.3%) growth, and with age data this increases 7.3% (95% CI:0.9-9.2%).   Fig. 7. The estimated population trajectories shown by the median and 95% credibility intervals of the different models, together with the abundance estimates and their 95% confidence intervals. eB, dB and sB are the exponential, density-regulated, and selection delayed models for Baffin Bay walruses. nW, dW and sW are the density-regulated (no age data), density-regulated, and selection-delayed models for West Greenland/Baffin Island walruses. nE and sE are the density-regulated (no age data) and selection-delayed models for East Greenland walruses.

NAMMCO
The selection-delayed model estimates a slightly higher growth rate of 9.3% (95% CI:8.3-10%). The increased estimate of the selection-delayed model is most likely caused by the positive trend over the tree abundance estimates, while the lower estimate for the density-regulated model is caused by its internal constraints that limit the range of possible growth rate estimates. With 2012 abundance estimates between 3,600 (95% CI:2,300-5,500) and 4,100 (95% CI:2,800-5,500) animals, the three models estimate a 2012 replacement yield between 210 (95% CI:50-310) and 390 (95% CI:220-610) walruses.

East Greenland
The two models for East Greenland walruses estimate quite different trajectories. The historical catch in East Greenland was characterised by relatively large takes in the late 19 th and early 20 th Century; the low catch history estimates that 268 walruses were taken in 1889, and that 573 walruses were taken between 1898 and 1908. Since that time takes have been low, except the period from 1925 to 1939 where a total of approximately 350 walruses were taken. Based on the abundance estimate from 2009, the density-regulated model estimates that these takes caused only small perturbations, with a maximal depletion to 80% in 1890. The population has since been increasing slowly and is now basically recovered with a current growth rate of 1.5% (95% CI:0.5-3.9%) in the absence of catches.
Owing to its inability to allow for population growth in populations at carrying capacity, the density-regulated model is unable to reconcile the growth rate of the age-structured data from Baffin Bay with the abundance and catch data from East Greenland. The three data components may instead be reconciled by the selection-delayed model. It estimates a large decline in the abundance from an initial equilibrium abundance of 1,200 (95% CI: 810-1,300) individuals in 1888 to only 18 individuals in 1957. Since then the population has increased steadily to a fully recovered population (depletion ratio of 1.1 [95% CI:0.62-2.7]) with 1,200 (95% CI:700-2,800) individuals in 2012. With an estimated annual growth rate of 3.3% (95% CI:1.6-4.3%), we obtain a replacement yield of 41 (95% CI:22-100) individuals in 2012.

Biological parameters
The ranges we have chosen for the various prior distributions were based on previous knowledge of the biology of Atlantic and Pacific walruses and estimates of their vital parameters. Although recognized as two geographically and taxonomically distinct subspecies (e.g. Fay 1985), the life history of Atlantic and Pacific walrus appears to be very similar. However, in both cases the segregation of different sex and age classes for most of the year, and the selective hunting pattern make it difficult to obtain unbiased samples for determining biological parameters (Fedoseev and Goltsev 1969;Fay 1982;DeMaster 1984). Our prior and posterior distributions of the biological parameters are generally in good agreement with estimates in other studies.
The sex ratio in walrus populations is not well known. On the assumption that walruses are polygynous, an adult sex ratio of 1 male to 3 females has been suggested in the population of Pacific walrus (Fay 1982;Fay et al. 1984;Sease and Chapman 1988). The sex ratio in the Pacific walrus changed from (M:F) 1:1.7 to 1:3 between 1960 and 1985 during a period of population increase according to Fay et al. (1997), a pattern that might actually be expected from density and frequency dependent natural selection (Witting 1997(Witting , 2000b. To our knowledge there are no indications of the sex ratio being severely skewed towards females in Atlantic walrus subpopulations, and we therefore assumed that it was even. This is in accordance with findings in Atlantic walrus from the North Water (Born 2001) and in Pacific walrus from the Bering Strait region (Fay 1982;Fay et al. 1984).
The natural mortality rate in walruses is not well established but is assumed to be low because productivity is low and longevity is relatively high (Fay 1982;Fay et al. 1997). Natural mortality rate has been estimated to be 3 to 5% for the entire population of the Pacific walrus (DeMaster 1984;Fay et al. 1989). Fay et al. (1994) suggested that natural mortality in adults was higher than 1% but probably not higher than 2% per year. A natural mortality rate of 1.5% per year was applied in simulations of Bering Sea walruses (Fay et al. 1997). Chivers (1999) used an estimate of survival in adults between 9 and 24 years of 0.98 (it then decreased until 40 years of age). We applied a uniform prior from 0.97 to 0.99 for adult mortality, and obtained a posterior estimate of 0.98 (95% CI:0.97-0.99).
Survival in calves until their second year of life (when suckling ceases) has been estimated to range between 0.5 and 0.9 (Fay et al. 1997). We applied a corresponding uniform prior between 0.7 and 0.95 for the annual survival rate for walruses in age class zero, and obtained a posterior estimates of 0.83 (95% CI:0.71-0.95) for the exponential model of Baffin Bay walruses, and estimates from 0.80 (95% CI:0.70-0.94) to 0.92 (95% CI:0.73-0.95) across all models.
Based on analyses of reproductive organs, rates of fecundity have been estimated at between 0.29 and 0.40 in walruses (Mansfield 1958;Fay 1982;Garlich-Miller and Stewart 1999;Born 2001). According to Mansfield (1958), the reproductive cycle of the female Atlantic walrus in Foxe Basin was basically biennial, but, to an unknown extent, older females may give birth at three or four year intervals. According to Fay (1982) female Pacific walruses tend to breed at 2year intervals or less often, and hence maximum fecundity was thought to be 0.5 (Fay et al. 1997). Our prior of maximum annual fecundity ranged between 0.35 and 0.65 calf/mature female/year, in order to set a statistically neutral prior around a maximal fecundity of 0.5. We obtained a posterior estimate of 0.52 (95% CI:0.38-0.64) for the exponential model of Baffin Bay walruses, and estimates ranging from0.45 (95% CI:0.35-0.64) to 0.65 (95% CI:0.54-0.65) across all models. Mansfield (1958) found that the age at first ovulation varied from 5 to 10 years, but that the majority of females in his Canadian sample became sexually mature at the age of 7. Born (2001) found that the average age at sexual maturity was 6 years in Atlantic walruses from the North Water and stated that attainment of sexual maturity in female Atlantic walruses is similar to that in the Pacific subspecies. By the age of 6 years, two thirds of the female Pacific walruses have ovulated at least once and by the 8th or 9th year practically all have ovulated (Fay 1982). Our prior for the age at maturity ranged between 5 and 9 years of age, and our posterior estimates ranged from 5.1 (95% CI:5.1-6.5) to 7.3 (95% CI:5.1-8.9).
The population birth rate (fraction of newborn in the total population) has been estimated to be 0.07 (Mansfield 1966) or 0.11 (Mansfield 1973) in Atlantic walrus and between 0.12 and 0.17 in Pacific walrus (Fedoseev and Goltsev 1969;Fay 1982). Instantaneous net growth rate of the population of Pacific walruses during the late 1950s to mid-1970s was estimated to be 0.067 (Tavrovski 1971;Sease and Chapman 1988), indicating a finite growth rate of about 7% per year for a population in a phase of growth under favourable environmental conditions with no food limitations. Chivers (1999) modelled an annual maximum growth rate of 8% but stressed that because survival rates are unknown for walruses, the growth rate of the model should not be considered an estimate of maximum growth rate for walruses.
We did not apply a prior to the growth rate, but a realized prior was generated from the priors on other parameters. With only a few, imprecise abundance estimates, our estimate of the population dynamic growth rate is determined by the age-structured data. This implies that our growth rate estimate will be positively correlated with our prior on the natural survival rate of adults (p); more exactly, we expect that if the median of our survival prior is 0.01 too high relative to the true natural survival rate, then, our estimate of the exponential growth rate would also be approximately 0.01 too high (Witting 2012). Hence, we sat the median of the survival prior to 0.98; the point estimate of adult survival obtained in most studies (see discussion above). And with a posterior survival estimate of 0.98 (95% CI:0.97-0.99) for Baffin Bay, we are confident that our 2012 growth rate estimate of 7.7% (95% CI:6.7-8.9%) for this population (assuming no harvest) is unlikely to be positively biased.

Population dynamics
The inclusion of the age-structured sample from the North Water hunt allowed us to estimate the exponential growth rate of the Baffin Bay population. The age-structure was also applied to the two other populations because it represents the only data that we have on the population dynamic growth rate of walruses in Greenland. These populations, however, may show other growth rates.
The uncertainty on the historical catch from the Baffin Bay population makes estimates of the historical trajectory uncertain. The trajectory of the exponential model from 1960 to the present ( Fig. 7; model eB), and the slow increase for the near future, may represent what we can estimate with confidence for this population.
The historical catches are much more certain for the West Greenland/Baffin Island population, where the largest uncertainty in relation to estimates of the historical trajectory rest on the assumed models and our lack of a population specific growth rate estimate. It is thus reassuring that the three models that differ in dynamics and growth rate assumptions all estimate a similar shaped trajectory with an initial decline from 1900 to the 1960s, followed by an increasing phase until 1993, and a second decline until the drop in the annual catches around 2005 caused the population to increase again. The major differences in the trajectories of the three models lie in the scaling to absolute abundance over time. Nevertheless the three estimates of current depletion are not that different ranging from 0.40 (95% CI:0.23-0.74) to 0.67 (95% CI:0.39-0.98).
East Greenland is like West Greenland/Baffin Island in the sense that the largest uncertainty lies in the applied model and the lack of a population specific growth rate estimate. But, unlike the West Greenland/Baffin Island case, the trajectory for East Greenland walruses depends quite strongly on the assumed dynamics. If the dynamics are strictly density-regulated with no selection response, we expect a flat trajectory with a recovered population and almost no current growth. Selection-delayed dynamics, on the other hand, estimate a recovered population that is currently increasing and was historically strongly depleted to around 2% in 1957.

Comparison to 2005 assessment
Apart from the dynamics of the selection-delayed model, our results for walruses in East Greenland are similar to the results in 2005 (Witting and Born 2005). The 2005 assessment suggested a carrying capacity of 1,600 walruses, and today we estimate an equilibrium abundance in the absence of harvest that ranges from 1,200 (95% CI:810-1,300) to 1,500 (95% CI:740-3,100) individuals.
The populations of Baffin Bay and West Greenland/Baffin Island walruses are estimated to be in much better shape than indicated in 2005. Given the catch history and the abundance estimate available in 2005, we projected extinction of West Greenland walruses around year 2000, and the population in the Baffin Bay was estimated to be depleted to 2% in 2010 given continued catches at the 1999 level. Averaging across the models in this paper, we now estimate a 2012 depletion of 0.52 (95% CI:0.30-0.9) for the West Greenland/Baffin Island population and 0.22 (95% CI:0.14-0.39) for the population in Baffin Bay. Furthermore, it is estimated that both populations have increased over the last half decade, and the increase should be maintained for the near future given current catch levels.
A main reason for the improved status assessment for West Greenland/Baffin Island walruses is a change in the abundance estimate. In 2005, we had only one old estimate of approximately 1,000 walruses in 1990, while we now have three independent estimates from 2007 to 2009 of approximately 3,000 individuals each. For Baffin Bay, however, we have basically the same abundance information available (a rough 1999 estimate of 1,500 individuals in 2005, versus an average of 1,500 individuals for the 2009 and 2010 estimates in the current assessment). But by fitting the assessment model to the age-structured data from the North Water, we obtained a much more solid estimate of the growth rate (the posterior of the msyr in 2005 was 2% [90% CI:0-7%], versus 7.7% [95% CI:6.5-8.8%] today) and thus we obtain a much more positive estimate of current status.