Table of Contents

  1. Introduction
  2. Methods
  3. R Code
  4. Application

Introduction

Alewives

Anadromous fish migrate from the ocean into freshwater to reproduce. Populations of anadromous fish species in Maine have been in decline for centuries, primarily due to fisheries exploitation, degraded water quality, and inaccessible spawning habitat (Saunders et al. 2006).

Species historically present along the northeast coast of North America included Atlantic salmon (Salmo salar), American shad (Alosa sapidissima), alewife (A. pseudoharengus), blueback herring (A. aestivalis), and sea lamprey (Petromyzon marinus). Of these species, alewives are the most abundant and have the highest potential for population restoration. Studies have suggested that this species has a remarkably high reproductive potential, meaning that a small number of returning adults can produce a large number of offspring (Gibson and Myers 2003).

Alewives are native to the east coast of North America and range from South Carolina to Newfoundland (Durbin et al. 1979; Figure 1). Alewives are largely iteroparous in the northern part of their range, meaning that adults return to the ocean after spawning and are able to spawn several years in a row (Bozeman & Van Den Avyle 1989). In contrast, a large percentage of spawning alewives die after completing their run in the southern part of their range. This life history approach is called semelparity.

Figure 1. Historic Range of Alewife Populations



Alewife Life Cycle

It is thought that alewives live in the ocean over the continental shelf until they reach sexual maturity after 3-5 years, when they migrate to freshwater to spawn (Davis & Schultz 2009). Adults prefer to spawn in slow moving water such as lakes or flowages (Bozeman & Van Den Avyle 1989). After hatching, juvenile alewives spend 2-7 months in freshwater, and at 23-100 mm total length they migrate to the ocean (Richkus 1975).


Figure 2. Alewife Life Cycle

Figure 2. Alewife Life Cycle


Alewives were historically numerous throughout New England, but populations have been in decline for decades thanks to water quality issues, the construction of dams, high rates of ocean mortality, changes in offshore food availability, and myriad of other influences (Hall et al. 2012). However, populations in many rivers in Maine have recently demonstrated increasing abundance, including the Penobscot and Kennebec Rivers which saw a rise from tens of thousands to 1.3 million and 3.5 million returning adults in 2017, respectively (DMR 2017). These increases followed projects that focused on fishway improvement and dam removals.



Dam Passage

While upstream migration can be affected by natural barriers as well as extreme high and low water levels (Meixler et al. 2009), it is primarily dam construction that has led to declines in alewife populations through a number of both direct and indirect effects. These man-made barriers can prevent fish from reaching high-quality spawning habitat either because fish don’t survive when trying to move past a dam or because the presence of the dam causes the fish to be delayed when moving upstream (Hall et al. 2010; Pess et al. 2014).


Effects of delays at dams

Upon reaching a dam, fish may spend a lot of time swimming around and looking for the fish passage structure that allows them to move upstream. This means that an individual fish may use a lot of energy just trying to get past a dam (Castro-Santos & Haro 2003; Jonsson et al. 2010). As alewives migrate both upstream and back downstream in a short period of time, they often don’t eat much while they are in freshwater and so have to rely on the energy they already have stored in their body during migration. If a fish is delayed at a dam, it may use most of its stored energy to move upstream and not have as much as it needs left over to spawn and migrate back downstream. This is especially a problem when alewives have to migrate upstream past multiple dams before reaching their spawning habitat.

This delay can affect whether adults survive their spawning run, and some fish may not be able to lay their eggs in high quality habitat. When eggs are laid in poor quality habitat, this can affect how large and how quickly juveniles are able to grow while they are in freshwater, which in turn can determine if they survive their migration back to the ocean (Noonan et al. 2012; Hall et al. 2010).


Effect of dams on natural stream flow

In addition to blocking fish from moving upstream or downstream, the presence of dams can alter the natural stream flow. Anadromous fish are well adapted to their particular environment, which is why when they are ready to spawn adult fish generally return to the same river where they were born (reference?). Sediment can build up behind a dam and be prevented from moving downstream (Gregory et al. 2014). This can affect the structure of the stream below the dam. In addition, building dams on a stream can also cause an increase in water temperature and can reduce seasonal high flows that are a natural occurrence (Gregory et al. 2014). This is a problem because many anadromous fish species use changes in temperature and flow as a cue to know when to move upstream or downstream.

Alteration to the natural stream flow can also prevent nutrients from moving downstream, which reduces freshwater productivity (Hall et al. 2012). Dams also affect freshwater species, and when impoundments form above a dam there is often a shift in the fish species that are present (Gregory et al. 2014; Jono’s paper?). For example, bass prefer slow-moving water and their numbers can be high in an impoundment. Bass eat both juvenile and adult alewife in large numbers (reference), and having more bass present can reduce the alewife population.


Fishways

Passage structures (which will be referred to as “fishways”) can differ in how well they work for moving fish past a dam. Many different types of structures have been built, and depending on size of the dam they have different lengths and steepness (Noonan et al. 2012; Cooke & Hinch 2013). Also, fishways may be easier for a fish to find and pass through at certain flow rates than at others, making success differ throughout the season (Meixler et al. 2009).

Many fishway types favor only one or two species, and it is hard to design a fishway that can pass all the fish that are migrating. Anadromous species can differ quite a lot in how well they can find the fishway entrance, how high they can jump, and how long they are able to swim without getting tired. This means that the type of fishway built at a dam can affect how well a species is able to recover because it determines how well an individual is able to move past a dam (Meixler et al. 2009; Bunt et al. 2012; Noonan et al. 2012). For example, alewives have a hard time swimming for long periods of time, and so a longer fishway will be harder for them to use.

In addition, each type of fishway has a daily limit in the maximum number of fish that can be moved upstream. At the peak of the spawning run, the number of fish that pass over the dam may be much lower than the number of fish that approach the dam in their effort to move upstream. This can delay fish during their migration, and as mentioned above can lead to an individual using a lot more energy to move upstream than if the dam wasn’t present.


The trouble with multiple dams

While passage efficiency through one dam may be high, many individuals have to pass through multiple dams to reach upstream spawning habitat. This has a additive effect that reduces overall passage efficiencies at a watershed scale (Brown et al. 2013). The passage rate at an individual dam may be high, as is seen in Figure 3 below, but the more dams a fish has to pass the lower the probability of success gets.

Figure 3. Effect of multiple dams on watershed-scale passage rates

Figure 3. Effect of multiple dams on watershed-scale passage rates


This same additive effect is seen when adults and juveniles move back downstream through multiple dams. Low downsteam passage efficiency can also hinder population growth and recovery. For iteroparous species such as alewives, the success of both upstream and downstream passage is important for population growth as it affects fecundity, which is the ability to produce a lot of offspring. When population fecundity is low, fewer eggs are being produced overall. Population fecundity can be reduced when adults attempt to come back for multiple years to spawn, but each year have a low probability of surviving migration both upstream and back downstream.

Figure 4 illustrates the effect that dams have through time on an alewife population. In this example, we are following one age group through three years of migration.

  • The first year (Age-3), we start with 1,000 fish in the age group. Due to the levels of passage at four dams, 500 of these fish make it to the spawning habitat.
  • Assuming no fish die due to spawning, these 500 fish then attempt to migrate back downstream, but only 328 make it through all four dams back to the ocean.
  • The second year (Age-4), the 328 fish in this age class again migrate upstream to spawn. We are assuming none of these fish died in the ocean during the year, which while not realistic serves to makes the example more straightforward.
  • Of these fish, 170 make it successfully to the spawning habitat, and 111 make it back to the ocean after spawning.
  • The third year (Age-5), these 111 fish migrate upstream, with only 57 successfully entering the spawning habitat.
  • Of these 57 fish, only 37 make it back to the ocean.

So as you can see, after only three years our age class has been reduced from 1,000 fish to 37 fish, and keep in mind that only half of these fish are female and are laying eggs.


Figure 4. Effect of multiple dams through time

Figure 4. Effect of multiple dams through time


Dam Passage and Age Structure

In the absence of dams, alewives will continue to spawn up until they are 8 or 9 years old. However, when dams are present alewife populations lose these older fish, which are usually larger and produce more eggs than younger fish (Leggett et al. 2004; Oldani et al. 2007). Downstream passage has more of an effect on the age structure of a spawning run than upstream passage, as can be seen in Figure 5. Each bar represents the abundance of an age class of spawners, going from fish that are 3 years old to 8 years old.

  • When looking across a row from left to right, downstream passage (DS) is held at the same percentage and upstream passage (US) is increased from 20% to 80% successful passage.
  • Increasing upstream passage increases overall abundance within each age group, but does little to affect the presence of older age groups.
  • The opposite occurs when moving down a column, where upstream passage is held constant and downstream passage increases from 20% to 80%.
  • As downstream passage increases, more age groups are present within the population, and the abundance of older age groups increases.
Figure 5. The effect of adult passage on age structure. DS = downstream passage, US = upstream passage. Passage rates are listed as percent success.

Figure 5. The effect of adult passage on age structure. DS = downstream passage, US = upstream passage. Passage rates are listed as percent success.


Deterministic population model

The alewife population nmodel was developed to compare theoretical spawner abundance between scenarios with different dam passage rates. Spawner abundance was calculated using a deterministic population model. This type of model defines inputs using averages applied to groups. A deterministic model is used to explore general trends and compare the results of scenarios when different average values are used as inputs. For example, we can compare the difference in spawner abundance at a dam when the average rate is 80% successful passage versus when it is 90% successful passage, as shown below:

