REPORT OF THE NAMMCO-ICES WORKSHOP ON SEAL MODELLING

To support sustainable management of apex predator populations, it is important to estimate population size and understand the drivers of population trends to anticipate the consequences of human decisions. Robust population models are needed, which must be based on realistic biological principles and validated with the best available data. A team of international experts reviewed agestructured models of North Atlantic pinniped populations, including Grey seal (Halichoerus grypus), Harp seal (Pagophilus groenlandicus), and Hooded seal (Cystophora cristata). Statistical methods used to fit such models to data were compared and contrasted. Differences in biological assumptions and model equations were driven by the data available from separate studies, including observation methodology and pre-processing. Counts of pups during the breeding season were used in all models, with additional counts of adults and juveniles available in some. The regularity and frequency of data collection, including survey counts and vital rate estimates, varied. Important differences between the models concerned the nature and causes of variation in vital rates (age-dependent survival and fecundity). Parameterisation of age at maturity was detailed and time-dependent in some models and simplified in others. Methods for estimation of model parameters were reviewed and compared. They included Bayesian and maximum likelihood (ML) approaches, implemented via bespoke coding in C, C++, TMB or JAGS. Comparative model runs suggested that as expected, ML-based implementations were rapid and computationally efficient, while Bayesian approaches, which used MCMC or sequential importance sampling, required longer for inference. For grey seal populations in the Netherlands, where preliminary ML-based TMB results were compared with the outputs of a Bayesian JAGS implementation, some differences in parameter estimates were apparent. For these seal populations, further investigations are recommended to explore differences that might result from the modelling framework and model-fitting methodology, and their importance for inference and management advice. The group recommended building on the success of this workshop via continued collaboration with ICES and NAMMCO assessment groups, as well as other experts in the marine mammal modelling community. Specifically, for Northeast Atlantic harp and hooded seal populations, the workshop represents the initial step towards a full ICES benchmark process aimed at revising and evaluating new assessment models.


