A model to illustrate the potential pairing of animal biotelemetry with individual-based modeling

Abstract

Background

Animal biotelemetry and individual-based modeling (IBM) are natural complements, but there are few published examples where they are applied together to address fundamental or applied ecological questions. Existing studies are often found in the modeling literature and frequently re-use small datasets collected for purposes other than the model application. Animal biotelemetry can provide the robust measurements that capture relevant ecological patterns needed to parameterize, calibrate, and assess hypotheses in IBMs; together they could help meet demand for predictive modeling and decision-support in the face of environmental change.

Results

We used an simple exemplar IBM that uses spatio-temporal movement patterns of 103 acoustic-tagged juvenile yearling Chinook salmon (Oncorhynchus tshawytscha), termed ‘smolts’, to quantitatively assess plausibility of two migratory strategies that smolts are hypothesized to use while migrating north through the plume of the Columbia River (United States of America). We find that model smolts that seek to maximize growth demonstrate movement patterns consistent with those of tagged smolts. Model smolts that seek to move quickly out of the plume region by seeking favorable currents do not reproduce the same patterns.

Conclusions

Animal biotelemetry and individual-based modeling are maturing fields of inquiry. Our hope is that this model description and the basic analytical techniques will effectively illustrate individual-based models for the biotelemetry community, and perhaps inspire new collaborations between biotelemetry researchers and individual-based modelers.

Background

Individual-based models (IBMs), which are also referred to as agent-based models, track discrete, autonomous individuals with static or dynamic state variables, attributes, and behavior, simulate interactions between them and their environment, and record the emergent effects on populations and ecosystem ecology [1, 2]. IBMs have been used in a wide variety of applied and theoretical studies within ecology, as well fields such as epidemiology, where representing individual variation was required or desired [2, 3]. For example, IBMs have been used to explore the effects of conservation and management interventions on habitat connectivity of bighorn sheep [4], and elucidate factors that drive population dynamics of Antarctic salps [5]. They are also emerging as a promising tool for predictive modeling and decision-making in the face of environmental change [6,7,8].

IBMs do present some limits and drawbacks, which have been well described. Developing an IBM may present researchers with a steep learning curve in several fields of science, as well as scientific programming [9]. The complexity of such IBMs can also make them difficult to describe and communicate [10], although a regularly updated protocol for describing IBMs in publications exists [11]. Applied IBMs, in particular, may require extensive research, or a large extant body of literature, in order to be parameterized [12]. The numerous bioenergetics parameters in the IBM presented here, for example, could be parameterized with literature values only because Chinook salmon (Oncorhynchus tshawytscha) have been the subject of extensive research investments by the government, academia, and non-profit sectors. While new remote sensing and animal tracking technologies will increase the volume of available data [13], data may remain a limitation for species and regions that are difficult to sample. Relatedly, fluid dynamic models of the ocean or atmosphere that may be used for studying migrations of birds or marine organisms are often available, but their spatio-temporal resolution and skill in remote or particularly complex physical environments may be limited. IBMs are often computationally intensive, and although the emergence of commercial cloud computing may overcome physical constraints on computation and storage, computational barriers may still present themselves as issues of cost. Beyond these implementation challenges, individual-based modeling remains a research topic in its own right, with current challenges identified by An et al. in transparency and reusability, validation, software and big-data, human decision-making, and theory [14].

Nonetheless, given that parameterizing, calibrating, and validating individual-based models, and assessing the hypotheses they are intended to address, require robust measurements that capture relevant ecological patterns [15, 16], animal biotelemetry, which focuses on the individual and has provided rich insights into survival [17], spatio-temporal movement [18], and habitat selection [19, 20], should be a natural complement to individual-based modeling, as noted by Byron and Burke [21]. There is an opportunity for expanding and improving the use of animal biotelemetry in IBMs. In five representative examples of the pairing of biotelemetry and IBMs, described below, two use biotelemetry data for parameterization and calibration and three extend the use of biotelemetry data to include assessing alternative hypotheses. Four of the models draw on small biotelemetry datasets, and all re-use extant biotelemetry data. Re-use can be fruitful, but may be limiting when biotelemetry data collections would ideally be designed to ensure an IBM can address the questions researchers wish to answer, not just ones they can answer with the available data.

Nabe-Nielsen et al. [22] and Liukkonen et al. [23] used biotelemetry data for IBM parameterization and calibration. Nabe-Nielsen et al. in a modeling framework for predicting the impacts of anthropogenic disturbances on animal movement and fitness, used statistical and visual techniques to calibrate a correlated random-walk submodel with data from a single dead-reckoning tagged North Sea harbor porpoise (Phocoena phocoena), as well as dispersal parameters from 25 satellite-tagged porpoises. Biotelemetry data were not used to corroborate the model. Similarly, Liukkonen et al. used satellite telemetry data to simulate movement patterns and characterize home range formation and spatial ecology of Saimaa ringed seal (Phoca hispida saimensis), with the stated intent to later extend the model to study the effect of changing environmental conditions and different conservation scenarios on population dynamics. They used satellite biotelemetry data from one seal to parameterize their model and verified the model with movement data from five other seals. They also found they needed to calibrate a distance-from-home parameter separately for each of the five seals used for verification, highlighting how larger telemetry samples collected for the purpose of developing such an IBM may be helpful.

