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 1,370 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 (Figure 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 (°C), current velocities (x- and y-; m s-1), 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-1), 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-hour steps. The model world origin was in the top-left grid cell, at position 3601 255000mE, 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 x 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 Figure 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 1,250 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.
Design Concepts
Eight of eleven 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 (50).
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:
- 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 (51).
- 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 (52), flow (41) (53), temperature (54), and salinity (55). Their presence in trawls is also positively correlated with chlorophyll concentration and plankton abundance (56), (57), (58).
- 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 (59), (60), (61).
- Smolts are typically found in the top 12 m of the water column (62) 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
where FLtag is the measured fork length of the tagged smolt at the time of tagging, Triv is the tagged smolt’s travel time (d) from release to plume entry, and 1.05 (mm d-1) is an observed daily growth rate in the Columbia River (60). 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 (63).
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 (64). 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) (65).
DB-22 contains flow data at 90-second 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 fifteen-minute 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 (62) 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:
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 (66) and parameterized in (67):
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 (68) (Table 1) for consumption (their eqn. 3), respiration (their eqn. 1), egestion and excretion (their eqn. 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.
Submodel D: Prey
A submodel representing declining prey availability with distance from shore was specified from Peterson et al.’s (69) 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_Cmax; (70), values of which are drawn from one of three beta distributions representing inshore waters (depth ≤ 50 m, beta(α=5, β=1), mean p_Cmax ~ 0.8), mid-shelf waters ( 50 m < depth ≤ 150 m, beta(α=5, β=4), mean p_Cmax ~ 0.6,) and outer shelf/open water (depth > 150 m, beta(α=.1, β=5), mean p_Cmax ~ 0.02). Although the assignment of p_Cmax is consistent with the biomass estimates in (69) and assumptions about the volume of water searched by a transiting salmonid (71), 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 d 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 (70). 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 150mm and their weight was 38.8 grams. 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. (72) (Figure 3).
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 (73), where the null hypothesis is that there is no difference in the distributions. Syrjala’s (73) test was calculated using the R package ecespa.