INTRODUCTION
The NAMMCO-ICES Joint Workshop on Seal Modelling (WKSEALS 2020) was convened on the basis of a recommendation from the 2019 meeting of the ICES-NAFO-NAMMCO Working Group on Harp and Hooded seals (WGHARP). WGHARP focuses on harp seal stocks in the Northeast (NE) Atlantic (i.e., the Barents Sea/White Sea and the Greenland Sea populations) and the Northwest (NW) Atlantic (i.e., populations in Canada and Greenland), as well as hooded seal populations breeding in Greenland and Canada. The working group is tasked with compiling and analysing data regarding these populations and where sufficient data are available, aims to perform stock assessments that could be used as basis for scientific advice, including catch options. During the 2019 meeting, WGHARP raised concerns regarding the existing population assessment models and their projections and recommended that "ICES and/or NAMMCO convene a workshop on population assessment models for seals in the North Atlantic to advance model development in the ways identified as required, before the next WGHARP" (ICES, 2019).
The limitations identified in the models included a poor fit to historical time-series harp seal pup count data in the Northeast Atlantic/Greenland populations, especially where substantial changes in pup counts occurred over a short time. The models' future projections, suggesting rapid population recovery, were thought to be questionable.
Based on the recommendation from WGHARP, the NAMMCO-ICES Joint Workshop on Seal Modelling (WKSEALS 2020) was therefore convened to: a) Compare approaches used for modelling populations of a range of northern pinniped species with similar life histories, focusing on species for which abundance estimates are primarily based on pup counts (see Table 1); b) Evaluate the appropriateness of data and methods to determine seal stock status and investigate methods for short term predictions, with the evaluation to include consideration of: • Data input required • Data sampling frequency • Comparison of modelling approaches, including deterministic and stochastic models • Comparison of model-fitting approaches • Model validation • Model convergence time • Further inclusion of environmental drivers, multi-species information, and ecosystem impacts for stock dynamics in the assessments and outlook, to the extent possible.
c) Develop recommendations for future improvements of the assessment methodology and data collection.
The outputs from this workshop will feed into an ICES benchmark for harp and hooded seals planned for 2022. ICES benchmarks are regularly organised with the aim to agree on independently reviewed assessment methodologies that are to be used in future update assessments. The result will be the 'best available' method that ICES advice can be based on.
Due to travel restrictions implemented to limit the spread of COVID-19, the workshop had to be transitioned to an online format. To accommodate working online across various time zones, three 2-hour plenary sessions were held over 4-6 November 2020. In between these plenary sessions, four 3-hour sub-group meetings were held, with times scheduled to allow participants on both sides of the Atlantic to join at least 2 of the sessions. A copy of the agenda and participant list for the workshop are available as Supplementary Files 1 & 2. To help the participants exchange information, code, and data related to their models, a GitHub site was established (https://github. com/ices-eg/wk_WKSEALS-2020/), where this information remains freely available, also for non-participants in the workshop.

CURRENT SEAL POPULATION MODELS
This workshop focused on models that describe the characteristics of: • harp seals in the Northwest Atlantic, the Barents Sea/White Sea and the Greenland Sea; • hooded seals in the Northwest Atlantic and Greenland Sea; and • grey seals in the UK, Canada/US, and the Wadden Sea.
These species share similar life histories and estimates of abundance are based on age-structured deterministic models fit to pup counts. The raw pup count data may require processing and interpretation because not all the seals born in a given year may be seen during a single survey count, and not all sites may be surveyed. These pre-steps may influence how the population models are designed to fit the processed data. The models also differ in how demographic parameters such as survival and fecundity are specified in the model structure, as well the type of data or processes that are included to describe variation in these vital rates. Furthermore, there are different approaches to fitting the models, which influence speed, inference, the use of informative priors, and model selection and validation. A summary overview of each of the seal population models compared at the workshop has been provided by developers and/or current users of these models and is presented below.

Summary by T. A. Øigård and M. Biuw
Assessments of the population status and hunting potential for harp and hooded seals in the Greenland Sea and harp seals in the Barents Sea/White Sea are carried out by WGHARP. The set of iterative catch equations for calculating the number of seals of age 1 year and older is given by: The number of 1 year old seals for year y is given by the number of pups the previous year minus the pup catch the previous year and a survival rate. This iterates through all age groups until a maximum age group A where seals of age A and older are pooled together. Data on age structured catch is not available. Catch data have historically been summarised annually by the number of pups and 1+ seals separately (see Figure 1) and is thus not fully age-specific. To obtain age structured catch, it is assumed that the age distribution in the catch is the same as the age proportion of the various age groups. Reproduction (number of seal pups born in year y) is given by:

=1
Here p is the age-specific maturity and F is the proportion of mature females that are pregnant.
The model is fitted to the following input data: • Number of pups from breeding-season surveys (~5-year interval) Priors on estimated parameters: the initial population size K, the pup mortality M0 and the adult mortality M1+ The model attempts to fit past population trajectory to survey estimates of pup production and can also be used to project the future population trajectory (for information, Figure 2 shows projections for the next 15 years). While 15 years is substantially longer than the 1-3 years traditionally used in fisheries assessments, this is justified by: 1) the longevity of seals, 2) the high mean age of first reproduction (about 4-6 years) and the fact that pup production surveys are only carried out every 5 years (minimum period). As part of the planned ICES benchmark meeting, the 15-year projection horizon that is currently used will be re-evaluated in light of model developments, and the harvest control rules that are based on these projections will also be reassessed. For the Barents Sea/White Sea harp seal population model, there appears to be some lack of fit ( Figure 2). There are a number of difficulties in parameterising this model. The problems are most clearly illustrated for the Barents Sea/White Sea harp seal stock, while the issues have been less problematic for the other two stocks managed using the same model, namely Greenland Sea harp and hooded seals. Most importantly, the model is parameterised with only three parameters (K, M0, and M1+). The mortalities are assumed to be constant for all years, and hence it is not possible for this model to capture rapid changes, such as those shown in the pup production estimates for the Barents Sea/White Sea stock. The main argument for not including more complexity and flexibility in the model has been the sparsity of available data.
While catch statistics are available since the mid-1940s, pup production estimates are only available since the 1990s for the Barents Sea/White Sea population, and all three stocks managed using this model are surveyed only every 5 years or so. Also, fecundity and maturity data are available infrequently. In addition, only early pregnancy data are available, so the model is not accounting for possible late-term abortions, which may reduce the natality rate. Natality (number of pups born) is probably overestimated by the current model, and consequently, mortality is probably overestimated as the model-fitting attempts to compensate in order to capture the overall population trend. Several projects have been initiated to improve the model. A stochastic model was developed by Øigård and Skaug (2015). This model attempted to model variation in fecundity as an autoregressive random process and was able, to a much higher degree, to capture the dynamics seen in the pup production data. However, the high uncertainty estimates in future projections produced by this alternative model resulted in WGHARP not accepting this as an appropriate assessment model at that stage. Recently, modifications of this model to include environmental and biological covariates as drivers of change in vital rates have been explored, with promising results. WGHARP agreed this could be a very useful approach (ICES, 2019). Work along these lines will continue to be pursued and was given support by the WKSEALS workshop.