Figure 6. Deterministic model

Figure 6. Deterministic model


What we cannot do with a deterministic model is make forecasts or predictions about the exact number of spawners that will be present in the river after a certain number of years. This is because a deterministic model does not produce an estimate of environmental variability. In order to estimate environmental variability within the model, we would need to run multiple scenarios using a range of values for each possible input, as shown in Figure 7.


Figure 7. Estimating environmental variability

Figure 7. Estimating environmental variability


These scenario outputs would then be used to estimate an average value for the population as well as variation around that average. This type of modelling approach needs a lot of data because we need a realistic range of values associated with each input. This amount of data was not available for alewife, particularly in the St. Croix River, which is why a deterministic modelling approach was used.

Both of these approaches are informative ways to model a population, as long as we keep in mind the assumptions associated with each. For the deterministic alewife population model, the following assumptions must be kept in mind when interpreting the output:


Model Assumptions

  • There was no annual environmental variability built into the model
    • Inputs are averages
  • All spawning habitat was considered to be of the same quality
  • All density-dependent mortality is included in the recruitment curve
  • In the absence of dams, fish distribute throughout the system according to habitat availability
    • But this will change and there will be several distribution options

Again, this model was developed to compare management scenarios, not to make detailed predictions about the exact number of alewife that will return.



Management Questions

So what types of questions can be answered using this alewife population model? Ones that make comparisons between scenarios, such as the ones listed below:

  1. How could improving passage at a specific dam increase total alewife abundance in the river?

  2. Where would passage improvements give us the largest increase in spawner abundance?

When given the option to improve passage at multiple dams, we can use this model to see which dam would be the best to focus on and give you the biggest bang for your buck.

  1. How does removing or building a dam increase or decrease alewife abundance?

Within the model, passage can be set to 100% success to represent a dam removal. This is essentially saying that the dam has no effect on fish moving upstream. On the other hand, setting passage to 0% success would represent an impassable barrier. A barrier will have a stronger or weaker effect on total spawner abundance depending on its location in the river.


Back to Top

Methods

Spawner abundance is calculated using a deterministic population model. This type of model defines inputs using averages for group values. This type of model is not meant to make forecasts or predictions, but rather to explore general trends and how outputs change given a new set of input values. One major assumption of this type of model is that there is no estimate of environmental variability around an input.

Inputs and data sources:

Inputs were based on the life history of anadromous alewife. Fish were moved through the life cycle using a stepwise annual progression. Stocks in the ocean entered freshwater to spawn and produce young of the year (YOY). A proportion of these survived to “graduate” into juveniles which then became part of the ocean population, completing the cycle (Figure 8). Estimated values from the ocean population from the previous year were then used in the current year to estimate inputs for the rest of the life cycle, and this was repeated for each year the model was run.

Figure 8. Simplified diagram of alewife life history.

Figure 8. Simplified diagram of alewife life history.


The model was organized so that inputs were estimated for alewife populations in the northernmost part of their range, as many were based on a meta-analysis that was conducted using populations from Rhode Island to Nova Scotia (Gibson 2004; Gibson and Myers 2003). These inputs are shown in the table below under “Forward-projecting population model”. This was an age-structured model, meaning that some values were estimated for each age class and combined to get a total for the population.

When exploring a particular river, parameters could then be derived for a specific alewife population using information that may be available such as mass-at-age, fecundity, and the probability of maturing. For the rest of the methods section, I will describe the calculations for alewife populations in general and then use the St. Croix River as an example alewife population. Data for each age class was taken from both the St. Croix Milltown trap alewife data (Fisheries and Oceans Canada et al. 1981-2016) as well as from the literature. For the online application, all population-specific values default to those taken from the St. Croix River. However, these values can be changed when setting up scenarios for different rivers within the state of Maine.


Table 1. Population inputs, including those taken from the literature and those estimated from the St. Croix (SC) Milltown Trap Data (1981-2016)

Table 1. Population inputs, including those taken from the literature and those estimated from the St. Croix (SC) Milltown Trap Data (1981-2016)

Back to Top

Mathematical Components of the Model

General equations: Forward-projecting population model for alewife

The forward-projecting population model consisted of a series of equations for each life history stage. Egg production for a given year t, \(N_{eggs,t}\), is calculated using the number of females that survived to spawn multiplied by the fecundity relationship as below: \[\Large N_{eggs,t} = \sum_{\alpha=3}^{8} [S_t(e^{-M_{spawn*0.25}})*0.5]*F_a\]

in which \(S_t\) is the total number of fish in the spawning population in year t (described below), \(M_{Spawn}\) is the probability of mortality associated with the spawning run prorated for 3 months (or 0.25), 0.5 is the assumed female:male ratio, and \(F_a\) is the fecundity relationship (Table 1). For this model, mortality associated with the spawning run was assumed to occur before spawning such that only females that survived contributed to egg production.

The number of juveniles produced in a year class was modelled as a density dependent process, which was characterized using a spawner-recruit (SR) relationship. The choice of SR curve can affect the dynamics of the recruitment rate as the spawning population increases (Elliott 1985; Needle 2002; Subbey 2014). There are many different types of SR curves used in fish population modelling (Hilborn and Walters 1992), and the Ricker curve has been previously used to explore the Alewife SR relationship (Tommasi et al. 2015). However, the Beverton-Holt (BH) curve was used for this model because Gibson (2004) found it provided a better fit to the data for eight Alewife populations in the northern part of their population range than did the Ricker, and the data available for these populations were not sufficient to fit a three parameter model.

An example of the BH curve is illustrated below in Figure 9. The slope of the curve is determined by the lifetime reproductive rate, or \(\alpha\). The curve approaches a plateau estimated by the asymptotic recruitment level, or \(R_{asy}\). At a small number of spawners, the slope of the curve is steep and each spawner has a high reproductive rate. This means each spawner is producing a large number of surviving juveniles. As the number of spawners increases and the plateau is approached, the per capita reproductive rate decreases. This means each spawner is producing fewer surviving juveniles because the carrying capacity of the system is being approached. As that capacity is approached, which in the figure is represented by the red dotted line, an increase in spawners does not results in more surviving juveniles.

Figure 9. Beverton-Holt spawner-recruit curve

Figure 9. Beverton-Holt spawner-recruit curve

The BH curve was used to model a density-dependent relationship in the population model by tying egg production to juvenile production (\(J_t\)) as follows: \[\Large J_t = \frac{\alpha N_{eggs,t}}{1 + \frac{\alpha N_{eggs,t}}{R_{asy}}}\]

Here, juvenile abundance was calculated for year t based on the total number of eggs for year t \((N_{eggs,t})\), the asymptotic recruitment level (\(R_{asy}\)), and the maximum number of juveniles given the average fecundity per unit mass at the origin of the SR relationship (\(\alpha\)).

The population of immature fish in the ocean was divided into age classes between Age-0 and Age-8 fish, each with an associated instantaneous mortality rate for fish in the ocean, \(M_{ocean}\), and probability of maturing at that age, \(m_a\) (Table 1). The ocean population was linked to juvenile production in freshwater by setting the number of Age-0 fish in the ocean population in year t, \(O_{0,t}\), equal to the number of juveniles produced in that year:

\[\Large O_{0,t} = J_t \]

Abundance of immature fish in other age classes was calculated by projecting the abundance forward using the mortality rates and maturity probabilities: \[\Large O_{a +1,t+1} = O_{a ,t}e^{-M_{Ocean}}(1-m_{a +1})\]

Immature fish between Age-2 and Age-7 also had a probability of maturing, which allowed them to enter the spawning run the next year (Ages 3-8), with survival occurring between spawning year classes. For age a and year t, first-time spawners \((S_{\alpha +1,t+1,0})\) and repeat spawners \((S_{a+1,t+1,p+1})\) that spawned p times previously were calculated as follows: \[\Large S_{a +1,t+1,0} = O_{a ,t}e^{-M_{Ocean}}(m_{a +1})\]

\[\Large S_{a +1,t+1,p+1} = \varphi S_{a ,t,p}e^{-[M_{Ocean}*0.75+M_{Spawn}*0.25]}+(1-\varphi )S_{a ,t,p}e^{-[M_{Ocean}]}\] Each year class had an associated probability of spawning \((\varphi)\), which was separate from the probability of maturing and allowed those individuals that did not successfully spawn to return to the ocean.

Mortality rates were used in several population equations, and total natural mortality for mature spawners was split into ocean mortality (\(M_{ocean}\)) and spawning mortality (\(M_{spawn}\)) in order to estimate carcass nutrient inputs. The instantaneous natural mortality rate for adults was reported as 1 by Gibson and Myers (2003a) and an average rate of 0.7 was reported by the ASMFC (2012), and so an instantaneous rate of 0.85 was used in this study. A range of interval spawning mortality was reported in both Kissil (1974) and Durbin (1979). The average of those reported values (45%) was used to calculate an instantaneous spawning mortality rate.