Bauduin et al. [24] and Ohashi and Sheng [25] extended the use of biotelemetry data to assess alternative hypotheses about the mechanisms underlying important ecological processes. Bauduin et al. used VHF-based biotelemetry data from 35 Atlantic-Gaspésie caribou (Rangifer tarandus caribou) tracked between 1998 and 2001 to parameterize their IBM model and assess alternative movement hypotheses, including random walk, biased correlated random walk, foray loops to reproduce caribou extra-range movement patterns, and caribou fidelity during mating season. Ohashi and Sheng assessed ten alternative swimming strategies, including various passive, random, and active directed movements, of juvenile Atlantic salmon (Salmo salar) migrating in the Gulf of St. Lawrence. They assessed their hypotheses by comparing travel times from juveniles tagged in 2009 and 2010 against travel times of simulated salmon implemented in an ocean circulation model with numerical particle tracking. Although the original dataset included 93 tagged juvenile Atlantic salmon released into the St-Jean River, travel time was measured from only three smolts that were subsequently detected at distant subarrays of acoustic receivers in the Gulf of St. Lawrence.

Morrice et al. [20] also extended the use of biotelemetry data to assess alternative hypotheses. They investigated potential migration pathways of yearling and subyearling Chinook salmon (Oncorhynchus tshawytscha) in the Columbia River estuary using alternative movement models; passive drift, random walk, biased correlated random walk, and taxis for yearling Chinook and passive drift, random walk, kinesis, and area search for subyearlings. These movement models, and the migration pathways that emerged, were assessed against detections of tagged yearling and subyearling Chinook. Morrice et al.’s [20] model is notable for the large sample of tagged animals used: 3880 yearling Chinook salmon and 4449 subyearling Chinook salmon tagged and released in 2010, and the small tags that made such large sample sizes possible are an example of the technological advancements identified by [13] as aiding individual-based modeling.

These five studies suggest that a fuller realization of the promise of IBMs and biotelemetry would include biotelemetry researchers and individual-based modelers working together at the outset of a research project. Premier et al. [26] allude to this when they observe that improvements in spatio-temporal resolution of telemetry technology could enable contemporary re-analyses to advance animal movement ecology and form the basis of future simulation studies. One reason why biotelemetry data are not more widely used in individual-based modeling may be limited cross-fertilization between these fields, as evidenced by the re-use of biotelemetry data in the preceding examples and publication of three of these papers in Ecological Modelling, a specialized journal whose audience may have limited overlap with the biotelemetry research community.

In an effort to highlight the application of animal biotelemetry to individual-based modeling for the biotelemetry community, we present here a simple IBM that draws on a large acoustic telemetry dataset and uses quantitative methods to evaluate hypotheses regarding the migration of juvenile yearling Chinook salmon (smolts). The motivation for this study is the decline in populations of Chinook salmon in the Columbia River basin, which have resulted in five major population groups being listed as threatened or endangered under the U.S. Endangered Species Act [27]. Chinook salmon migration studies, such as this one, are one facet of long running research efforts to understand and address the causes of their decline.

The Columbia River basin drains nearly 671,000 square kilometers of land, entering the Pacific Ocean at the border between the USA states of Oregon and Washington (Fig. 1) [28]. The river creates a dynamic coastal river plume that juvenile yearling Chinook salmon, termed smolts, must navigate during the earliest stage of their marine migration. Previously, the plume was hypothesized to be a refuge from predators and a rich feeding ground [29], but more recent research suggests that the plume is predator-rich, resulting in lower survival in the plume than adjacent estuarine and ocean habitats [30,31,32,33,34,35]. Due to the perceived importance of the early marine period of salmon migration, there has been a persistent interest in the interaction between the Columbia River plume environment and juvenile salmon migration and survival, particularly as to how migration and survival may be affected by river dynamics [35,36,37,38].

Fig. 1

Regional map. Dashed line delineates the model world boundaries (distortion due to re-projection). Solid black lines mark the subarrays of acoustic receivers at Astoria and Willapa Bay. Blue contour lines mark the 200-m and 500-m isobaths

Our model uses spatio-temporal movement patterns of 103 acoustic-tagged smolts to quantitatively assess two migratory strategies smolts could use while migrating north through the plume of the Columbia River (USA). These strategies include (a) maximizing growth, or (b) selectively using local currents to minimize plume occupancy time, and thus exposure to predation in the plume region. A quantitative assessment of the patterns in cross-shelf distribution of in silico smolts against cross-shelf distribution from the telemetered smolts reveal that the strategy of maximizing growth is plausible, but not the selective use of coastal currents to speed passage through a predator-rich region. This model description, and the basic analytical techniques, illustrate individual-based models for the biotelemetry community and can serve as a starting point for collaborations between biotelemetry researchers and individual-based modelers.

Methods

Model introduction