Summary by M. Hammill, G. Stenson and A. Mosnier
Northwest Atlantic harp seal population abundance is estimated by fitting a model to independent estimates of the total pup production, and reproductive rates observed for seals 8 years of age and older (referred to as 8+) (Hammill, Stenson, Doniol-Valcroze, & Mosnier, 2015). The model assumes that density-dependent factors are affecting demographic parameters that can be described using a theta-logistic curve acting on reproductive rates and juvenile survival. The value of theta is 2.4. It is also assumed that the sex ratio is 1:1. The input data consist of periodic estimates of pup production from aerial surveys and mark-recapture studies, annual removals (reported catch, struck and loss, and by-catch), and annual reproductive rate data (1952-present) (Hammill, Stenson, Mosnier, & Doniol-Valcroze, 2021;Stenson, Buren, & Koen-Alonso, 2016;Stenson, Buren, & Sheppard, 2020;Stenson & Upward, 2020).
In an earlier version, the model estimated adult mortality rate and fixed juvenile mortality as three times the adult rate (Hammill et al., 2015). However, in most large mammal populations adult mortality rates tend to be low and show little inter-annual variation, while juvenile mortality rates are thought to be more variable. In a newer version of the model, adult mortality is fixed at 0.03. The model is fitted by adjusting the starting population size, first year (juvenile) mortality rates and a model estimate of carrying capacity (K) to minimise the sum of square differences between model estimates of pup production and reproductive rates and observed values for these parameters. As the model has evolved, the least squares fitting process has been given a somewhat Bayesian flavour as certain limits have been defined for the fitted parameters, for example, values for K are restricted to values between 7 and 15 million (Hammill, Stenson, & Kingsley, 2011).
Environmental variables are also incorporated into the model including a measure of ice related mortality that acts on pup survival prior to animals entering the water to forage independently and a comprehensive environmental index (CEI) that is multiplied by K (Hammill et al., 2021;. This is done because it is expected that K is not constant.
Instead, there will be fluctuations in K in response to interannual changes in environmental conditions, which will affect productivity and possibly young of the year survival. The CEI allows for some of this environmental uncertainty to be incorporated into the model.
The model also has a projection component to provide advice on future population trends. The uncertainty in reproductive rates is projected forward by assuming that changes in reproductive rates are driven by density dependent relationships or the model draws from a sample of reproductive rates that have been observed over the last 5 or 10 years. There is a provision to incorporate into the model environmental factors that will also affect productivity, but this component has not been switched on due to uncertainty into how to forecast foraging conditions beyond the next couple of years. Ice related mortality varies in the model by drawing from a sample of ice mortality rates that have been observed over the last decade.
In the most recent assessment, ice-cover projections from a climate change model for the Newfoundland-Labrador coast were included to provide insights into the long-term outlook for the stock (Han, Colbourne, Pepin, & Xie, 2015;Han, Ma, Long, Perrie, & Chasse, 2019). Future plans are to modify the model into a true Bayesian model implemented using STAN.

Summary by G. Stenson and M. Hammill
Northwest Atlantic hooded seals pup in three areas: the southern Gulf of St. Lawrence ('Gulf'), off southern Labrador/northeast Newfoundland ('Front') and in the Davis Strait. The Front is the largest area for pupping, accounting for ~90% of pupping in 2005. The most recent assessment of the population occurred in 2006 (Hammill & Stenson, 2006). The estimate of total abundance and population trends since 1960 was obtained from an age structured model (programmed in Excel) based upon the existing harp seal population model used at the time. Estimation is carried out using least squares model fitting (see Northwest Atlantic harp seal model).
To date, the model has incorporated reported catches from Canada (1952Canada ( -2005, Greenland (1955Greenland ( -2002 and the Norwegian moulting hunts , as well as age specific pregnancy rates (constant across years) based upon back calculations of pregnancy rates from females collected in the whelping patches between 1979 and 2003. Reported catches have been adjusted to account for seals killed but not landed or recorded using the same estimates applied to harp seals. Increased mortality of young of the year hooded seals (m=0.25) has been assumed to occur in 1981, a year with extremely poor ice and in 2005 when a mid-March storm resulted in the documented death of a large number of this age class.
Four estimates of pup production at the Front are currently available (1984, 1990, 2004 and 2005), but unfortunately, the only year in which all three pupping areas were surveyed was 2005. Therefore, the model has been fit to the estimates of pup production at the Front alone, as well as to the total population by making assumptions about pup production in areas where there were missing estimates.
Since the last assessment, little new data have become available. Coltman et al. (2007) found that Northwest Atlantic and Greenland Sea hooded seals are not distinguishable genetically although WGHARP has continued to treat them separately for management purposes. Frie, Stenson and Haug (2012) re-examined the reproductive data and found that the mean age of pregnancy increased by almost 2 years between 1956 and 1995. There also appeared to be a decline in pregnancy rates since the 1990s, which they attributed to ecosystem changes rather than population mediated density dependence.
Hooded seal pups at the Front have also been counted on images obtained during harp seal pup production surveys in 2012 and 2017 and are currently being analysed to determine if it is possible to obtain newer estimates of pup production for this area. Unfortunately, however, there are only a few new samples that can be used to estimate recent reproductive rates.
Plans are underway to update the population model for Northwest Atlantic hooded seals to a format that incorporates additional ecological indicators, similar to that being used for harp seals.

Summary by L. Thomas, C. Fagard-Jenkin and F. Empacher
UK grey seal population size is estimated using a Bayesian statespace model that is described in detail by Thomas et al. (2019), with latest results summarised by Thomas (2020). Input data are: 1) estimates of pup production, derived from aerial surveys of breeding colonies and aggregated into four regions, for the years 1984-2010, and biennially thereafter; and 2) estimates of total population size in 2008 and 2014, from aerial surveys of haul-out sites with proportion hauled out calculated from telemetry data. These data are fitted to a population dynamics model that tracks seven age classes of seal (age 0, 1, 2, …, 6+) in each of the four regions. The model optionally includes density dependence in pup survival or fecundity, and allows for densitydependent movement of recruiting females between regions.
Only the density dependent survival model is currently used, and the form of the density dependent function is: Where Φ , , is pup survival in region at time , Φ is maximum pup survival, 0, , −1 is pup production in region at time − 1 and and are model parameters related to carrying capacity and the shape of the density dependence function.
The population model incorporates demographic stochasticity but not environmental stochasticity and does not currently have any covariates affecting the demographic parameters. The population model does not include the adult male component of the population, which is included by a sex ratio parameter. Informative prior distributions are used on the demographic parameters and are essential for adequate model fit (i.e., the input data sources alone are not sufficient to accurately estimate model parameters).
Model fitting is via a particle filtering (also called sequential importance sampling) algorithm implemented in customwritten C code. This currently takes 1-2 days for model running and post-processing. Ongoing research aims to speed up model fitting in three ways: 1) revising the algorithm using a mixed particle filtering/MCMC strategy (e.g., particle MC); 2) implementing the algorithms using GPU (Graphics Processing Unit) computing (work by Fagard-Jenkin); 3) using fast approximations based on the Kalman filter (within an MCMC algorithm).