The ocean mortality rate was then calculated as the difference between total mortality and spawning mortality rates (Table 1). Ocean mortality was iteratively adjusted based on age-class proportions and probabilities of spawning such that the product of \(M_{ocean}\) and \(M_{spawn}\) was equal to the total mortality rate. \(M_{ocean}\) was applied as an annual rate to both immature fish and the small percentage of individuals that did not successfully enter the river to spawn \((1-\varphi)\). Only the ocean mortality rate was used to project their abundance forward, whereas a higher mortality rate associated with spawning, \(M_{spawn}\), was included for fish that spawned successfully. \(M_{ocean}\) was prorated to 9 months and \(M_{spawn}\) to 3 months to reflect the timing of the Alewife spawning run.
The total number of fish in the spawning population in year t \((S_t)\) was calculated as: \[\Large S_t = \sum_{a,p}S_{a ,t,p}\varphi\] The number of spawners was used to calculate egg production, thereby closing the loop for calculating population dynamics associated with each portion of the Alewife life history.

Back to Top

Site-specific equations: Parameters derived for the St. Croix River

The equations in the forward-projecting population model required the derivation of multiple parameters. Morphometric information collected from the St. Croix River (1981-2016) was used as inputs for mass and fecundity (Table 1; St. Croix Milltown Trap 1981-2016). For mass-at-age, an average mass was calculated for each age class combining males and females. Maturity rates were averaged from Gibson and Myers (2003a, Table 1).

Spawner-recruit parameter derivation was taken from the Alewife model developed by Gibson (2004) based on multiple Alewife populations in northeastern North America. Parameters of the SR curve were adjusted based on the total habitat size in the St. Croix River.

Four parameters were used in the calculation of the egg deposition from the number of spawners. The probability of spawning was kept constant for all ages at 95%, the sex ratio was assumed to be 50%, and spawning mortality was calculated for 3 months as described above. Fecundity slope and intercept were calculated using a linear regression with average mass-at-age and corresponding average gonad mass-at-age recorded from the St. Croix River (St. Croix Milltown Trap 1981-2016).

For each female with a recorded ovary mass, the total egg mass was calculated by subtracting a spawned gonad mass from an unspawned gonad mass with the assumption that the former represented the mass of just the organ itself. The spawned gonad weight was an average mass from 14 downstream migrants from the St. Croix River. A total egg count for each fish was then calculated by multiplying the total egg mass by 7,890 \(eggs*g^{-1}\) (Kissil 1974). This led to an average of ~130,000 eggs per female (range: 5,050 - 305,659 eggs). The total number of eggs per fish was then regressed against mass to get a slope and intercept that was used to estimate average egg production for each age class. This was done to account for differences in fecundity with the age of the fish when calculating the total number of eggs produced in a given year. A linear regression function provided the best fit to the data (\(R^2\) = 0.66; y = 871.72x - 50916) when compared with an exponential (\(R^2\) = 0.4711), a logarithmic (\(R^2\) = 0.63), and a power function (\(R^2\) = 0.51), though the difference was small.

Juvenile production involved an estimate of density-dependent survival using a Beverton-Holt SR relationship. Two parameters were derived for this equation in addition to egg deposition: \(R_{asy}\) and \(\alpha\), neither of which are available for the St. Croix River Alewife population. These two parameters were both calculated from the results of a meta-analysis of the dynamics of Alewives based on eight populations in the northern part of their distribution ranging from Rhode Island, USA to Nova Scotia, Canada (Gibson 2004). Following the approaches of Myers et al. (1999) and Myers et al (2001), Gibson (2004) standardized the data prior to analysis in order to produce probability distributions for the maximum lifetime reproductive rate \((\tilde{\alpha})\), and the asymptotic recruitment levels in terms of the spawning population size \((\tilde{R}_{asy})\). For his analysis, Gibson defined the age-of recruitment as Age-3. For our analysis, the relationship was rescaled from Age-3 recruits and SSB to juveniles and egg production. \(R_{asy}\) was the asymptotic recruitment level in terms of number of juveniles, and was calculated as follows:

\[\Large R_{asy} = \frac{\frac{\tilde{R}_{asy}}{SPR_{F=0}}}{(e^{-M_{Ocean}*3})}\] \(M_{ocean}\) is multiplied by 3 because the age of recruitment was defined as 3 years in Gibson’s meta-analysis (2004). \(\tilde{R}_{asy}\) was divided by the spawning biomass per recruit in the absence of fishing mortality \((SPR_{F=0})\). This value represented the rate at which Age-3 recruits produce spawners throughout their lives (Gibson 2004) and can be calculated for a specific population as follows:

\[\Large SPR_{F=0} = \sum_{a =3}^{8}SS_a W_a\] in which \(SS_a\) = spawning stock for a given age class and \(W_a\) is the average mass for each age class. Each year class contribution reflected i) probability of maturity, ii) cumulative adult mortality and iii) juvenile mortality such that:

\[\Large SS_3 = m_3\] \[\Large SS_4 = SS_3e^{-M_{adult}}+(1-m_3)e^{-M_{juv}}m_4\] \[\Large SS_8 = SS_7e^{-M_{adult}}+(1-m_3)(1-m_4)(1-m_5)(1-m_6)(1-m_7)e^{-5M_{juv}}m_8\]

in which \(M_{adult}\) = instantaneous mortality rate of mature fish, \(M_{juv}\) = instantaneous mortality rate of immature fish, and \(m_a\) = probability that immature fish alive at age a will mature at that age (Gibson 2004; Table 1). For this paper, an average SPR was used from Gibson’s meta-analysis (\(SPR_{F=0}\): 0.357 kg/recruit; 2004).

The alpha value (\(\alpha\)) is the slope of the origin for the SR curve, and was calculated similarly. Gibson (2004) provided a probability distribution for the maximum lifetime reproductive rate \((\tilde{\alpha})\), expressed in units of spawners/spawner, was first divided by SPR in the absence of fishing mortality to calculate the number of Age-3 recruits per unit spawner biomass and then by \(M_{ocean}*3\) to convert the units of recruitment to number of juvenile fish:

\[\Large \alpha = \frac{\frac{\tilde{\alpha}}{SPR_F=0}}{(e^{-M_{Ocean*3}})}\] A second standardization was required to change from units of spawner biomass to the number of eggs. \(SSB_t\) was calculated for each year of the model based on the number of fish in the spawning population as follows:

\[\Large SSB_t = \sum_{a = 3}^{8}\sum_{p=1}^{5}S_{a ,t,p}W_a\] in which for each year t, \(S_{a,t,p}\) is the number of fish of age a that have spawned p times and \(W_a\) is the mass at age a (Gibson and Myers 2003b). The model was run with a habitat size of 4.047x106 \(m^2\) for 300 years to allow the population to stabilize, and outputs then were used to iteratively estimate the \(\alpha\) value that described the slope at the origin of the SR curve that related egg and juvenile production.

Survival and mortality rates are prorated once fish mature and join the spawning population. It is assumed that spawners spend nine months in the ocean and three months in freshwater. (put in equations for instantaneous rates and explain them)

Back to Top

Dam Passage and Associated Variables


Spawner Distribution

Distribution is calculated based on habitat size, with the assumption that in the absence of any effects of dam passage, fish migrate to where the majority of spawning habitat is located. Habitat size is calculated as the acreage of lake and pond habitat between each set of dams, as well as what is available above the dam furthest upstream in the system. For example, the St. Croix River has four dams on the main stem. There are three habitat units between each pair of dams, and a fourth habitat unit that is located upstream of dam 4 (Figure 10). The area of these habitat units is shown in Table 2.

Figure 10. Habitat units and dams on the St. Croix River

Figure 10. Habitat units and dams on the St. Croix River


Habitat availability was used as a proxy for motivation for spawners to move upstream. As spawners moved upstream, the number that stayed in each habitat unit was calculated using the proportion of habitat available within that particular unit out of the total amount of upstream habitat left to spawners.

For example, on the St. Croix River a total of 60,847 acres was available as spawning habitat. The first habitat unit was 252 acres, which was 0.4% of the total acreage available in the river. That meant that 0.4% of the total spawners stayed in Habitat Unit 1 to spawn, and 99.6% of spawners continued upstream to the next habitat unit. From the second habitat unit, there was 60,595 acres of upstream spawning habitat avilable, and 1.9% of it was located in the second habitat unit. This meant that of the spawners that entered Habitat Unit 2, 98.1% continued upstream to the third habitat unit. In the St. Croix River, most of the available spawning habitat is in the last two habitat units. For spawners that entered the third habitat unit, 38.1% of fish stayed within that unit to spawn and 60.9% moved upstream to the fourth habitat unit. As the Habitat Unit 4 was the last one in the river, all remaining spawners stayed within this unit.


Table 2. Area of habitat units in the St. Croix River with percent of total river-wide spawning habitat represented by each unit.

Habitat Unit Acres Percentage of Total Upstream Habitat Available % of Spawners moving upstream from HU
1. Milltown 252 0.4% 60,847 99.6%
2. Woodland 1174 1.9% 60,595 98.1%
3. Grand Falls 23,212 38.1% 59,421 60.9%
4. Spednic 36,209 59.5% 36,209 0%
Total 60,847


Dam Passage

There were three types of passage rates associated with each dam. These included:

  • Adult upstream passage
  • Adult downstream passage
  • Juvenile downstream passage

Passage rates were defined as the percentage of fish that successfully passed a dam out of those the total number that approached the dam. These rate ranged between 0% (where no fish passed the dam) to 100% (where all fish passed the dam).


