Model Introduction
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. The river creates a dynamic coastal river plume that juvenile yearling Chinook salmon (Oncorhynchus tshawytscha), ‘smolts,’ must navigate during the earliest stage of their marine migration (17). Previously, the plume was hypothesized to be a refuge from predators and a rich feeding ground (18), but more recent research suggests that the plume is predator-rich and with lower survival in the plume than adjacent estuarine and ocean habitats (19), (20), (21), (22), (23), (24). Due to the perceived importance of the early marine period, 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 (25), (26), (27), (24).
Determining the effects of changing 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 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 (28).
The second strategy, navigating to minimize opposing current while migrating north, would 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. (29) 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 (30), and may exploit counterflowing eddies (31). Juvenile salmon are believed to similarly exploit turbulent flow during their downstream migration (32), (33).
Individual-based models are a cost-effective means of evaluating which strategy produces 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 (34) and uses simulation output (salinity, temperature, and current) from the Virtual Columbia River model developed by the Center for Coastal Margin Observation and Prediction. The IBM doesn’t 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 detections of acoustic tagged juvenile yearling Chinook (Oncorhynchus tshawytscha) and oceanographic simulations from the 2009 spring outmigration. This dataset was chosen because tagged smolts were released from early-April through May, covering a wide range of spring ocean conditions.
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. Smolts detected on both the Astoria and Willapa Bay sub-arrays (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 (29) and (35).
Individual-based Model Description
The model was implemented in NetLogo and is described here using the ‘Overview, Design concepts, and Details,’ or ODD, protocol of Grimm et al. (36).
Purpose
The purpose of the model was to compare two hypothesized strategies smolts may use during their northward migration through the Columbia River plume, i.e., do smolts select habitat to maximize growth, or attempt to minimize opposing currents to more quickly transit through plume region?
Entities, state variables, and scales
There were two entities in the model, ocean cells and in-silico smolts. Each ocean cell had six state variables: salinity, temperature, current velocities (x- and y-), water depth, and a representation of feeding conditions. The two types of in-silico 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, weight, optimal swimming speed, heading, and a binary variable indicating whether the smolt was in estuarine (PSU < 27) or ocean (PSU ≥ 27) waters.
Simulations 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.5km 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 1. At each time step, eight smolts of each type 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 at time step 1250. At each step, the ocean cells updated their salinity, temperature, current velocities, and prey availability. Subsequently, each smolt calculated its weight- and temperature-dependent optimal swimming speed, evaluated whether it was in estuarine water (cell salinity < 27 PSU) or marine water (cell salinity ≥ 27 PSU) and moved accordingly (see below). Once all smolts executed their move, the model stepped forward.
Design Concepts
Pearcy (37) is widely credited with the hypothesis 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 and wind-driven currents and other oceanographic processes (38). Survival in this region is presumably affected by the migratory strategy adopted by smolts, and insight into the manner that 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 (39).
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:
- 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 (40).
- 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 (41), flow (30) (42), temperature (43), and salinity (44). Their presence in trawls is also positively correlated with chlorophyll concentration and plankton abundance (45), (46), (47).
- 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 (48), (49), (50).
- Smolts are typically found in the top 12 m of the water column (51) and vertical migrations are assumed not to significantly affect horizontal progress, allowing the model to be collapsed into 2 spatial dimensions.
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 (52).
Initialization
The model was initialized with flow environment variables, water depth, and the simplified representation of the prey field specified from (53). Initial smolt fork lengths were assigned by drawing randomly from fork lengths of tagged smolts detected at Astoria and Willapa Bay in 2009. The 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 smolts travel time (d) from release to plume entry, and 1.05 (mm d-1) is an observed daily growth rate in the Columbia River (49).
Input Data
The flow environment data 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 (CMOP) 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 (54).
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 (51) finding that smolts in the plume region are found in the upper 12 m of the water column.
Submodels
Submodel 1: 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). The model was fitted to the data using the nonlinear least squares function provided in the R stats package (52).
Submodel 2: Optimal Swimming Speed
Optimal swimming speed was calculated using the formula described in (55) and parameterized in (56):
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 3: Bioenergetics
In bioenergetics models, growth over time is estimated using a simple mass balance equation, growth = consumption – (respiration + egestion + excretion). This submodel uses the Wisconsin bioenergetics equation sets (57); Table 1) for consumption (eqn. 3), respiration (eqn. 1), egestion and excretion (eqn. 2). The equations are parameterized with values from the literature on Pacific salmon bioenergetics (Table 2). Prey energy density, Qf, is a single global parameter. All smolts use this submodel to grow. The growth-seeking model smolts use it to identify and move at a metabolically optimal speed towards habitat (grid cells) where growth opportunities are greatest.
Sensitivity analysis revealed the bioenergetics submodel to be sensitive to two parameters, prey energy density and proportion of maximum consumption, which is consistent with the findings of similar analysis reported in (55) and (58) (following (59)), as well as (60). When the bioenergetics submodel was validated by simulating migrations in the model world with bioenergetics parameters from Table 2, model smolts grew at rates consistent with a range of reported growth rates for juvenile yearling Chinook observed in the Columbia River plume region (61).
Submodel 4: Prey
A submodel representing declining prey availability with distance from shore was specified from Peterson et al.’s (53) description of three cross-shelf zones. In this sub-model, prey is represented to the model fish through the proportion of maximum consumption (bioenergetics parameter ‘p_Cmax’; (62), values of which are drawn from one of three distributions representing inshore waters (depth ≤ 50 m, mean ‘p_Cmax’ = 0.8), mid-shelf waters ( 50 m < depth ≤ 150 m, mean ‘p_Cmax’ = 0.6) and outer shelf/open water (depth > 150 m, mean ‘p_Cmax’ = 0.02). Although the assignment of ‘p_Cmax’ is consistent with the biomass estimates in (53) and assumptions about the volume of water searched by a transiting salmonid (63), they are a characterization, rather than a replication, of actual feeding conditions.
Submodel 5: Movement rules
In estuarine waters (PSU<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 six hours after first exposure to marine waters (PSU≥27). The timing of the switch in movement rules was calibrated using the median plume residence times of the model smolts, defined here as the median time between release and first detection on the Willapa Bay subarray. Calibration experiments indicated that a 6 h delay provided the closest match between observed plume residence times (median = 4.3 d) and plume residence time of growth maximizing smolts (median = 4.2 d) and opposing current minimizing smolts (median = 4.1 d). Model results were insensitive to variations in the switch time of ±50%.
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 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. In the estuary, smolts were not permitted to swim east (upriver) of the model boundary.
Submodel 6: 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 meters. 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 2011 season were recorded, but removed prior to analyzing the cross-shelf distributions against observed values.
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 (64), where the null hypothesis is that there is no difference in the distributions. Syrjala’s (64) test was calculated using the R package ecespa, and summary circular statistics using circular.