Wadden Sea Grey Seals
Summary by G. Aarts (based on Brasseur et al., 2015) During the workshop, the model described in Brasseur et al. (2015) was presented.
Grey seals were first observed breeding in the Dutch Wadden Sea in 1985, after centuries of absence, and the breeding colony there is now the largest on the European continent. Brasseur et al. (2015) describe the changes in grey seal numbers and their geographical expansion, and estimate how these processes were influenced by immigration from other colonies. Yearly counts of hauled out animals were carried out between 1985 and 2013, monitoring three different periods of the seals' annual cycle.
Using priors determined for the UK population, a Bayesian demographic model was fitted to pup numbers to estimate the population parameters driving the growth. This included immigration of subadults into the breeding population, which contributed to an average growth rate in the pup counts of 19% per annum, much higher than expected in a closed population. This immigration may account for approximately 35% of the total annual growth. In addition, at least 200 grey seals from the UK visit the area temporarily.
Recovery of the population in the Netherlands occurred more than 50 years after grey seals were protected in the UK. Therefore, conservation and management actions responding to population changes in long living marine mammals may require several decades.
One aspect of the population model that is potentially valuable for other population assessments, is the direct integration of raw pup (and population counts) into the fitted population model: where ( ) is the pup presence probability, describes the pup birth and describes the pup departure. Pup birth and departure were defined as logistic functions that contain the following parameters: mean date of birth ( ℎ ), mean pup presence duration ( ) (the period that pups are recognisable as pups), variability in pup birth 1 and finally the parameter 2 that allows for the observed forward shift in mean birth date.
The results of this model fitted to pup data are presented in Figure 3.
The model was fitted in WinBUGS and all details can be found in Brasseur et al. (2015).

Summary by M. Hammill, C. den Heyer, G. Stenson and A. Mosnier
Grey seals in eastern Canada have shown remarkable recovery from approximately 15-20,000 in the early 1960s to a recent estimate of 424,000 animals in 2016. Growth rates have slowed, but the population still appears to be increasing at an annual rate of approximately 4% per year.
Abundance is estimated using a modified version of the model developed for harp seals and is estimated by fitting a model to independent estimates of the total pup production and reproductive rates observed for seals 8 years of age and older (referred to as 8+) (Hammill et al. 2015). The model assumes that density-dependent factors are affecting demographic parameters that can be described using a theta-logistic curve acting on reproductive rates and juvenile survival. The value of theta is assumed to be 2.4. It is also assumed that the sex ratio is 1:1. The input data consist of periodic estimates of pup production from aerial surveys and mark-recapture studies, annual removals (reported catch, struck and lost, by-catch and scientific research) and annual reproductive rate data (1960present) (den Heyer & Bowen, 2017Hammill, Gosselin, & Stenson, 2017;Hammill, den Heyer, Bowen, & Lang 2017). Grey seals are managed as two separate herds. Separate models are fit respectively to pup production data from the Gulf of St Lawrence and from the Scotian Shelf (coastal Nova Scotia and Sable Island) (den Heyer, Lang, Bowen, Hammill, Gosselin et al. 2017). A single reproductive rate data set is shared by both models (Hammill, den Heyer et al., 2017).
In an earlier version, the model estimated adult mortality rate and fixed juvenile mortality as three times the adult rate. However, in most large mammal populations, adult mortality rates tend to be low and show little inter-annual variation, while juvenile mortality rates are thought to be more variable. A branding program initiated on Sable Island around 1985 and continued intermittently has provided information on adult and juvenile survival rates and shown that juvenile mortality rates are approximately 15 times adult female mortality rates. The model was adjusted with a new multiplier for juvenile rates and the model continues to estimate adult mortality rates. The mark-recapture analysis of the branded resighted animals has also shown that male mortality rates are higher than those of females (den Heyer & Bowen, 2017). As a result, the most recent estimate of abundance was adjusted post-hoc assuming a male to female ratio of 0.7M to 1 female (based on a markresight analysis suggesting lower survival rates in males) (den Heyer, Bowen, & McMillan, 2014;den Heyer & Bowen, 2017;Hammill, den Heyer et al., 2017). Land-breeding tends to be the norm for this species, but historically grey seals pupped primarily on the pack ice in the southern Gulf of St Lawrence. However, since the turn of the century, the amount of ice in the gulf in January and February has declined. In some years, mortality has been high and this has been incorporated into the model formulation. The response of grey seals has been a shift by breeding females to have their pups on small, isolated islands. There is a process underway to develop a new model to describe the dynamics of grey seals in eastern Canada. Instead of using a theta-logistic curve density-dependent, changes in pup mortality are described using a Beverton-Holt relationship. In this model, density dependent mortality of young of the year is considered either as a function of overall abundance or as a function of first-year abundance. The model was originally fitted using a least squares approach (see Northwest Atlantic harp seal model). The new model (coded in TMB, i.e., using ML to fit the model) estimates the sex and age-specific survival and age at maturity by fitting to the brand resighting data. It also fits the pup production time series well. The 2-sex model will be useful for estimating the role of grey seals in the ecosystem.