Adult Upstream Passage

Adult upstream passage rates were applied at each habitat unit as spawners moved up the river. This ensured that a cumulative effect was taken into account when multiple dams were present. Multiple variables were calculated separately for each habitat unit, and the abundance of fish moving to the next habitat unit was based on the previous unit.

For each habitat unit, the following were calculated:

  1. The number of fish that approached the dam associated with that habitat unit
  2. The number of fish that successfully pass the dam based on the passage rate applied. This is a percentage of the fish that approached the dam.
  3. The number of fish that stay in that habitat unit to spawn. This is based on the proportion of habitat available out of the total spawning habitat found in the river.
  4. The number of fish that continue upstream to approach the next dam. Like #3 above, this is also based on the proportion of habitat available within a particular habitat unit.
  5. The number of fish that successully pass the next dam in the river based on the applied passage rate for that dam.
  6. The number of fish that did not successfully pass the next dam in the river. These fish are “moved” back downstream and their total added to the number of fish that stayed to spawn in the preceding habitat unit.


Figure 11. Variables calculated for each dam and its associated habitat unit

Figure 11. Variables calculated for each dam and its associated habitat unit


Adult and Juvenile Downstream Passage

The calculations for downstream passage were less complicated as we did not have to calculate the number of fish staying within a habitat unit. We did, however, still need to calculate a cumulative effect based on where a fish started migrating downstream. This was done by taking the total number of fish leaving each habitat unit and multiplying by the passage rates for each dam those fish move past. The following table illustrates how these calculations were completed, using an example with four dams. For this example, we are assuming that each of the four dams has a passage rate of 0.9 or 90% successful passage.


Table 3. An example of calculations for downstream passage. Each dam has a passage rate of 0.9 (\(Dam_n\))

Habitat Unit Number of Fish Leaving HU Passage Equation Cumulative Passage Rate Number of Fish Successfully Entering Ocean
4 1000 \(Dam_4*Dam_3*Dam_2*Dam_1\) 0.66 660
3 1000 \(Dam_3*Dam_2*Dam_1\) 0.73 730
2 1000 \(Dam_2*Dam_1\) 0.81 810
1 1000 \(Dam_1\) 0.9 900


Passage Calculations for the St. Croix River

The History of Dam Construction on the St. Croix

Though Europeans started building small mill dams on the St. Croix River in the late 1700s, the first to span the entire river was built at the head of tide in 1836 (Fletcher & Meister 1982). This run of the river dam, located at what is now Milltown (Figure 1.1), had no fishway, and anadromous fish runs in the St. Croix River sharply declined after it was built (Atkins 1889; Decker 1967; Fletcher & Meister 1982). Although a fishway was installed in 1960 (pool and weir, height = 7.3 m), it was thought to be ineffective. Effective passage was provided in 1981 when an improved pool-and-weir fishway with resting pools was constructed (Flagg 2007), revising the over century long decline of the historically large anadromous fish runs in the St. Croix River.

Three other dams were built on the main stem of this river and are currently owned by the United States company Woodland Pulp, LLC. Woodland Dam was originally built in 1905, Grand Falls Dam in 1912, and Vanceboro Dam in 1836 (Flagg 2007). Neither Woodland nor Grand Falls dams were built with a fishway. It wasn’t until 1964 that denil fishways were constructed at both Woodland Dam (227 m) and Grand Falls Dam (183 m) (Decker 1967). When the fishway at Woodland Dam was constructed, it was the longest denil in the eastern United States (Decker 1967). Vanceboro Dam was constructed at the outlet of Spednic Lake, and it raised the water level of this natural lake. The fishway at Vanceboro Dam is a vertical slot, which likely has the highest passage rate of the four main-stem dams (Decker 1967).

Estimating “Baseline” Passage Rates

Baseline upstream passage efficiency was estimated based on escapement information from 1984-1988 recorded at the first three dams (St. Croix Milltown Trap). While these calculations were useful, the majority of the management questions that can be addressed by this model will focus on the dynamics associated with changes in passage rather than baseline values because data were limited for this river. Escapement values were used to calculate an average proportion of fish within each habitat unit (Table 5.1). This proportion was then used to calculate baseline passage probability for each dam (Table).

Passage at dam 1 was initially set at 1.0 to calculate escapement probabilities, but this is not a realistic rate. However, we have no information on the number of fish that approach this dam and therefore cannot calculate passage efficiency based on collected data. Because of this, the passage efficiency at Milltown was estimated based on the fishway type (Noonan et al. 2012). In addition, escapement was not recorded at dam 4, and so the proportion of fish was estimated based on the slope of the relationship for the previous three dams. This proportion was used to determine passage probability (Table).


Table 4. Recorded escapement levels from 1984-1988 and associated proportions entering Habitat Units 2 and 3.

Year Escapement Dam 1 Escapement Dam 2 Prop Entering HU2 Escapement Dam 3 Prop Entering HU3
1984 152,900 78,000 0.51 65,000 0.425
1985 368,900 93,000 0.252 87,000 0.236
1986 1,984,720 1,300,000 0.655 625,000 0.315
1987 2624,700 930,000 0.354 800,000 0.305
1988 2,590,750 1,004,200 0.388
Average 0.432 0.32


Table 5. Calculated baseline upstream passage probability for Dams 1-4 in the St. Croix River based on the average proportion of fish entering Habitat Units as well as best estimates for dams with no historic data available.

Dam Average Proportion Entering HU Equation Baseline Passage Probability
1. Milltown (1) changed to 0.6
2. Woodland 0.4 =0.4/1 0.4
3. Grand Falls 0.3 =0.3/(0.4x0.1) 0.75
4. Spednic 0.2 =0.2/(0.75x0.4x1) 0.67


Back to Top

R Code

Input Variables

There are several variables that need to be defined before running the model. These variables relate to population demographics as well as simulation-specific inputs. While many of these variables are kept constant between model simulations, it is possible to change their values.


These variables represent alewife population demographics that are generally kept constant for each simulation

Each variable is listed, followed by its label within the R code

  1. For each age class:
    • The probability of maturing (ProbMature)
    • The probability of spawning (ProbSpawn)
    • Mass-at-age (MassAtAge)
    • Used to calculate the distribution of fish for the first year the model is run: -The initial abundance within each age class of fish both in the ocean and within the spawning population (InitialOceanDist, InitialSpawnDist)
      • Proportion of spawners to total alewives within the population (PropSpawners)
  2. To calculate mortality rates:
    • Instantaneous annual mortality rate of fish in the ocean (InstOcean)
    • Instantaneous annual mortality rate of fish in the spawning population (InstSpawn)
    • The residence time of mature fish in freshwater in proportion of months in the year (SpawnFWRes)
    • The residence time of mature fish in the ocean in proportion of months in the year (SpawnOceanRes)
    • The residence time of juveniles in proportion of months in the year (JuvFWRes)
  3. To calculate spawner-recruit relationship:
    • alpha (alpha)
    • Recruits per acre, originally calculated as Age-3 recruits per acre (RPA)
    • YOY per acre, back-calculated from recruits per acre and juveniles mortality rates (YOYPA)
  4. Other population demographics:
    • Slope of the fecundity relationship (FecSlope)
    • y-intercept of the fecundity relationship (FecIntercept)
    • Sex ratio in males:females (MFRatio)

Actual values used in the model are listed below:

#For each age class:
ProbMature <- c(0.35, 0.51, 0.955, 1.0, 1.0, 1.0)
ProbSpawn <- rep(0.95, times = 6)
MassAtAge <- c(144, 186, 209, 244, 277, 353)
InitialOceanDist <- c(0.60, 0.24, 0.093, 0.0238, 0.0076, 0.000423, 0.000148, 0.0000535, 0.0000193)
InitialSpawnDist <- c(0.33877, 0.31471, 0.22514, 0.08143, 0.02937, 0.01059)
PropSpawners <- 0.0378

#Mortality rates:
InstOcean <- 0.648365
InstSpawn <- 2.3913
JuvFWRes <- 0.66667
SpawnFWRes <- 0.25
SpawnOceanRes <- 0.75

#Spawner-recruit:
alpha <- 0.0019
RPA <- 582.659
YOYPA <- floor(RPA / ((1-0.35095) * (1-0.4771) * (1-0.4771)))

#Other population demographics:
FecSlope <- 871.72
FecIntercept <- 50916
MFRatio <- 0.50


These variables represent simulation-specific variables that can be changed to answer a particular question

Values listed are examples from the St. Croix River

  1. Dam passage
    • The number of dams present in the river (ndams)
    • A list of the upstream passage rates for adults at each dam in the river (UpPassAdult)
    • A list of the downstream passage rates for adults at each dam in the river (DnPassAdult)
    • A list of the downstream passage rates for juveniles at each dam in the river
    • Initial number of spawners entering the river in the first year of the simulation (Escapement)
  2. Other variables
    • Size in acres for each available spawning habitat unit within the river (HabitatSize)
      • This is used to calculate the distribution of spawners in the river, which is defined according to habitat availability (Dist_US)
    • Number of years the simulation will run (nYears)
#Dam passage:
ndams<-4
UpPassAdult <- c(0.6, 0.4, 0.75, 0.67)
DnPassAdult <- c(0.9, 0.9, 0.9, 0.9)
DnPassJuv <- c(0.9, 0.9, 0.9, 0.9)
Escapement <- 1000