Determining the effects of changing Columbia River plume conditions on smolt migration and survival requires some understanding of the strategies smolts may employ during this phase of their migration. Two conceivable strategies are (a) to select habitat that maximizes growth or (b) navigate so as to minimize opposing currents during northward migration. The first strategy would be expected to allow smolts to more quickly reach a size large enough to reduce availability to size-selective predators while contributing to the growth necessary for them to survive their first winter at sea. This strategy is also consistent with the critical size, critical period hypothesis which posits that the adult return rate of salmon is influenced first by predation-driven mortality at ocean entry, then by starvation-driven mortality during the following winter that results from a failure to build minimum energy reserves [39].

The second strategy, navigating to minimize opposing current while migrating north, could speed northward passage, reducing the period of exposure to predation in the plume region (assuming that survival rates were lower in the plume than farther north). Brosnan et al. [40] demonstrated that survival of groups of yearling Chinook released to migrate through the plume is negatively related to their travel time, so this strategy could be beneficial. This behavior has also been observed in the field; stereovideographic studies in tributaries of the Fraser River reveal that adult sockeye salmon (Oncorhynchus nerka) migrating upriver selectively travel in lower current [41], and may exploit counterflowing eddies [42]. Juvenile salmon are believed to similarly exploit turbulent flow during their downstream migration [43, 44].

Individual-based models are a cost-effective means of evaluating the plausibility of these strategies by examining whether they produce results consistent with observed migration patterns. The model presented here was implemented in NetLogo v.5.2, a popular Java-based software for individual-based modeling [45] and uses numerical model output (salinity, temperature, and current) from the Virtual Columbia River simulation, which was developed by the Center for Coastal Margin Observation and Prediction. The IBM does not test the effects of different strategies on plume survival, but rather which, if any, of the plausible strategies reproduces migratory patterns observed in acoustic telemetry studies of juvenile yearling Chinook—a first step in addressing the survival question.

Tagging and telemetry

The model draws on oceanographic simulations and detections of acoustic tagged juvenile yearling Chinook released from early-April through May 2009. VEMCO V7-2L acoustic tags were surgically implanted in 1370 hatchery-raised Columbia River Basin juvenile yearling Chinook smolts. Acoustic receivers were deployed in subarrays across the river at Astoria Bridge and across the continental shelf at Willapa Bay to detect the passage of tagged smolts (Fig. 1). Smolts detected on both the Astoria and Willapa Bay subarrays (N = 103) were used to calculate plume residence time, and the cross-shelf distribution at Willapa Bay; the relatively small sample size of 103 tracked fish is a consequence of mortality losses and the need for each included smolt to be detected on both subarrays. Greater detail on the tagging procedures and acoustic array used to track tagged smolts can be found in [40] and [46].

Individual-based model description

The model was implemented in NetLogo and the model description follows the ODD (Overview, Design concepts, Details) protocol for describing individual- and agent-based models [10], as updated by [11].

Purpose and patterns

The purpose of the model is to evaluate two hypothesized strategies smolts may use during their northward migration through the Columbia River plume. Specifically, do smolts select habitat to maximize growth, or attempt to minimize opposing currents to more quickly transit through the plume region? We use a quantitative test to compare patterns in the cross-shelf distribution of two types of model smolts, each using one of these hypothesized strategies, against the distribution of acoustic-tagged smolts detected on an array of cross-shelf receivers.

Entities, state variables, and scales

There were two entities in the model, ocean cells and smolts. Each ocean cell had five state variables: salinity (dimensionless), temperature (℃), current velocities (x- and y-; m s⁻¹), depth (m), and a representation of feeding conditions (proportion of maximum consumption). Salinity and depth are used for navigating to the ocean and preventing smolts from ‘beaching’ or crossing dry land. Current velocities are used for habitat selection and movement. Feeding conditions are used for habitat selection and growth. Temperature is used for growth and swimming speed.

The two types of virtual smolts are defined by their migratory strategy, maximizing growth (MaxGro) or minimizing the current opposing their northward migration (MinCur). Each smolt, regardless of strategy, had five state variables: fork length (mm), weight (g), optimal swimming speed (cm s⁻¹), heading, and a binary variable indicating whether the smolt was in estuarine (salinity < 27) or ocean (salinity ≥ 27) waters. These salinity thresholds are derived from Horner Devine’s [47] conceptual model of the Columbia River plume, which defines coastal waters with salinity ≥ 27 as the final mixing level for estuarine and ocean waters.

The simulation initiated on April 16th, 2009 and ran for 70 days in 1-h steps. The model world origin was in the top-left grid cell, at position 3601 255000 mE, 369500mN (Oregon State Plane coordinates, North American Datum of 1927) and the model extended 100 km east of the origin and 200 km south in a grid with 0.5 km × 0.5 km cells. Grid cell size is the approximate distance that a 130 mm smolt (the smallest size tagged) would swim in one hour at 1 body length per second. This model world encompassed the lower Columbia River estuary and plume region, including the arrays of acoustic receivers at the Astoria Bridge and Willapa Bay.

Process overview and scheduling

A line diagram illustrating the model processes can be found in Fig. 2. At each time step, eight MaxGro and eight MinCur smolts were introduced into the model. To permit the last smolts introduced into the model sufficient time to migrate past Willapa Bay, the introduction of model smolts ended after 1250 time steps (0200 on June 7th, 2009). At each step, the ocean cells updated their salinity, temperature, current velocities, and representation of feeding conditions. Subsequently, each smolt calculated its optimal swimming speed, evaluated whether it was in estuarine water (cell salinity < 27) or marine water (cell salinity ≥ 27), moved according to its decision rules (see below), and grew. The model stepped forward once all smolts moved and grew.