COMPARISON OF MODELS
A comparison of the different models is presented in overview form in Tables 2 and 3 below.

IDENTIFIED PROBLEMS WITH CURRENT MODELS & PROPOSALS FOR SOLUTIONS
The WKSEALS workshop held 4 breakout sessions to discuss the current data needs and availability for the different models, and to work with the models given the data and code provided on the shared Github repository.
All of the breakout groups discussed various approaches one might take to improve on current seal models, and/or any additional data that may be required to inform that approach. In terms of data, the discussions included additional life history data that may improve model performance but focussed more on various datasets on environmental conditions that may affect various life history parameters. Significant discussion was had regarding what ecological processes may improve performance if they could be incorporated into existing models, and an assessment of the availability of data on these processes, as well as the possibility to use proxies or assumptions to capture the relevant processes where data are not available.
Given that the workshop was arranged based on a recommendation from WGHARP, the group specifically discussed potential advancements to the existing assessment model for harp and hooded seals in the North Atlantic, providing suggestions for alternative model formulations and the inclusion of additional data. In addition though, the grey seal assessment models used in the UK, Canada/USA and the Netherlands were also discussed, both in terms of the potential to inform the North Atlantic assessment models, but also to identify potential improvements of the other models.

Grey seals
The following is a summary of key points made during the breakout discussions, with particular relevance for harp and hooded seals: 1.
Regarding data collection for harp and hooded seals, discussion was had regarding the timing of sampling and whether this may be improved upon. The timing of sampling was specifically discussed in relation to reproductive parameters of harp seals (and hooded seals to a lesser extent), and the question of how representative ovary samples taken during the moult are in terms of assessing actual natality in that season (1-2 months prior to samples being collected). This depends partly on how late term abortions are reflected in the ovary-based pregnancy rates in the most recent reproductive cycle. Regressing large corpora lutea/early corpora albicantia with luteinized tissue of firm texture as well as some connective tissue are considered signs of pregnancy in the preceding season, but it is not known how late an abortion may occur and still leave an ovarian corpus of this description.

2.
Including migration in the models was discussed in relation to various populations and for different reasons.
For the Barents Sea/White Sea harp seals, one way of accounting for the discrepancy between the dramatic decline in pup production and the still high fecundity is to assume movement into new pupping areas. It was suggested that looking into the use of satellite imagery to scan for (new) pupping areas outside traditional areas could be a valuable way to investigate this possibility. The CAPS project (Censusing Animal Populations from Space) and its work on Antarctic ice seal populations was mentioned as an illustrative example of how this may be approached in practice. The potential for using satellite imagery to identify areas of poor ice quality and therefore the potential for high pup mortality was also discussed.