#Other variables:
HabitatSize <- c(252, 1174, 23212, 36209)
nYears = 100

#Distribution calculations using habitat size:
HSRev <- rev(HabitatSize)
Totals <- cumsum(HSRev)
TT <- rev(Totals)
Dist_US <- (1-(HabitatSize / TT))

Back to Top

Running simulations

This population model was written to be reactive, meaning that for a particular simulation the structure of the output will change depending on the level of variables used as input. This type of reactive model takes a lot of time and effort up front, but once it is finished it is quite simple to change the value of input variables when running a new simulation without having to change anything in the code.

For example, we can define the number of dams we want in a particular simulation by specifying ndams. If ndams is defined as 4, there are four separate output lists automatically created that are associated with each of the dams and populated with values. If we change ndams to 5 in the next simulation, 5 separate output lists are automatically created. If we were using a non-reactive version, each time a simulation was run with a different number of dams we would have to go in manually and tell R how many output tables to create.

The model was written as a function where each of these inputs could be defined at the beginning of a simulation, as shown below. This allows us to easily compare the outputs from multiple simulations.

Full_Model <- function (ProbMature, ProbSpawn, alpha, RPA, FecSlope, FecIntercept, MFRatio, ndams, 
                        UpPassAdult, DnPassAdult, DnPassJuv, HabitatSize, InstOcean, InstSpawn, 
                        JuvFWRes, SpawnFWRes, SpawnOceanRes, MassAtAge,nYears, Escapement, 
                        InitialOceanDist, PropSpawners, InitialSpawnDist, Dist_US){


1. Initial Calculations

Calculations were conducted witin the function for several inputs that may have changed between simulations. These included cumulative downstream passage rates (upstream rates are calculated later in the function), annual ocean mortality rates, the annual survival rate of spawners, and the initial population distribution.

#Passage rates
DnSurvAdultHU <- vector(mode = 'numeric', length = ndams)
DSYOYSurvHU <- vector(mode = 'numeric', length = ndams)
DSYOYSurvHU <- cumprod(DnPassJuv)
DnSurvAdultHU <- cumprod(DnPassAdult)
  
#Ocean Mortality
IntervalOceanMort <- (1-exp(-InstOcean))
YOYOceanMort <- (1-exp(-InstOcean*JuvFWRes))
# 
# #Spawning Survival
IntervalSpawnFWSurv <- (exp(-InstSpawn * SpawnFWRes))
IntervalSpawnOceanSurv <- (exp(-InstOcean * SpawnOceanRes))
# 
# #Initial population distribution
Initial <- Escapement / PropSpawners


2. Create empty lists and dataframes

Empty lists and dataframes were created to store all of the data from a simulation. Annual outputs for the ocean portion and the spawning portion of the population were kept in separate dataframes. A column was created for each age class, and a row for each year the simulation was run. Keep in mind that for each simulation, year 1 is popualted with the calculated initial values, and year 2 is when the forward-projection starts.

#Define number of rows to include:
  Years <- c(1:nYears)

##Ocean Population##
  #Specify columns for each age class:
  OceanAge <- c(1:9)
  
  #Create empty list of lists to store data:
  list_of_OceanAbund <- list()
  list_of_OceanMort <- list()
  list_of_OceanSurv <- list()
  for (k in (OceanAge)) {
    list_of_OceanAbund[[k]] <- vector(mode = 'numeric', length = nYears)
    list_of_OceanMort[[k]] <- vector(mode = 'numeric', length = nYears)
    list_of_OceanSurv[[k]] <- vector(mode = 'numeric', length = nYears)
  }
  
  #Label columns:
  Ocean_Abund <- data.frame(Years, list_of_OceanAbund)
  colnames(Ocean_Abund) <- c("Year", "OceanAbundYOY", "OceanAbundYr1", "OceanAbundYr2", "OceanAbundYr3", 
                             "OceanAbundYr4", "OceanAbundYr5", "OceanAbundYr6", "OceanAbundYr7", 
                             "OceanAbundYr8")
  Ocean_Mort <- data.frame(Years, list_of_OceanMort)
  colnames(Ocean_Mort) <- c("Year", "OceanMortYOY", "OceanMortYr1", "OceanMortYr2", "OceanMortYr3", 
                            "OceanMortYr4", "OceanMortYr5", "OceanMortYr6", "OceanMortYr7", "OceanMortYr8")
  Ocean_Surv <- data.frame(Years, list_of_OceanSurv)
  colnames(Ocean_Surv) <- c("Year", "OceanSurvYOY", "OceanSurvYr1", "OceanSurvYr2", "OceanSurvYr3", 
                            "OceanSurvYr4", "OceanSurvYr5", "OceanSurvYr6", "OceanSurvYr7", "OceanSurvYr8")


##Spawning Population##
  #Specify columns for each age class: 
  SpawnAge <- c(1:6)
  
  #Create empty lists:
  list_of_Abund <- list()
  list_of_mort <- list()
  list_of_surv <- list()
  list_spn_ocean <- list()
  for (k in SpawnAge){
    list_of_Abund[[k]] <- vector(mode = 'numeric', length = nYears)
    list_of_mort[[k]] <- vector(mode = 'numeric', length = nYears)
    list_of_surv[[k]] <- vector(mode = 'numeric', length = nYears)
    list_spn_ocean[[k]] <- vector(mode = 'numeric', length = nYears)
  }
  
  #Label columns:
  SpawnAbund <- data.frame(Years, list_of_Abund)
  colnames(SpawnAbund) <- c("Year", "SpawnYr3", "SpawnYr4", "SpawnYr5", "SpawnYr6", "SpawnYr7", "SpawnYr8")
  
  SpawnMort <- data.frame(Years, list_of_mort)
  colnames(SpawnMort) <- c("Year", "SpawnMortYr3", "SpawnMortYr4", "SpawnMortYr5", "SpawnMortYr6", 
                           "SpawnMortYr7", "SpawnMortYr8")
  
  SpawnSurv <- data.frame(Years, list_of_surv)
  colnames(SpawnSurv) <- c("Year", "SpawnSurvYr3", "SpawnSurvYr4", "SpawnSurvYr5", "SpawnSurvYr6", 
                           "SpawnSurvYr7", "SpawnSurvYr8")
  
  Spawn_Ocean <- data.frame(Years, list_spn_ocean)
  colnames(Spawn_Ocean) <- c("Year", "Spawn3_Ocean4", "Spawn4_Ocean5", "Spawn5_Ocean6", "Spawn6_Ocean7", 
                             "Spawn7_Ocean8", "Total Death")


3. Create empty lists to hold habitat-specific data

Empty lists were created to hold the output associated with each dam in the simulation, and all variables were calculated for each age class. These variables were calculated according to the distribution rules outlined above and for each habitat unit included the number of spawners that: + Entered a habitat unit (based on passage rates) + Stayed in a habitat unit to spawn + Attempted to move upstream to the next habitat unit

In addition, survival and mortality rates were calculated. The number of eggs produced for each age group of spawners was also calculated to later estimate juvenile production.

#Make list to hold the next set of lists
results <- list()
#habitat-specific lists: create a list for each dam
for (b in 1:ndams){
list_of_EnteringHU <- list()
list_of_StayHU <- list()
list_of_MovingToNextHU <- list()
list_of_SpawningMortHU <- list()
list_of_SpawningSurvHU <- list()
list_of_EggProdHU <- list()
list_of_DSMortHU <- list()
list_of_DSSurvHU <- list()
for (k in SpawnAge){
  list_of_EnteringHU[[k]] <- vector(mode = 'numeric', length = nYears)
  list_of_StayHU[[k]] <- vector(mode = 'numeric', length = nYears)
  list_of_MovingToNextHU[[k]] <- vector(mode = 'numeric', length = nYears)
  list_of_SpawningMortHU[[k]] <- vector(mode = 'numeric', length = nYears)
  list_of_SpawningSurvHU[[k]] <- vector(mode = 'numeric', length = nYears)
  list_of_EggProdHU[[k]] <- vector(mode = 'numeric', length = nYears)
  list_of_DSMortHU[[k]] <- vector(mode = 'numeric', length = nYears)
  list_of_DSSurvHU[[k]] <- vector(mode = 'numeric', length = nYears)
}
  
#create dataframes from lists and column names:

EnteringHU <- data.frame(Years, list_of_EnteringHU)
colnames(EnteringHU) <- c("Year", "EnterHUYr3", "EnterHUYr4", "EnterHUYr5", "EnterHUYr6", "EnterHUYr7", 
                          "EnterHUYr8")  

StayHU <- data.frame(Years, list_of_StayHU)
colnames(StayHU) <- c("Year", "StayHUYr3", "StayHUYr4", "StayHUYr5", "StayHUYr6", "StayHUYr7", "StayHUYr8")

MoveToNextHU <- data.frame(Years, list_of_MovingToNextHU)
colnames(MoveToNextHU) <- c("Year", "MoveToNextHUYr3", "MoveToNextHUYr4", "MoveToNextHUYr5", "MoveToNextHUYr6",
                            "MoveToNextHUYr7", "MoveToNextHUYr8")

SpawnMortHU <- data.frame(Years, list_of_SpawningMortHU)
colnames(SpawnMortHU) <- c("Year", "SpawnMortHUYr3", "SpawnMortHUYr4", "SpawnMortHUYr5", "SpawnMortHUYr6", 
                           "SpawnMortHUYr7", "SpawnMortHUYr8")

SpawnSurvHU <- data.frame(Years, list_of_SpawningSurvHU)
colnames(SpawnSurvHU) <- c("Year", "SpawnSurvHUYr3", "SpawnSurvHUYr4", "SpawnSurvHUYr5", "SpawnSurvHUYr6", 
                           "SpawnSurvHUYr7", "SpawnSurvHUYr8")

EggProdHU <- data.frame(Years, list_of_EggProdHU)
colnames(EggProdHU) <- c("Year", "EggProdHUYr3", "EggProdHUYr4", "EggProdHUYr5", "EggProdHUYr6", 
                         "EggProdHUYr7", "EggProdHUYr8")

DSMortHU <- data.frame(Years, list_of_DSMortHU)
colnames(DSMortHU) <- c("Year", "DSMortHUYr3", "DSMortHUYr4", "DSMortHUYr5", "DSMortHUYr6", "DSMortHUYr7", 
                        "DSMortHUYr8")

DSSurvHU <- data.frame(Years, list_of_DSSurvHU)
colnames(DSSurvHU) <- c("Year", "DSSurvHUYr3", "DSSurvHUYr4", "DSSurvHUYr5", "DSSurvHUYr6", "DSSurvHUYr7", 
                        "DSSurvHUYr8")
  
results[[b]]<-list(EnteringHU = EnteringHU, StayHU = StayHU, MoveToNextHU = MoveToNextHU, 
                   SpawnMortHU = SpawnMortHU, SpawnSurvHU = SpawnSurvHU, EggProdHU = EggProdHU, 
                   DSMortHU = DSMortHU, DSSurvHU = DSSurvHU)
}

#Create separate data frames for recruits (BH) and total eggs produced (TotalEggs)
temp <- ndams+1
BH <- data.frame(matrix(ncol = temp, nrow = nYears))
colnames(BH) <- c("Year",(lapply(1:ndams, function(x){paste0("BH", x)})))
BH[,1] <- c(1:nYears)

TotalEggs <- data.frame(matrix(ncol = temp, nrow = nYears))
colnames(TotalEggs) <- c("Year", (lapply(1:ndams, function(x){paste0("TotalEggs",x)})))
TotalEggs[,1] <- c(1:nYears)


4. Define initial values within the framework

Now that we can created the tables, we can start putting numbers in them. Initial values need to be calculated at the start of each simulation. These values are used to define the abundance of fish for each variable in each age class at Year 1, and were primarily determined using age-structured data from the St. Croix River. These initial values are used to calculate variables like survival and mortality for each age class as the model moves forward into Year 2. As such, these initial values are defined in the code differently than the values for Year 2 onward.

##These are the values in row 1, which are defined differently than in rows 2:nrow

#First create dataframes to hold initial values:
#Ocean Population:
Year <- c("YOY", "Yr1", "Yr2", "Yr3", "Yr4", "Yr5", "Yr6", "Yr7", "Yr8")
Initial_Ocean_Pop <- (data.frame(Age = Year, OceanAb = floor(Initial * InitialOceanDist)))

#Spawning population:
Spawn_Year <- c("SpawnYr3", "SpawnYr4", "SpawnYr5", "SpawnYr6", "SpawnYr7", "SpawnYr8")
Initial_Spawning_Pop <- data.frame(SpawnAge = Spawn_Year, SpawnAb = floor((Escapement * InitialSpawnDist)))

#Set initial parameters for each row in ocean population and spawning population (THIS IS YEAR 1), and calculate all 
#further values based on these abundances so can carry forward to next year

####Initial values for Ocean Population####
#Ocean Abundance:
Ocean_Abund[1,]$OceanAbundYOY <- Initial_Ocean_Pop[1,2]
Ocean_Abund[1,3:10] <- Initial_Ocean_Pop[2:9,2]

#Ocean Mort:
Ocean_Mort[1,]$OceanMortYOY <- floor((Ocean_Abund[1,]$OceanAbundYOY * YOYOceanMort))
Ocean_Mort[1,3:10] <- floor((Ocean_Abund[1,3:10] * IntervalOceanMort))

#Ocean Survival:
Ocean_Surv[1,]$OceanSurvYOY <- floor((Ocean_Abund[1,]$OceanAbundYOY * (1-YOYOceanMort)))
Ocean_Surv[1,]$OceanSurvYr1 <- floor((Ocean_Abund[1,]$OceanAbundYr1 * (1 - IntervalOceanMort)))
Ocean_Surv[1,4:10] <- floor((Ocean_Abund[1,4:10] * (1-IntervalOceanMort) * (1 - ProbMature)))

####Spawning Population####
#First define initial spawner abundance so can calculate habitat-specific abundances
SpawnAbund[1,2:7] <- Initial_Spawning_Pop[,2]

##Next we need to determine habitat-specific abundances##
####Habitat-Specific equations for row 1####
if (ndams == 1){
    results[[1]]$EnteringHU[1,2:7] <- floor((SpawnAbund[1,2:7] * ProbSpawn * UpPassAdult[1]))
    } else {
      #This calculates the number of fish entering HU1
      results[[1]]$EnteringHU[1,2:7] <- floor((SpawnAbund[1,2:7] * ProbSpawn * UpPassAdult[1]))
      }
for (c in 1:ndams){
    if (c!=ndams){
      #This calculates the number of fish that move upstream to attempt to pass the next dam
      results[[c]]$MoveToNextHU[1,2:7] <- floor((results[[c]]$EnteringHU[1,2:7] * Dist_US[c]))
      #This calculates the number of fish that successfully pass the dam to enter the next HU
      results[[c+1]]$EnteringHU[1,2:7] <- floor((results[[c]]$MoveToNextHU[1,2:7] * UpPassAdult[c+1])) 
      #This calculates the number of fish that stay within a give HU 
      #(based on habitat and those that can't pass the next dam)
      results[[c]]$StayHU[1,2:7] <- floor(((results[[c]]$EnteringHU[1,2:7] - 
                                              results[[c]]$MoveToNextHU[1,2:7]) + 
                                       (results[[c]]$MoveToNextHU[1,2:7] - 
                                          results[[c+1]]$EnteringHU[1,2:7])))
      } else {
        #This calculates the number of fish staying in the final HU (as fish no longer move upstream)
        results[[c]]$StayHU[1,2:7] <- floor(((results[[c]]$EnteringHU[1,2:7])))
          }

#Now define spawning mortality and spawn survival:
    results[[c]]$SpawnMortHU[1,2:7] <- floor((results[[c]]$StayHU[1,2:7] * (1-IntervalSpawnFWSurv)))
    results[[c]]$SpawnSurvHU[1,2:7] <- floor((results[[c]]$StayHU[1,2:7] * (IntervalSpawnFWSurv)))

#Define egg production:
    results[[c]]$EggProdHU[1,2:7] <- floor((((results[[c]]$StayHU[1,2:7] * IntervalSpawnFWSurv * ProbSpawn) * 
                                              MFRatio) * ((FecSlope*MassAtAge) - FecIntercept)))

#Define downstream adult mortality:
    results[[c]]$DSMortHU[1,2:7] <- floor((results[[c]]$StayHU[1,2:7] *
                                            (1-(IntervalSpawnFWSurv * DnSurvAdultHU[c]))))

#Define downstream adult survival:
    results[[c]]$DSSurvHU[1,2:7] <- floor((results[[c]]$StayHU[1,2:7] * IntervalSpawnFWSurv 
                                           * DnSurvAdultHU[c]))  

}

#Define total egg production and recruitment for first row for each Habitat unit:
temp <- ndams+1
TotalEggs[1,2:temp] <- lapply(1:ndams, function(x){rowSums(results[[x]]$EggProdHU[1,2:7])})
#Define initial row of YOY produced for each HU:
BH[1,2:temp] <- floor((alpha * TotalEggs[1,2:temp]) / 
  (1 + ((alpha * TotalEggs[1,2:temp]) / 
  (HabitatSize * YOYPA))))

####Calculate Total Spawning mortality, survival to next age class, and Spawn_Ocean for row 1####
#Spawning Mortality:
Spn <- data.frame(matrix(nrow = 6, ncol = ndams))
for (d in 1:ndams){
  Spn[d] <- unlist(results[[d]]$DSSurvHU[1,2:7])
  }
SpawnSurv[1,2:6] <- floor((rowSums(Spn[1:5,1:ndams])*IntervalSpawnOceanSurv))
SpawnSurv[1,7] <- floor(rowSums(Spn[6,1:ndams]))

#Fish that mature but do not spawn (small percentage of spawners)
Spawn_Ocean[1,2:7] <- floor((SpawnAbund[1,2:7] * IntervalSpawnFWSurv * IntervalSpawnOceanSurv * 
                               (1 - ProbSpawn)))

#Mortality related to spawning
SpawnMort[1,2:7] <- floor(SpawnAbund[1,2:7] - Spawn_Ocean[1,2:7] - SpawnSurv[1,2:7])

####Maturing Fish####
#first create dataframe for maturing to spawning population:
list_of_maturing <- list()
for (k in (SpawnAge)) {
  list_of_maturing[[k]] <- vector(mode = 'numeric', length = nYears)
  }

Maturing <- data.frame(Years, list_of_maturing)
colnames(Maturing) <- c("Year", "Mature2", "Mature3", "Mature4", "Mature5", "Mature6", "Mature7")

#Define first row for maturing fish:
Maturing[1,2:7] <- floor(Ocean_Abund[1,4:9] * ProbMature * (1 - IntervalOceanMort))


5. Calculate values for the remaining years in the simulation

Now that the initial row has been defined, calculations can be made for each remaining year in the simulation according to the formulas defined in the code. Many of the calculations for a particular year were based on the value of a particular variable in the previous year.

####Now need to define the rest of the rows (Years 2 to nYears)####
for (l in 2:nYears) {
  
  #First define recruitment based on previous year
  Recruit <- floor(data.frame(mapply(`*`, BH[,2:temp], DSYOYSurvHU)))
  Recruit$TotalRecruits <- Recruit %>% select(contains("BH")) %>% rowSums(na.rm = TRUE)
  
  #Define Ocean Population based on previous year's abundances:
  Ocean_Abund[l,]$OceanAbundYOY <- floor(Recruit$TotalRecruits[l-1])
  Ocean_Abund[l,]$OceanAbundYr1 <- Ocean_Surv[l-1,]$OceanSurvYOY
  Ocean_Abund[l,]$OceanAbundYr2 <- Ocean_Surv[l-1,]$OceanSurvYr1
  Ocean_Abund[l,]$OceanAbundYr3 <- Ocean_Surv[l-1,]$OceanSurvYr2
  Ocean_Abund[l,]$OceanAbundYr4 <- (Ocean_Surv[l-1,]$OceanSurvYr3 + Spawn_Ocean[l-1,]$Spawn3_Ocean4)
  Ocean_Abund[l,]$OceanAbundYr5 <- (Ocean_Surv[l-1,]$OceanSurvYr4 + Spawn_Ocean[l-1,]$Spawn4_Ocean5)
  Ocean_Abund[l,]$OceanAbundYr6 <- (Ocean_Surv[l-1,]$OceanSurvYr5 + Spawn_Ocean[l-1,]$Spawn5_Ocean6)
  Ocean_Abund[l,]$OceanAbundYr7 <- (Ocean_Surv[l-1,]$OceanSurvYr6 + Spawn_Ocean[l-1,]$Spawn6_Ocean7)
  Ocean_Abund[l,]$OceanAbundYr8 <- (Ocean_Surv[l-1,]$OceanSurvYr7 + Spawn_Ocean[l-1,]$Spawn7_Ocean8)

  #Ocean Mort:
  Ocean_Mort[l,]$OceanMortYOY <- floor((Ocean_Abund[l,]$OceanAbundYOY * YOYOceanMort))
  Ocean_Mort[l,3:10] <- floor((Ocean_Abund[l,3:10] * IntervalOceanMort))
 
  #Ocean Survival:
  Ocean_Surv[l,]$OceanSurvYOY <- floor((Ocean_Abund[l,]$OceanAbundYOY * (1-YOYOceanMort)))
  Ocean_Surv[l,]$OceanSurvYr1 <- floor((Ocean_Abund[l,]$OceanAbundYr1 * (1 - IntervalOceanMort)))
  Ocean_Surv[l,4:10] <- floor((Ocean_Abund[l,4:10] * (1-IntervalOceanMort) * (1 - ProbMature)))
  
  #Define SpawnAbund Based on previous year's abundances:
  SpawnAbund[l,]$SpawnYr3 <- Maturing[l-1,]$Mature2
  SpawnAbund[l,]$SpawnYr4 <- Maturing[l-1,]$Mature3 + SpawnSurv[l-1,]$SpawnSurvYr3
  SpawnAbund[l,]$SpawnYr5 <- Maturing[l-1,]$Mature4 + SpawnSurv[l-1,]$SpawnSurvYr4
  SpawnAbund[l,]$SpawnYr6 <- Maturing[l-1,]$Mature5 + SpawnSurv[l-1,]$SpawnSurvYr5
  SpawnAbund[l,]$SpawnYr7 <- Maturing[l-1,]$Mature6 + SpawnSurv[l-1,]$SpawnSurvYr6
  SpawnAbund[l,]$SpawnYr8 <- Maturing[l-1,]$Mature7 + SpawnSurv[l-1,]$SpawnSurvYr7

  #Now define habitat-specific values
    if (ndams == 1){
    results[[1]]$EnteringHU[l,2:7] <- floor((SpawnAbund[l,2:7] * ProbSpawn * UpPassAdult[1]))
      } else {
      results[[1]]$EnteringHU[l,2:7] <- floor((SpawnAbund[l,2:7] * ProbSpawn * UpPassAdult[1]))
        }
    for (c in 1:ndams){
      if (c!=ndams){
        results[[c]]$MoveToNextHU[l,2:7] <- floor((results[[c]]$EnteringHU[l,2:7] * Dist_US[c])) 
        results[[c+1]]$EnteringHU[l,2:7] <- floor((results[[c]]$MoveToNextHU[l,2:7] * UpPassAdult[c+1]))
        results[[c]]$StayHU[l,2:7] <- floor(((results[[c]]$EnteringHU[l,2:7] - 
                                                results[[c]]$MoveToNextHU[l,2:7]) + 
                                               (results[[c]]$MoveToNextHU[l,2:7] - 
                                                results[[c+1]]$EnteringHU[l,2:7])))
        } else{
          results[[c]]$StayHU[l,2:7] <- floor(((results[[c]]$EnteringHU[l,2:7])))
          }

      #Now define spawning mortality and spawning survival:
      results[[c]]$SpawnMortHU[l,2:7] <- floor((results[[c]]$StayHU[l,2:7] * (1-IntervalSpawnFWSurv)))
      results[[c]]$SpawnSurvHU[l,2:7] <- floor((results[[c]]$StayHU[l,2:7] * IntervalSpawnFWSurv))
      
      #Define egg production:
      results[[c]]$EggProdHU[l,2:7] <- floor((((results[[c]]$StayHU[l,2:7] * IntervalSpawnFWSurv * ProbSpawn) 
                                               * MFRatio) * ((FecSlope*MassAtAge) - FecIntercept)))
      
      #Define DS Adult Mortality (I think this may be define incorrectly):
      results[[c]]$DSMortHU[l,2:7] <- floor((results[[c]]$StayHU[l,2:7] * 
                                            (1-(IntervalSpawnFWSurv * DnSurvAdultHU[c]))))
      
      #Define DS Adult Survival:
      results[[c]]$DSSurvHU[l,2:7] <- floor((results[[c]]$StayHU[l,2:7] * IntervalSpawnFWSurv * 
                                            DnSurvAdultHU[c]))  
    }

  #Total Eggs and Recruits out of ndams loop b/c don't need to cycle through, 
  #but still in nYears loop because need to by calculated by row!   
  TotalEggs[l,2:temp] <- lapply(1:ndams, function(x){rowSums(results[[x]]$EggProdHU[l,2:7])})
  
  #Define initial row of YOY produced for each HU:
  BH[l,2:temp] <- floor(((alpha * TotalEggs[l,2:temp]) / 
                           (1 + ((alpha * TotalEggs[l,2:temp]) / 
                                   (HabitatSize * YOYPA)))))
   
  #Define mortality as the sum of each age class from each habitat unit:
  #First create dataframes to hold information from each age class:
  StayYr3 <- data.frame(matrix(ncol = ndams, nrow = nYears))
  StayYr4 <- data.frame(matrix(ncol = ndams, nrow = nYears))
  StayYr5 <- data.frame(matrix(ncol = ndams, nrow = nYears))
  StayYr6 <- data.frame(matrix(ncol = ndams, nrow = nYears))
  StayYr7 <- data.frame(matrix(ncol = ndams, nrow = nYears))
  StayYr8 <- data.frame(matrix(ncol = ndams, nrow = nYears))
  
  #Then select data for each age class from "StayHU":
  for (x in 1:ndams){
    unlist(results[[x]]$StayHU)
    StayYr3[,x] <- select(results[[x]]$StayHU, matches("StayHUYr3"))
    StayYr4[,x] <- select(results[[x]]$StayHU, matches("StayHUYr4"))
    StayYr5[,x] <- select(results[[x]]$StayHU, matches("StayHUYr5"))
    StayYr6[,x] <- select(results[[x]]$StayHU, matches("StayHUYr6"))
    StayYr7[,x] <- select(results[[x]]$StayHU, matches("StayHUYr7"))
    StayYr8[,x] <- select(results[[x]]$StayHU, matches("StayHUYr8"))
  }
  
  #Now calculate survival:
  #First create dataframes for each age class:
  SpawnYr3 <- data.frame(matrix(ncol = ndams, nrow = nYears))
  SpawnYr4 <- data.frame(matrix(ncol = ndams, nrow = nYears))
  SpawnYr5 <- data.frame(matrix(ncol = ndams, nrow = nYears))
  SpawnYr6 <- data.frame(matrix(ncol = ndams, nrow = nYears))
  SpawnYr7 <- data.frame(matrix(ncol = ndams, nrow = nYears))
  SpawnYr8 <- data.frame(matrix(ncol = ndams, nrow = nYears))
  
  #Then select data from each age class from "DSSurvHU":
  for (x in 1:ndams){
    unlist(results[[x]]$DSSurvHU)
    SpawnYr3[,x] <- select(results[[x]]$DSSurvHU, matches("DSSurvHUYr3"))
    SpawnYr4[,x] <- select(results[[x]]$DSSurvHU, matches("DSSurvHUYr4"))
    SpawnYr5[,x] <- select(results[[x]]$DSSurvHU, matches("DSSurvHUYr5"))
    SpawnYr6[,x] <- select(results[[x]]$DSSurvHU, matches("DSSurvHUYr6"))
    SpawnYr7[,x] <- select(results[[x]]$DSSurvHU, matches("DSSurvHUYr7"))
    SpawnYr8[,x] <- select(results[[x]]$DSSurvHU, matches("DSSurvHUYr8"))
  }
  
  #Then calculate Survival:
  SpawnSurv[l,2] <- floor((rowSums(SpawnYr3[l,])*IntervalSpawnOceanSurv))
  SpawnSurv[l,3] <- floor((rowSums(SpawnYr4[l,])*IntervalSpawnOceanSurv))
  SpawnSurv[l,4] <- floor((rowSums(SpawnYr5[l,])*IntervalSpawnOceanSurv))
  SpawnSurv[l,5] <- floor((rowSums(SpawnYr6[l,])*IntervalSpawnOceanSurv))
  SpawnSurv[l,6] <- floor((rowSums(SpawnYr7[l,])*IntervalSpawnOceanSurv))
  SpawnSurv[l,7] <- floor((rowSums(SpawnYr8[l,])*IntervalSpawnOceanSurv))
  
  #Calculate the number of fish that don't spawn but return to the ocean
  Spawn_Ocean[l,2:7] <- floor((SpawnAbund[l,2:7] * IntervalSpawnFWSurv * 
                              IntervalSpawnOceanSurv * (1 - ProbSpawn)))

  #Calculate spawning mortality
  SpawnMort[l,2:7] <- floor(SpawnAbund[l,2:7] - Spawn_Ocean[l, 2:7] - SpawnSurv[l, 2:7])
    
  #Maturing Fish:
  Maturing[l,]$Mature2 <- floor((Ocean_Abund[l,]$OceanAbundYr2 * ProbMature[1] * (1 - IntervalOceanMort)))
  Maturing[l,]$Mature3 <- floor((Ocean_Abund[l,]$OceanAbundYr3 * ProbMature[2] * (1 - IntervalOceanMort)))
  Maturing[l,]$Mature4 <- floor((Ocean_Abund[l,]$OceanAbundYr4 * ProbMature[3] * (1 - IntervalOceanMort)))
  Maturing[l,]$Mature5 <- floor((Ocean_Abund[l,]$OceanAbundYr5 * ProbMature[4] * (1 - IntervalOceanMort)))
  Maturing[l,]$Mature6 <- floor((Ocean_Abund[l,]$OceanAbundYr6 * ProbMature[5] * (1 - IntervalOceanMort)))
  Maturing[l,]$Mature7 <- floor((Ocean_Abund[l,]$OceanAbundYr7 * ProbMature[6] * (1 - IntervalOceanMort)))

}


6. Define outputs

The outputs are define as a list that can be called and used in the application to create figures and tables.

List_of_Outputs <- list(Recruitment = BH, YOYSurv = Recruit, Maturing = Maturing, Ocean_Abund = Ocean_Abund, 
                        Ocean_Mort = Ocean_Mort, Ocean_Surv = Ocean_Surv, Spawn_Ocean = Spawn_Ocean, 
                        SpawnAbund = SpawnAbund, SpawnMort = SpawnMort, SpawnSurv = SpawnSurv,HUs=results)

return(List_of_Outputs)
}