Fig. 2

Process overview diagram. The origin is at the top left blue box entitled “Start”. Purple boxes delineate the process, orange boxes delineate the algorithm for ensuring smolts do not cross or move to land, and white boxes indicate submodels

Design concepts

Eight of 11 individual-based model design concepts identified in [11] appear in this model, including basic principles, emergence, adaptation, objectives, prediction, sensing, stochasticity, and observation. Learning, interaction, and collectives are not included.

The basic principle of this model, which underpins the model design concepts of [11] applied here, is the notion widely credited to Pearcy [48] that the year class strength of returning salmon is established during the early marine life history, including the period of plume residency. Yearling Chinook departing the Columbia River must travel through the river plume, the dynamics of which are affected by hydropower-regulated river discharge, tides, wind-driven currents, buoyancy, and other oceanographic processes [49]. Survival in this region may be affected by the migratory strategy adopted by smolts, and insight into the manner in which individual smolts interact with the plume environment while migrating, revealed by patterns in their distribution across an acoustic array, could be used to evaluate the impact of changing oceanographic and river conditions.

The cross-shelf distribution of smolts during their northward migration is the emergent characteristic of this model, resulting from a combination of the prey field and ocean currents. The prey field drives the choice of habitat that MaxGro smolts orient towards, and, along with temperature, mediates the swimming speed of MaxGro and MinCur smolts via their growth. The ocean current is a physical forcing variable that affects both MinCur and MaxGro smolts and drives the choice of habitat that MinCur smolts orient towards.

In this model, smolts adapt to changes in themselves and their environment by selecting a weight and temperature-dependent swimming speed and orienting towards habitat consistent with their migratory strategy. Implementing these behaviors requires the following assumptions about sensing:

• 1. Smolts are assumed to employ negative rheotaxis (orientation in the direction of current) to navigate through the model estuary and into the plume. This assumption is grounded in a finding that upregulation of the hormone thyroxine in response to changing light conditions stimulates negative rheotaxis, leading smolts to travel downstream [50].

• 2. Smolts are assumed to have a compass sense and the ability to detect gradients in temperature, salinity, current, and prey availability. The specific mechanisms by which smolts sense these variables are not modeled, but field and laboratory observations indicate that smolts can sense gradients in the orientation and intensity of magnetic fields [51], flow [41, 52], temperature [53], and salinity [54]. Their presence in trawls is also positively correlated with chlorophyll concentration and plankton abundance [55,56,57].

• 3. Smolts are assumed to begin a directed northward migration shortly after reaching the ocean, an assumption that is strongly supported by decades of U.S. and Canadian coastal trawl surveys [58,59,60].

• 4. Smolts are typically found in the top 12 m of the water column [61] and vertical migrations are assumed not to significantly affect horizontal progress, allowing the model to be collapsed into 2 spatial dimensions.

Grimm et al. [11] distinguish adaptive behaviors as direct or indirect objective-seeking. In the former, individuals make decisions based on alternatives ranked according to some measure of how they contribute to meeting an objective. In the latter, they simply follow rules that reproduce observed behaviors. Orientation by negative rheotaxis in the estuary and compass sense (ocean) are rules that reproduce observed patterns and are characterized as indirect objective-seeking. Orientation toward habitat, to maximize growth (MaxGro smolts) or minimize opposing current (MinCur smolts), are direct objective-seeking behaviors whose proximate objective is to grow rapidly or speed through a predator-rich region, respectively, and ultimately to survive and return to reproduce. Orientation along these gradients is an implicit prediction per [11] and there is no explicit prediction in the model.

Stochasticity appears in the model at several points. First, each smolt initialized in the model is assigned a fork length drawn randomly from the starting fork lengths of 2009 tagged smolts detected at Astoria and Willapa Bay in 2009. The assigned fork length (FL) at tagging was adjusted for growth between release and detection at Astoria by

$$ {\text{FL = FL}}_{{{\text{tag}}}} + 1{\text{.05 * }}T_{{{\text{riv}}}} , $$

where FL_tag is the measured fork length of the tagged smolt at the time of tagging, T_riv is the tagged smolt’s travel time (d) from release to plume entry, and 1.05 (mm d⁻¹) is an observed daily growth rate in the Columbia River [59]. Second, the direction smolts iteratively search to avoid beaching or crossing dry land is a random binomial process with equal probability of turning left or right. Finally, the proportion of maximum consumption (p_Cmax) is assigned at each model step by submodel D (see below), which draws on beta distribution sampling.

At every time step, the unique identification number, type, size, weight, local flow variables, heading, optimal swimming speed, and model coordinates (including latitude, longitude, and the corresponding grid cell) of each smolt were recorded. All model output analyses were conducted in R [62].

Initialization

The model is initialized with current, salinity, temperature and water depth assigned from simulation input data, and prey field from submodel D (see below). Smolts are initialized at the Astoria array with fork lengths drawn from input data. Their initial weight is calculated from length in submodel A and swimming speed from submodel B (see below).