3.
The group discussed the high estimated adult mortality in the current Barents Sea/White Sea harp seal model and how this is most likely connected to the model compensating for overly high fecundity estimates. The current model for Barents Sea/White Sea harp and hooded seals estimates mortality (rather than using data or estimating mortality independently) due to lack of survival data or population age structure data for harp and hooded seals. A key challenge, which is shared with many seal models, is therefore that the model only estimates constant mortality (M0 and M1+) over time, which does not allow for annual variations, or mass mortality events. This means that as long as fecundity is treated as fixed data (with high values that do not take late-term abortion into account), mortalities are pushed upward to compensate. Different approaches for working with the fecundity data were therefore further discussed. 4.
The current model for the NE Atlantic harp and hooded seals is deliberately 'stiff' (due to an historical desire from WGHARP not to fit an overly flexible model to sparse data). It is therefore unable to reproduce the rapid changes in pup production that have been observed for the Barents Sea/White Sea population. The population modelling shows results that are sensitive to fecundity rate and a realistic variation in fecundity could account for the observed pup production changes, but it is not certain that this is the correct explanation. The group therefore discussed the value of examining if this can be linked to changes in the Barents Sea ecosystem. Approaches to this could include: a) incorporating environmental data into the model that explain potential late-term abortions, and link this environmental relationship via a condition index to the fecundity data; or b) further develop a stochastic model that includes a late term abortion parameter as a random effect. Other suggested additions would be to test models with timevarying mortality rates (primarily for adults). One example of such a model was provided by Lars Witting towards the end of the workshop, which illustrated how a dramatic drop in adult mortality (or emigration) resulted in a vastly improved fit to the pup production data. 5.
Fecundity may be overestimated from samples, because possible late-term abortions may not be detected based on examination of corpora albicantia in the ovaries. To inform modelling approaches for the Barents Sea/White Sea harp seals that could incorporate abortion rates, it was suggested that an investigation of Canadian samples might be carried out to explore relationships between spring and winter body condition and changes in fecundity. This could help inform variation in fecundity in the Barents Sea/White Sea population, where fecundity is measured from spring samples. 6.
Exploratory work could be conducted outside of the population modelling process, to investigate statistical relationships between environmental data, seal condition, and pregnancy rates for the Barents Sea/White Sea seal species. 7.
In determining which ecosystem variables may be relevant to consider and incorporate in the models, seal diets were highlighted as crucial indicators. Looking at changes in distribution/abundance etc. of prey by examining the existing data available from ongoing fish and invertebrate surveys was proposed as one potentially useful approach. It was noted that there are other groups and experts already looking into changes in the NE and NW Atlantic and performing comparative work across them (e.g., the COARC project and work within ICES Working Groups). It was noted that there are a number of candidate environmental data sets that could be considered for inclusion as drivers of vital rates. As one example, the annual Barents Sea Ecosystem surveys carried out by IMR and the Polar branch of VINRO (PINRO, named after N.M. Knipovich) (e.g., see Protozorkevich & van der Meeren, 2020) have produced a unique dataset describing changes in this marine ecosystem over time. This dataset could therefore be further explored to provide covariate inputs into mortality and/or fecundity estimates. The group emphasised that in working in this direction, the intention would not be to do a full integrated ecosystem assessment, but rather to identify potentially relevant factors and where information may be available. It was decided that one breakout group would consider these datasets in more detail.
Although the focus was primarily on harp and hooded seals in the North Atlantic due to terms of reference for the workshop and its connection to WGHARP, the group also discussed possible future work and collaborations related to grey seals.

1.
The group specifically noted that work is in progress under OSPAR to expand the UK grey seal population model to cover all European populations and the relevance of including migration in such a multi-stock model.

2.
Given that one disadvantage of WinBUGS/Bayesian MCMC methods is that they can be quite slow, developing a TMB implementation using ML to improve model fitting time for the European grey seal model under development was discussed as valuable.

3.
There was also discussion of implementing the current Norwegian grey seal model using approaches developed for the Dutch Wadden Sea model. This exercise was seen as an important first step towards implementing a Europe-wide grey seal assessment model. Ultimately, it was suggested that it may be desired to fit such a model using the existing Bayesian Importance Sampling approach by Thomas et al. (2019), but the speed of optimisation possible with the current Maximum Likelihood approach may provide a useful alternative for scenario testing and model prototyping.

4.
Applying a similar approach to the inclusion of density dependent effects on juvenile mortality for the Wadden Sea grey seal model as is done in the UK grey seal model was also noted as valuable developmental work to pursue.

EXPLORATORY WORK ON IMPROVEMENTS PERFORMED DURING WKSEALS 2020
The workshop breakout groups were able to perform some collaborative pilot work on the identified needs for a short time during the workshop. This included exploration of model specification, data processing to allow for models to be run with different data sets, and model implementation. The initial results of these pilot efforts to advance model improvements are summarised below, together with the recommendations for the further development of these ideas that were reported to the final plenary session of the workshop.

Data & Model Exchange:
Work was done towards investigating a data and model exchange between Canada and Norway for the NE and NW Atlantic harp and hooded seal models. There was some progress made towards converting the datasets for insertion into the other models and the aim is to have this data/model exchange tested before the planned ICES benchmark meeting in December 2022. This work will therefore be continued and will be primarily undertaken by Arnaud Mosnier, Tor Arne Øigård and Martin Biuw, with assistance from Alejandro Buren, Mike Hammill and Garry Stenson.