7. Call the function

Now that all of the equations in the model are defined within a function, we can run a simulation by calling the function. We define the name of our simulation (here I called it “Test1”), then call our function (which is names “Full_Model”) with specific inputs defined. To change the value of these inputs, we can define new values for the variables listed within the parentheses.

Test1 <- Full_Model(ProbMature, ProbSpawn, alpha, RPA, FecSlope, FecIntercept, MFRatio, ndams, UpPassAdult, 
                    DnPassAdult, DnPassJuv, HabitatSize, InstOcean, InstSpawn, JuvFWRes, SpawnFWRes, 
                    SpawnOceanRes, MassAtAge, nYears, Escapement, InitialOceanDist, PropSpawners, 
                    InitialSpawnDist, Dist_US)

Back to Top

Application

The online application was created in the R package shiny. This application was designed to allow the user to change the value of inputs for the alewife population model, run a simulation, and receive output in the form of figures and tables.

#These are the packages needed to run the application from RStudio:
library(shiny)
library(shinyWidgets)
library(shinycssloaders)

#This is the link to the application that is stored on GitHub:
runGitHub("alewifepopmodel","moctezumaii")


1. The first page of the application contains links to background information. The YouTube tutorial is a step-by-step example of how to use the application. The PowerPoint presentation is tailored for a general audience and has background information about how the alewife population model was set up and how the online application can be used.