Input data

The flow environment variables and depth for each ocean cell at each time step was derived from DB-22, a database of Virtual Columbia River (VCR) model simulations. The VCR was built by the Center for Coastal Margin Observation and Prediction using SELFE, an open-source, community-supported code designed for the effective simulation of 3D baroclinic circulation in the Columbia River estuary and plume that uses semi-implicit finite-element/volume Eulerian–Lagrangian algorithms to solve the Navier–Stokes equations on unstructured triangular grids [63]. Published skill assessments have demonstrated that the VCR models reliably reproduce the characteristics (current, salinity, temperature) and dynamics of the Columbia River estuary and plume [49, 64].

DB-22 contains flow data at 90-s intervals from the surface to the seafloor and recorded for multiple depths. To keep the computational demand tractable, flow data at 4 m, 8 m, and 12 m for each cell in the model were extracted from DB-22 in 15-min intervals and then averaged and reformatted in R to create single hourly values readable by NetLogo. The choice of depth intervals is based on Emmett et al.’s [61] finding that smolts in the plume region are found in the upper 12 m of the water column.

Submodels

Submodel A: length–weight regression

Length–weight conversions were made using the regression model:

$$ W{\text{ = e}}^{{{ - 14}{\text{.075 }}}} {\text{ * FL}}^{{{3}{\text{.514}}}} , $$

where W is weight (g) and FL is fork length (mm). This is an empirical model fitted to Columbia River basin hatchery-origin yearling Chinook smolt length–weight data collected by NOAA researchers trawling at three transects in the Columbia River plume region (Columbia River, Grays Harbor, and Willapa Bay) in May and June 2008–2011 (C Morgan, NOAA, personal communication, 2013).

Submodel B: optimal swimming speed

Optimal swimming speed was calculated using the formula described in [65] and parameterized in [66]:

$$ U_{{{\text{opt}}}} {\text{ = ACT * }}W^{0.13} * {\text{e}}^{{{0}{\text{.0405*}}T}} , $$

where T is temperature and a state variable of the cell and W, weight, is a smolt state variable. ACT is a parameter from the bioenergetics submodel described below.

Submodel C: bioenergetics

In bioenergetics models, growth is estimated using a simple mass balance equation, growth = consumption–(respiration + egestion + excretion). This submodel uses the Wisconsin bioenergetics equation sets [67] (Table 1) for consumption (their Eq. 3), respiration (their Eq. 1), egestion and excretion (their Eq. 2). The equations are parameterized with values from the literature on Pacific salmon bioenergetics (Table 2). These equations are used to grow smolts at each time step and calculate growth potential in each cell at each time step to inform movement rules of the MaxGro smolts.

Table 1 Bioenergetics equations

Table 2 Bioenergetics submodel parameters

Submodel D: prey