Hooded Seals in the Greenland Sea:
The pup production of hooded seals in the Greenland Sea has shown a steady decline over the last couple of decades. During the workshop, an agestructured selection-delayed population dynamic model was developed by Lars Witting to try and reconcile the decline in pup production with estimated trends in the pregnancy rate and age of reproductive maturity.
The selection delayed dynamics are based on an ecoevolutionary model that integrates density dependence with regulation from the across generational responses to the selection pressure of the density dependent interactive competition in the population (for further explanation/ clarification see Witting, 2013). The dynamics is typically damped cyclic, with the estimated trajectory for the hooded seal in the Greenland Sea resembling the bottom phase of a cycle. This trajectory might reflect a damped cycle that was initiated by the historical exploitation in the area, and/or by a historical increase in competition following a deteriorating habitat. However, catches from earlier historical time periods are needed in order to investigate this further.
Harp Seals in the Barents Sea/White Sea: Jay ver Hoef explored fitting the Barents Sea/White Sea harp seal data to a simple 2age Leslie matrix model, using MCMC sampling to estimate parameters that are stochastic from year to year. Each parameter had a common prior across years and there were also priors on starting values for adults and pups, as well as a prior that controls density dependent decay in fecundity. The 2age Leslie matrix model can be obtained from more detailed age-and sex-structured Leslie matrix models if they are collapsed properly, the details of which can be found in the appendix of Boveng, ver Hoef, Withrow, & London (2018). The model was able to capture the drop in harp seal pup production, but is not a realistic model for the data, as there are better priors that can be developed from auxiliary data. The code runs quickly for such a small data set, even though there are 228 parameters, by using batch sampling. There is very little survey data, which is available on pups only. The results of the model depend a lot on the priors and many more prior scenarios can be tried with the R package MCMCharp, which can be found at https://github.com/jayverhoef/MCMCharp. The attraction of a full Bayesian model is that samples are drawn from the full joint posterior distribution. The posterior distribution of pup and non-pup abundance is derived from the Leslie matrix model. Moreover, the posterior distribution of intrinsic growth, obtained from the first eigenvalue of the Leslie matrix model, is easily obtained if the joint posterior distribution of Leslie matrix parameters is available.
Lars Witting also provided a working paper that used a Bayesian age-structured and density regulated model to examine the potential reason for the crash around 2005 in the Barents Sea/White Sea stock of harp seals. The modelling attempted to explain the crash from 4 different combinations of timedependent variation in adult survival and reproduction.
The first model version assumed that the crash around 2005 was due to a short period with increased adult mortality. This model did a good job of fitting to the pup counts, but did not track the drop in fecundity, with the low 2006 estimate falling below the confidence intervals of the model. The second version assumed that the crash was due to a decline in fecundity and that the pregnancy data reflects this. This model tracked the pregnancy data, but could not track the pup production counts. The third model assumed that the crash was due to a decline in both fecundity and survival, and that the pregnancy data describes the decline in fecundity. This model did a good job of tracking both the pregnancy data and the pup counts.
A decline in adult survival is, however, perhaps not the most expected response of the population. The fourth model thus assumed that the crash was due to a decline in fecundity only, with the early years of the time series having no abortion or early pup mortality (pup mortality that occurs before survey counts are made) and it was assumed that these processes became important around the time of the crash. With no explicit parameter for abortion rate, abortion was included in the fecundity term within the model. The fecundity of the model was thus not expected to correspond with the observed late-pregnancy rates, and these were therefore removed from the likelihood function. As shown in Figure 4, this model explains the trajectory of the pup counts, as well as the early pregnancy rate data, and the estimated fecundity is significantly lower than the pregnancy rates in 2006 and 2011 due to increased abortion. The high pregnancy rate in 2018 is also explained by the model, suggesting that there is no longer elevated abortion in the population. Figure 4. Harp seals in the Barents Sea/White Sea: Projected medians and 90% credibility intervals for Rep+Abo (bat) -the 4th version of the Bayesian age-structured and density regulated model, which assumed that the crash was due to a decline in fecundity only, with the years with early data having no abortion or early (before count) pup mortality, while the latter were significant around the crash. Data from ICES (2019).

Restructuring to Address Fecundity Data:
Model equations for the NE Atlantic harp seal model were examined, and some restructuring was suggested to better address the issues related to the fecundity data. A stochastic model for fecundity could be included into a model for natality that includes potential for late term abortions, informed by including environmental drivers/covariates on annual fecundity rate. Population size could be one such driver, thus including density dependence into the model.
An autoregressive process is represented by a constant correlation coefficient and stochastic 'error' term which is Normally distributed around 0.
In year = −1 + This is then incorporated into a logistic expression for pregnancy rate = exp( + + 1 1 + ⋯ ) 1 + exp ( + + 1 1 + ⋯ ) Where covariate effects can be included as additional terms in the linear predictor, e.g., 1 might be an environmental variable, or population density of seals if density dependence is to be included.
A multiplier is included to scale the estimates of pregnancy with an additional effect of abortion during later pregnancy (after the moult), giving an estimate of natality = Environmental Covariate Data: One break out group developed a list of climatological and biological data that may influence vital rates in pinniped populations, and which could potentially be used as covariate inputs to proposed models. An excel sheet overview of environmental data that are/may be relevant for different harp and hooded stocks was collated, including repositories where these environmental data are available (provided in Supplementary File 3). This table could be enhanced to make a distinction between purely historical datasets and those for which future projection modelling results are available.
During model development, it will be important to consider how useful the inclusion of environmental drivers can be for projections of the effects of management decisions for seal populations, given that future levels of these drivers are unknown. For some environmental variables, predictive modelling results are available to suggest future directions of change (e.g., IPCC for climate variables and ice cover). For others (e.g., abundance of relevant fish stocks) knowledge of future levels is limited and seal population projections might then need to be carried out across a range of plausible scenarios of future change.
TMB Implementation: During the workshop, the Wadden Sea grey seal model was 'ported' from fitting in WinBUGS to TMB. One advantage of the TMB approach is that it is much faster (1-2 seconds, instead of 3 hours), however it does not provide the full posterior distribution of parameter values. Parameter estimates appeared to be somewhat different between the two fitting approaches, so further investigation may be needed to resolve the cause of this and decide on the correct approach.
In addition to the trials and early developmental work that was done by participants within the limited timeframe of the workshop, it is worth noting again that the model information and datasets are available on the GitHub site so others are welcome to continue experimenting.