2. When the application is started, the user can either choose to use the pre-set values for the St. Croix River or “create”" their own river.


3. If the user chooses to create their own river, the application can accomodate up to 10 dams. The user can specify the name of each dam and the acreage of available spawning habitat between each pair of dams and above the furthest dam upstream. Once these are defined, the user can click “Load” to move to the next page.


4. On the next page, each of the dams will be displayed. Associated with each dam will be upstream and downstream passage rates for spawning adults, and downstream passage rates for juveniles. The value for passage can be changed using the slider. The number of years that the simulation is ran can also be defined. The user simply has to define the values they want to explore and click “Run”. When using the preset values for the St. Croix River, clicking on “Reset Values” sets all passage parameters back to estimated “baseline” rates.


5. Once the user clicks “Run”, the simulation will run for a few seconds. Results take longer to load when the model runs for more years. Results are displayed as a figure showing the theoretical spawner abundance for the habitat unit associated with each dam. In addition, the values calculated in the simulation are listed in tables that are displayed below the figure. These results are divided into outputs associated with juvenile abundance, those associated with the ocean population, and those associated with the spawning population. These tables can also be downloaded using the buttons located above the displayed table.



6. For each simulation, the total number of spawners can be stored and displayed in a comparative graph and table. The table lists the passage rates associated with each simulation being compared, as well as the total theoretical abundance of spawners. The figure is a barplot that shows the total theoretical spawners abundance associated with each simulation being compared. To start this feature, the user needs to click the button that says “Start or reset”. After running a simulation, the user can click the button that says “Save total for comparison”, and the table and figure are automatically updated.