A submodel representing declining prey availability with distance from shore was specified from Peterson et al.’s [68] description of three cross-shelf zones. In this submodel, prey is represented to the model fish through the proportion of maximum consumption (bioenergetics parameter p_C_max [69], values of which are drawn from one of three beta distributions representing inshore waters (depth ≤ 50 m, beta(α = 5, β = 1), mean p_C_max ~ 0.8), mid-shelf waters (50 m < depth ≤ 150 m, beta(α = 5, β = 4), mean p_C_max ~ 0.6,) and outer shelf/open water (depth > 150 m, beta(α = 0.1, β = 5), mean p_C_max ~ 0.02). Although the assignment of p_C_max is consistent with the biomass estimates in [68] and assumptions about the volume of water searched by a transiting salmonid [70], they are a characterization, rather than a replication, of actual feeding conditions.

Submodel E: movement rules

In estuarine waters (salinity < 27), all model smolts align with the current in their cell (negative rheotaxis). The time between exposure to marine waters and the switch to one of the two plume migration strategies occurred 7 h after first exposure to marine waters (salinity ≥ 27). In marine waters, smolts set their heading towards the cell north of their position with maximum growth opportunity or minimum southward flowing current that it could reach in one hour at its optimal swimming speed. Smolts moved to the terminal point of a vector that was the addition of their movement vector (heading, optimal speed) and the current vector in the cell they occupied. If the smolt’s final position was outside the north, south, or west bounds of the model, the smolt was considered to have permanently emigrated, otherwise it grew according to the bioenergetics submodel. At the eastern boundary of the model, which could be reached in the estuary, smolts remain at the boundary until conditions result in westward travel.

A simple algorithm was applied at each movement to prevent a model smolt from ‘beaching’ itself or crossing dry land to reach water. Each smolt’s path and final destination were verified wetted. If the destination or path included dry land, smolts first searched in 45 degree heading increments for a clear path, then reduced their swimming speed and searched again. This process iterated until a clear path was found, or the smolt would move to an adjacent cell offering greatest growth or minimum south current. The direction of the heading change was chosen randomly, but only once so that heading increments were sequential.

Submodel F: simulated detections

Detections of model smolts on the Willapa Bay subarray receivers were simulated by evaluating smolt paths (Euclidean lines drawn between smolt positions at each model step) to determine if they overlapped any of the detection zones centered on each receiver. The estimated detection range of the acoustic tags was 400 m. For each smolt path that intersected a receiver’s detection radius, the receiver number, smolt number, and time of detection were recorded. Model detections on virtual receivers corresponding to those that were lost during the 2009 season were recorded, but removed prior to analyzing the cross-shelf distributions against observed values.

Bioenergetics model validation

Two simulations of fish growth were performed to ensure the bioenergetics model reasonably simulated observed growth of yearling Chinook smolts at ocean entry. First, growth was simulated over 45 days at a range of fixed sea surface temperatures (8, 10, 12, 14, and 16 °C) using the bioenergetics parameters in Table 2 and a fixed proportion of maximum consumption, 0.6, derived from [69]. To evaluate growth under dynamic conditions, 200 MaxGro smolts were introduced into the model (190 at random ocean locations and 10 at random lower estuary locations) and the model was run for 45 days. In both simulations, the initial length of the smolts was 150 mm and their weight was 38.8 g. Virtual smolt growth in both simulations were compared with juvenile yearling Chinook smolt growth rates observed in the Columbia River plume region by Tomaro et al. [71] (Fig. 3).

Fig. 3

Mean growth rates of model smolts (diamonds) compared with growth rates for juvenile yearling Chinook observed in the Columbia River plume region (circles) reported in Tomaro et al. [71]

Movement rule switch time calibration

The timing of the switch in movement rules from negative rheotaxis was calibrated by running the model with switch times of 4, 5, 6, 7, and 8 h and selecting the switch time where median plume residence times of the model smolts, defined here as the median time between release and first detection on the Willapa Bay subarray, most closely matched the median travel time of tagged smolts transiting from the Astoria subarray to the Willapa subarray.

Model output analysis

Model output analyses, including simulated detections, were conducted in R. The cross-array distribution of detections of tagged smolts and model smolts was compared using a modified Cramer von-Mises test [72], where the null hypothesis is that there is no difference in the distributions. Syrjala’s [72] test was calculated using the R package ecespa.

Model results

Virtual smolt growth in the bioenergetics model validation simulations and the individual-based model fell within range of juvenile yearling Chinook smolt growth rates observed in the Columbia River plume region and, in each case, were closer to the higher range of growth rates observed by Tomaro et al. [71] (Fig. 3). The movement rule switch time experiments revealed that a 7 h delay provided the closest match between plume residence time of MaxGro smolts (median = 4.38 d) and MinCur smolts (median = 4.29 d) and observed plume residence times (median = 4.29 d). Results of the cross-shelf distribution tests were not sensitive to changes in the movement rule switch time.

At a significance level of 0.05, there was no difference in the distributions of the MaxGro smolts and tagged smolts (p-value = 0.45; Fig. 4). There was a significant difference in the distributions of the MinCur smolts and tagged smolts (p-value = 0.00009; Fig. 4). MaxGro smolts exhibit a more diffuse pattern of cross-shelf distribution when compared with MinCur smolt tracks, which appear more bifurcated with pronounced nearshore and offshore tracks (Fig. 5).

Fig. 4

Cross-shelf distribution of detections on the Willapa Bay subarray of MaxGro (top panel, orange bars) and MinCur (bottom panel, yellow bars) smolts overlaid on detections of live, tagged smolts (both panels, white bars). Results in the top left corner of each panel are from a modified Cramer von-Mises [72] test for different distributions. A bi-modal distribution, also visible in Fig. 5 is evident. The MinCur smolts exhibit a more inshore distribution that may reflect use of a northward flowing nearshore plume current

Fig. 5

Heat map showing MaxGro and MinCur smolt paths. The solid black line marks the array of acoustic receivers at Willapa Bay. Axes are Oregon State Plane North coordinates (North American Datum of 1927)

The model plume was oriented southwest during the early and mid-points of the simulation, and northward and alongshore otherwise. Migratory paths of the model smolts were sensitive to coastal currents. 78% of model smolts migrated across a line extending from the shore to the detection radius of the outermost acoustic receiver on the Willapa Bay subarray (e.g., ‘on-shelf’), while the remaining 22% of ‘off-shelf’ smolts passed west of the outermost receiver’s detection range, or exited the model at one of the boundaries before migrating north of the line of receivers (Fig. 6). The ‘off-shelf’ model smolts experienced a mean current magnitude and direction of 1.0 m s⁻¹ (SD = 0.19) at 255 degrees True (SD = 19.5), while the ‘on-shelf’ smolts experienced a mean current of 1.2 m s⁻¹ (SD = 0.48) at 182 degrees True (SD = 68.2). As (Fig. 7) illustrates, the mean current direction of the ‘on-shelf’ model smolts tends to mask a strongly bi-modal pattern with east and west components. The initialization time of the ‘off-shelf’ smolts is also strongly bi-modal. Most ‘off-shelf’ smolts entered the model in either the first several days, or in the second month, whereas the initialization time of on-shelf smolts is more uniformly distributed (Fig. 7). The model plume was oriented southwest during the early and mid-points of the simulation, corresponding to the periods when off-shelf smolts entered the model (Fig. 7), and northward and alongshore otherwise.

Fig. 6

Heat map showing paths of ‘on-shelf’ smolts migrating across a line extending from shore to the detection radius of the outermost acoustic receiver on the Willapa Bay subarray and ‘off-shelf’ smolts either passing west of the outermost receiver’s detection range or exiting the model at one of the boundaries before migrating north of the subarray. Axes are Oregon State Plane North coordinates (North American Datum of 1927)

Fig. 7

Mean magnitude and direction of transport (top) and model initiation time (bottom) of ‘on-shelf’ smolts (left) migrating across a line extending from the shore to the detection radius of the outermost acoustic receiver on the Willapa Bay subarray and ‘off-shelf’ smolts (right) either passing west of the outermost receiver’s detection range or exiting the model at one of the boundaries before migrating north of the subarray

The differences in current experienced between MagGro and MinCur model smolts in the ‘off-shelf’ and ‘on-shelf’ groups were negligible. The ‘off-shelf’ MaxGro smolts experienced mean current of 1.03 m s⁻¹ (SD = 0.19) at 255 degrees True (SD = 20.5) and the MinCur smolts 0.98 m s⁻¹ (SD = 0.19) at 255 degrees True (SD = 18.9). The ‘on-shelf’ MaxGro smolts experienced mean current of 1.17 m s⁻¹ (SD = 0.48) at 189 degrees True (SD = 71.4) and the MinCur smolts 1.24 m s⁻¹ (SD = 0.48) at 174 degrees True (SD = 63.7).

Discussion of the individual-based model

The results from the individual-based model suggest that smolts may pursue a strategy of maximizing growth upon beginning their northward migration. Under the critical size-critical period hypothesis, which posits that adult returns are affected first by predation at marine entry, and then starvation during the first winter, this strategy could improve the probability of survival if smolts grow to exceed the gape size of predators and attain sufficient energy reserves to avoid winter starvation. Interestingly, while incorporating a representation of coastal feeding conditions acted to contain the MaxGro smolts on the shelf, 19% of the MaxGro smolts still migrated outside the bounds of the subarray at Willapa Bay. As Pearcy [48] notes, this could have implications for early marine survival. Columbia River basin smolts driven off-shelf at ocean entry might experience reduced predation pressure, particularly from seabirds nesting near the river mouth [30, 35]. Their feeding opportunities would likely also be reduced [55, 68], but directed northward migration would lead them to quickly regain the shelf environment, which arcs westward north of the Columbia River. Outmigrating yearling Chinook would be placed to take advantage of this offshore transport and potentially reduce their risk of predation since the necessary conditions occur in early spring when the transition to the upwelling season is beginning and winds are strong.

Trawl surveys and acoustic receiver subarrays terminate near the shelf break because few or no fish are caught in trawls near the shelf break. Results from this model suggest that there could be a distribution of catches/detections that would not be detected under the current sampling regime. Burke et al.’s [73] individual-based model of salmon migration revealed a similar pattern, although they describe using an Ornstein–Uhlenbeck process that oriented smolts towards the historic centers of mass of juvenile salmon sampled across the shelf in an attempt to correct it. Unfortunately, we can’t distinguish whether our IBM is missing a direct objective-seeking behavior that would further constrain model smolts to the shelf, which is captured by Burke et al.’s indirect objective-seeking behavior, or if the pattern that emerges in both models exists and hasn’t been observed.

While it might seem that the selective use of alongshore current by model smolts employing the MinCur strategy would speed passage out of the plume region and enhance survival, this behavior proved to reduce travel time to Willapa Bay by only a few hours relative to the MaxGro strategy. Applying the survival equation used in [40] to estimate survival to the Willapa Bay subarray using median travel times suggests that survival would be virtually indistinguishable between the growth maximizing smolts and current minimizing smolts (0.59 v 0.60, respectively) that reach the array. While this assumes that predation pressure on both smolt types is similar during their transit to the Willapa Bay subarray, it suggests there may be little benefit to the active use of favorable currents in the plume region.

In addition to the potentially limited benefit of selective use of coastal currents in the plume region, there are additional reasons why this pattern may not extend beyond the river. First, the high density of salmon in the river raises the costs of competing for limited food resources and salmon may restrict their feeding and perform the river migration on a limited energy budget [50, 74]; but see [75] for evidence that smolt feeding increases in the lower Columbia River. Thus, there may be strong pressure to minimize the cost of migration in the river. Conversely, the early marine environment represents a richer feeding ground where a delay in feeding could result in smolts failing to meet the energy requirements to survive their first winter at sea [39].

Second, the confined river environment provides turbulent flow and relative motion cues from the riverbed and similar features that salmon demonstrably respond to [41, 52]. Beyond the tidal plume, these cues are weakened. During the spring, current velocity in the lower Columbia River is approximately 1 m s⁻¹ (and can be much greater) and exceeds 3 m s⁻¹ in the tidal outflow at the river mouth, whereas ambient coastal current velocity is approximately 0.1 m s⁻¹ [36, 76]. Shear turbulence is likely reduced in the coastal ocean, and smolts are not positioned to use the seabed as a frame of reference [61]. Conceivably, they can detect the fine-scale changes in the intensity of turbulence and/or their motion relative to celestial objects and the earth’s magnetic field. This would allow them to take advantage of small, favorable gradients in coastal current. However, they would require very sensitive absolute and differential thresholds to these cues, which have not been determined experimentally.

An additional outcome of the model was that negative rheotaxis was sufficient to guide smolts out of the estuary. The initial expectation was that the tidally reversing flow of the Columbia River would significantly delay, or prevent, progress towards the ocean, whereas McInerney [77] reasoned that preference for increasing salinity could serve as an orientation mechanism for smolts in the estuary. However, whether salmon use salinity to orient remains unresolved [21] and the two mechanisms were compared during model development. Both produced similar migratory patterns, although orientation by negative rheotaxis resulted in smolts initially traveling further offshore under west transport conditions than did orientation by salinity. Negative rheotaxis was selected for the final model because it has a known biochemical basis wherein increasing day length stimulates the upregulation of thyroxine, which in turn causes smolts to orient with the current [51].

The successful transition from river to ocean requires that smolts complete a number of changes (i.e., smoltification). The precise trigger, or series of triggers, during this complex process that result in the behavioral transition to directed northward migration in the ocean have not been determined. It may be driven by hormonal changes prompted by environmental cues, much as upregulation of thyroxine triggers downstream migration [50]. We attempted to capture the transition via a simple threshold-delayed response model, but a mechanistic understanding of the process might explain why smolts are occasionally detected migrating at least as far south as Cascade Head, OR, before turning and swimming north [78]. Additional model development could be used to elucidate the factors that result in this behavior, as well as distinguish the means by which smolts navigate out of the estuary.

Model smolt paths suggest that migratory paths through the plume region are influenced by the winds that cause the Columbia River plume to adopt either a north alongshore orientation or southwesterly offshore orientation [47]. When the plume is oriented north and along the coast, smolts may begin migrating north nearer the coast, and further offshore when the plume is oriented toward the southwest. Figures 4, 5 and 6 show distinct nearshore and offshore model smolt paths on the shelf that suggest the strength and orientation of the plume affect how far offshore smolts are as they migrate north. This pattern is more distinctive in the MinCur smolts, likely because they seek habitat defined by currents that persist and vary over longer time scales, and thus take a narrower path, while the MaxGro smolts are more diffusely spread across the shelf as they migrate north seeking habitat that changes stochastically at each model step. The distribution of the MaxGro smolts matches that of tagged smolts, which indicates that, in practice, a smolt’s distance offshore may quickly become independent of where they started migrating north after exiting the plume. This finding suggests that timing the release of hatchery-raised smolts, or regulating discharge to influence plume dynamics, may only have limited effect on smolt survival if there is strong relationship with their migratory path in the plume.

This highlights an important limitation of this IBM, which is that predation is not modeled. Thus, while we can evaluate plausibility of the two strategies for migrating through the plume based on their ability to reproduce a migratory pattern, we cannot model which strategy results in higher survival. This was a deliberate choice based on limited availability of predation-related data, a common limitation of IBMs, but there are approaches that could be applied to this model to include predation in future. Railsback et al. [79] proposed a condition-based movement rule for stream fish in IBM’s, termed expected survival, where fish move to habitat where the probability of surviving non-starvation mortality risks, multiplied by the probability of surviving starvation risk over a pre-determined time horizon, is greatest. Yearling Chinook in the marine environment appear to respond to feeding conditions during northward migration, so applying Railsback et al.’s [79] expected survival, using the first-winter critical period as the time horizon, may be a fruitful development of the IBM presented here. This would require collaboration among researchers who have collected data on salmon predators, salmon prey, and early marine survival to develop a more accurate representation of predator and prey fields. Fortunately, such data are now becoming available (e.g., [35]). Railsback et al. [79] also proposed expected reproductive maturity, an extension of expected survival wherein model fish also consider how close they are to their size of first reproduction. This rule could be considered in future model developments, but may be more appropriate for modeling later life stages (or other species) where model fish are not seeking to reach a critical size to survive their first winter at sea.

Conclusion

The foregoing model is intended to illustrate the pairing of animal biotelemetry data and individual-based models by highlighting simple ways that biotelemetry data can be used to better inform the development of an IBM. Conditions appear ripe for a productive pairing of individual-based modeling and animal biotelemetry. First, although the body of literature describing IBMs implemented with biotelemetry data is relatively small and relies largely on repurposing of existing small datasets (but see [20] for re-use of a much larger biotelemetry dataset), it is a promising sign of active interest among individual-based modelers in using biotelemetry data. Second, national and international data assembly, archiving, and access efforts such as the U.S. Animal Tracking Network (atn.ioos.us), the Ocean Tracking Network (oceantrackingnetwork.org), and the International Movebank animal tracking database (www.movebank.org) will make more data discoverable and accessible. This will create opportunities to aggregate datasets for robust IBM parameterization, calibration, and hypothesis testing. Still, the fact that animal biotelemetry data are expensive to collect and require experience and specialized skills for most elements of array design and deployment, tagging, and data processing and analysis will still require partnerships for prospective projects. Finally, as individual-based modeling and animal biotelemetry mature, many of the challenges of implementing a successful biotelemetry-supported IBM activity have been identified. Researchers now have access to a wealth of methodological information in specialized journals and textbooks [1, 80] and a multitude of proven software tools for modeling and analysis.