CONCLUDING REMARKS
Across species, models and areas, there are inevitably benefits and trade-offs involved in different model-fitting approaches that need to be considered during model selection and development, including fitting speed (where MLE approaches are superior) and taking account of the potential for inference based on model outputs (where Bayesian methods are strong). In addition, there may be an issue of ease-of-use and model update and maintenance, where publicly available frameworks such as TMB and WinBUGS/JAGS, STAN and others may be very useful. For the provision of scientific basis for management advice, WKSEALS participants suggested that a combination of model-fitting methods, including both a fast-fitting approach and then a slower Bayesian option with access to the complete posterior distribution for detailed inference, may be most beneficial.
It was also concluded that by incorporating environmental covariates in the population models, there is the potential both to improve our understanding of important processes, and to improve the fit of models to data. Covariates could include both abiotic and biological information (e.g., ice cover, abundance of relevant fish stocks). It was also noted that it could be relevant to consider including other apex predator species (such as large predatory fish, whale and bird species), which could be seen as indicators of overall ecosystem state, or as competitors for grey, harp and hooded seals. Environmental covariates might be associated with seal population trends in different ways in different areas and there was general agreement that modelling should continue to make allowance for this.
A concern was raised that WGHARP provides the basis for scientific advice based on projections of seal population trends, and this raises questions of how to handle the high levels of uncertainty involved in environmental projections and how to present management options on the basis of this. It was recommended to update Supplementary File 3 (collating information about environmental data sets) to include information about the availability of predictive modelling results for the environmental variables listed, and the time scale and spatial scale/resolution for which predictions may be available.
The group also concluded that working to include environmental covariates in population modelling without moving into full integrated ecosystem assessment was important. Trying to understand the whole ecosystem was recognised as beyond the scope of working groups like WGHARP and may not offer significant help with population predictions either. For population assessment, the interest is more in the use of ecosystem proxies for data in the model. Improving relations between the ICES integrated ecosystem assessment groups (such as WGIBAR for the Barents Sea) and WGHARP may be a valuable way to build relations and exchange knowledge/expertise.

PATHWAYS FORWARD
Revised models for harp and hooded seals in the Northeast Atlantic (Barents Sea/White Sea & Greenland Sea Harp Seals & Greenland Sea Hooded Seals) are to be developed on the basis of the discussions and experiments from WKSEALS 2020. These will be further advanced throughout 2021 for presentation, testing and validation at the ICES benchmark meeting in 2022. To support this aim, further work on the pilot initiatives identified and begun at WKSEALS is planned. Work will also continue on the useful proposals and initial work done to improve and expand the grey seal models, although this will occur outside the frame of the planned benchmark meeting for harp and hooded seals. The group also intends to advance both lines of work through broadening and building the network of seal modelling expertise with regular correspondence, meetings, and shared coding efforts (using the GitHub repository established for the workshop, which will remain accessible through the ICES system).
The specific objectives for model development on harp and hooded seals agreed at WKSEALS include: 1.
Run the Barents Sea/White Sea harp seal data through other models, particularly the Canadian assessment models. Since the Canadian models work with more frequent fecundity estimates, this will require the development of an approach to address missing years of fecundity data in the Barents Sea/White Sea harp seal data sets. In the current NE Atlantic assessment models, this is done by simple linear interpolation between periods of observed fecundity. A similar approach can be used as a first approach, unless the Canadian model is able to deal with missing data in a similar way to the Bayesian approaches demonstrated by Jay ver Hoef and Lars Witting.

2.
Investigate the potential to improve on estimates of harp seal abortion rates, using Canadian samples to improve on methodology, which could then be applied to Barents Sea/White Sea samples. 3.
Modify model structure for the Barents Sea/White Sea harp seal models to include a 'late' abortion term that scales estimates of fecundity to produce more realistic estimates of natality. This should be based on the existing stochastic version of the model (presented in Øigård & Skaug, 2015). 4.
Test a model for the Barents Sea/White Sea populations of harp and hooded seals that also allows interannual variation in mortality, to account for potential mass mortality or emigration events. 5.
Modify model structure for both the Barents Sea/White Sea harp and hooded seal populations to include effects of density dependence on vital rates. 6.
Explore the potential to modify model structure for the Barents Sea/White Sea harp and hooded seal populations to include the effects of environmental drivers on seal vital rates/carrying capacity.

Recommendations
• Strengthen relations, communication and collaboration between WGHARP and ICES regional integrated assessment and ecosystem modelling communities. • Support the future development of this seal population modelling community within ICES and NAMMCO, potentially reaching out to a wider group of marine mammal population modellers.

•
Build links between the ICES marine mammal working group and the NAMMCO coastal seals working group by presenting the WKSEALS outputs at their meetings.

•
Organise an ICES/NAMMCO workshop on seal population modelling at the biennial marine mammal conference of the Society for Marine Mammalogy in Florida in 2021, or alternatively, arrange an online event to continue the work and network building.