To make correct comparisons, simulations must use the same number of years!



Back to Top

Example Question and interpretation:


Question: How might removing Milltown dam in the St. Croix River affect the recovery of the alewife population?

For this question, we will set passage at Milltown Dam (Dam #1) to 100% for adult upstream, adult downstream, and juvenile downstream rates. This means that 100% of fish are able to pass this dam, which is representative of a dam removal. This is the only change we will make. Passage rates at the other three dams will be held at a “baseline” value.

Inputs

  1. Milltown Dam:
  • Adult Upstream: 60% vs. 100%
  • Adult Downstream: 90% vs. 100%
  • Juvenile Downstream: 90% vs. 100%
  1. Woodland Dam:
  • Adult Upstream: 40%
  • Adult Downstream: 90%
  • Juvenile Downstream: 90%
  1. Grand Falls:
  • Adult Upstream: 75%
  • Adult Downstream: 90%
  • Juvenile Downstream: 90%
  1. Spednic
  • Adult Upstream: 67%
  • Adult Downstream: 90%
  • Juvenile Downstream: 90%
  1. Years: 100
  2. Initial Number of Spawners: 1000 fish


To address this question, we will first simulate the alewife population in the St. Croix using the “baseline” passage rates. Then we will change passage rates at Milltown dam to 100% for all three categories. Finally, we will compare the results from the two scenarios.


Here are the results from the first simulation using the baseline passage rates at each dam for the St. Croix River. Theoretical spawner abundance is the number of fish entering the river, which is before spawning mortality is applied.


Simualtion 1: Baseline passage rates


And here are the results from the second simulation where passage at Milltown Dam is set to 100% for all three rates.


Simulation 2: Milltown Dam Removed


We can look at the total theoretical spawner abundance associated with each simulation. As you can see in the following figure, when Milltown is removed the total number of spawners entering the river is much higher.



Interpretation

When the passage rates at Milltown Dam are set to 100%, the total theoretical number of spawners entering the river and entering each habitat unit increases. However, the river-wide distribution pattern for spawner abundance in each habitat unit remains the same because passage at Woodland Dam is still quite low and most of the spawning habitat is located upstream. The majority of spawners are still found within the first habitat unit despite its small size because of low passage rates at Woodland Dam. Removing the dam at the head of tide helps reduce the cumulative effect of both upstream and downstream passage for fish migrating throughout the river. This means that regardless of where fish distribute within the river, a larger number of both spawners and juveniles are surviving their downstream migration. In the simulation, this affects the age distribution of spawners and increases the number of older fish within the population. This in turn increases population fecundity, as older fish generally produce more eggs because of their larger body size. All of these effects combined lead to a much larger total spawner abundance in the simulation where Milltown is removed than in the simulation where Milltown is present and is passage is set at the estimated baseline rates.



Back to Top