Targeting river operations to the critical thermal window of fish incubation: Model and case study on Sacramento River winter‐run Chinook salmon

Allocating reservoir flows to meet societal and ecosystem needs under increasing water demands and climatic variability presents challenges to resource managers. Often, rivers have been regulated to meet flow and temperature compliance points or mimic historical patterns. Because it is difficult to assess if this approach is efficient, process‐based models are being used to design river operations. This paper describes a model for fish incubation survival based on the premise that mortality from thermal stress occurs over a critical window (CW) of embryo development. A model for the embryo CW based on metabolic studies of development is combined with density‐dependent and background mortalities to describe salmonid survival from egg fertilization to fry detection downstream. The model is calibrated with a two‐decade dataset of Sacramento River winter‐run Chinook salmon egg‐to‐fry survival. The effects of temperature exposure over a range of CWs were explored. Based on statistical and biological support, two alternative CWs were identified for temperature control: the entire incubation period and a short duration window prior to hatching. Survival under different CW assumptions and temperature control operations were simulated with an internet‐accessible form of the model. The analysis indicated that under years of limited cold‐water resources, targeting water releases to the CW prior to hatching would yield the highest incubation survival.

McCloud and Pit rivers in the spring and summer, were displaced to the warmer Sacramento River habitat by the construction of Shasta Reservoir (1938Reservoir ( -1945 (Fisher, 1994;Williams, 2006) (Figure 1). In most years, mimicking their historical conditions has been possible with cool water releases from the reservoir during the spring and summer (NMFS, 2009). (Seager et al., 2015) and with climate warming it has become increasingly difficult to control the river temperature (Beer & Anderson, 2013;CDWR, 2021;Hallnan, Saito, Busby, & Tyler, 2020;Thompson et al., 2011). In this paper we develop a mechanistic model for the effect of temperature on salmon incubation and survival. We calibrate the model with a two-decade time series of SRWRC data and then evaluate impacts of alternative Shasta

| METHODS
At the core of the system is a mechanistic model of the effects of temperature on the development and survival of individual embryos, which is expressed in terms of a critical window (CW) of temperature.
Applying the CW framework at a habitat scale, the system predicts the effects of river temperature on population-level egg-to-fry survival. Thus, in this framework, reservoir operations shift from mimicking historical temperature patterns to operations that target coldwater releases to the transient appearance of embryo CWs across the riverine habitat. We first develop a model of a CW based on embryo metabolism and then use a simplified CW form to describe a population of embryos incubating across a river habitat.

| Embryo model
Several studies of freshwater and marine fish species assume that the CW for population viability includes all phases of embryogenesis, that is, fertilization to fry emergence (Dahlke, Wohlrab, Butzin, & Pörtner, 2020;Martin et al., 2017). However, physiological studies indicate that fish thermal sensitivity is most important at the end of organogenesis when internal organs and eyes are formed, and the embryo is about to hatch. Studies suggests that mortality associated with hatching involves hypoxia (Del Rio, Davis, Fangue, & Todgham, 2019;Greig, Sear, & Carling, 2007).
We focus on the duration of hypoxia as a surrogate for the level of mortality. Prior to hatching, the metabolic oxygen demand is proportional to the embryo mass while the cutaneous supply through membrane diffusion is proportional to egg surface area (Martin et al., 2020). As the mass-to-area ratio increases with age, the oxygen demand increases more rapidly than supply and at some point, demand can exceed supply resulting in hypoxia. This transition likely occurs over a few days because the oxygen demand increases exponentially while the supply increases linearly (Rombough, 1994). Upon hatching, branchial respiration augments the diffusive supply (Wells & Pinder, 1996) and the embryo can recover quickly from hypoxia (Rombough, 2007) even though the surrounding oxygen level remains unchanged. Also, because the gill surface area increases as the square of body mass, oxygen delivery can increase faster than demand after hatching. Additionally, alevin movement can enhance redd ventilation and decrease distance to the gravel-water interface (Dill & Northcote, 1970). After hatching oxygen sensitivity is largely independent of temperature in salmonids (Rombough, 1986;Wood, Clark, Elliott, Frappell, & Andrewartha, 2019). Furthermore, SRWRC alevin are tolerant of temperatures up to 16.8 C, which is generally higher, has been observed historically in salmonid habitats (Myrick & Cech, 2004;USFWS, 1999). Thus, metabolic information suggests a threshold-like onset and cessation of hypoxia is likely associated with hatching. We designate this to the hypoxia CW.
Because the respiration rate and hatch age are both functions of temperature, the CW duration can be modeled. First, the salmon embryo oxygen demand rate R (μgO 2 /hr) is logR ¼ a R þ b R logy þ c R logT where y is age (days post fertilization), T is temperature ( C), and the coefficients obtained from laboratory studies are a R = À16.897, b R = 2.873, and c R = 3.840 (Rombough, 1994). Rearranging the equation, the CW onset, which is the age at which R first exceeds the oxygen supply rate F, is Next, the hatching age, designating the end of the CW, can be approx- Þwhere ATU H = 515 C Á d is the accumulated temperature units (ATU) ( C Á d) required for hatching and c H is an empirical adjustment factor (Alderdice & Velsen, 1978), which here is set to zero. The CW of hypoxia duration is then This equation does not predict the mortality rate, only days of hypoxia. Nor does it imply that the temperature T fully defines the window duration. Studies indicate that temperature stress weeks before hatching, in the gastrulation phase (Irvine, 2020;Mueller et al., 2015) could alter the efficiency of respiration and cellular development, which in turn could affect the coefficients in Equation (1).

| Habitat scale egg-to-fry model
At the population level, survival between egg fertilization to fry detection downstream is described with three processes: temperature-dependent mortality (TDM) occurring during egg incubation; density-dependent mortality (DDM) characterizing the competition between spawners for suitable habitat; and background mortality (BGM) mainly characterizing predation on fry after emergence and prior to detection downstream (Table 1).
To describe TDM, first denote β y as the daily mortality rate associated with temperature where y is age and fry emerge at y = E. Thermal survival is then V ¼ exp P E y¼0 β y . The daily mortality rate changes with both temperature and developmental phase and is unknown. However, the rate over a specific CW can be expressed using a modification of a model (Martin et al., 2017) in which the thermal mortality rate depends on the differential between the daily temperature T y and a critical temperature T crit as Δ y ¼ max T y À T crit ,0 ð Þwhere temperature at age y is specific to redd i and depends on redd location and date of fertilization d such that the daily temperature for redd i is T y ¼ T diþy .
The CW is now characterized by a square wave that is zero everywhere except within the interval Y -δ to Y where Y is the end date of the CW. The mortality rate coefficient within the window is b δ and thermal survival in redd i is then Thus, thermal mortality only occurs within the CW and is defined by the embryo age Y i at the CW end, the window duration δ, the temperature over the window T y,i , the critical temperature T crit , and the rate coefficient b δ . Because b δ is constant and characterizes the mortality rate per unit of Δ y , the window width and b δ are then related as where α is the intrinsic mortality rate coefficient that, as we argue in Appendix S1B, is proportional to the oxygen flux to the embryo prior to hatching.
At the end of the CW the embryo is no longer susceptible to thermal mortality. In the Embryo model, this transition is defined by the accumulated temperature for hatching, that is, ATU H . In the Egg-to-Fry model, the ATU CW endpoint is a free fitted population-level parameter, ATU Y . Consequently, the embryo age at the end of the window is different for each redd i as set by The ATU at the beginning and middle of the CW are defined where T δ is the mean temperature over the CW.
The density-dependent mortality follows the Beverton-Holt model (Beverton & Holt, 2012), where the effect of population density on mortality depends on a habitat carrying capacity D expressing the density of spawners at which survival is half the maximum level as where ρ i designates the density of redds surrounding redd i and is calculated for each river reach segment by dividing the number of redds in a reach by its length. Correspondingly, D has units of redds/km.
Note the mortality determined from Equation (6) is not a function of redd deposition date, temperature or flow.
Any additional mortality not associated with temperature or spawner density is designated background mortality and the corresponding background survival, B is a free fitted parameter that is assumed to be constant for all years and redd locations.
The total population level incubation survival in year j is the product of the three survivals as where A j is the total number of redds produced in year j.

| Calibration strategy
The model parameters [δ, b δ, ATU Y , T crit , B, D] were estimated by fixing δ and estimating the free parameter set P = [b δ, ATU Y , T crit , B, D] by minimizing the logit objective function T A B L E 1 Parameter definitions for Embryo and Egg-to-Fry models  where S obs,j is survival in year j and S j is estimated from Equation (7) as implemented in the Egg-to-Fry model available at ) (Appendix S1A). S j was estimated with observed daily temperature, T i,j , and spawning date, d i,j . P was fit using the derivative-free optimization Nelder-Mead simplex nmkb in the dfoptim package of the R© programming language (R Core Team, 2018) where parameter search limits were set.
In setting search limits, our goal was not to find a global fit but to identify parameter sets for biologically plausible CWs, specified by duration δ and endpoint ATU Y . To this end, we first estimated P across the possible range of δ (1-79 days) using parameter limits: 0.01 < b δ < 2, 0.4 < B < 0.6, 50 < D < 100, 11 < T crit < 13, and 100 < ATU Y < 947. Initial search values were taken as midpoints of ranges except for b δ which from Equation (3) was set as b δ = 2/δ. The fitted ATU Y represents a CW endpoint for the associated δ, The window start ATU 0 and middle ATU M were then approximated with Equation (5)

| Model application
We use the model to evaluate the effects of alternative Shasta Reservoir operations on the survival of SRWRC from egg fertilization to fry passage at RBDD (Figure 1). In the spring and summer, the reservoir operation goal is to maintain cool temperatures during embryo incubation, but the control ability depends on the available volume of cold water in the reservoir prior to fish spawning.
Consequently, the US Bureau of Reclamation plans the operations based on the cold-water pool and total reservoir volume expected on May 1. Other factors considered include the flows required to avoid dewatering the redds and meet water delivery obligations for downstream users. Additionally, the assessment uses long-term weather forecasts and current reservoir levels to assess the probability of refilling the reservoir the following year. Using this information, the planned operation is set into one of four Tiers (Table 2) which are designed to best allocate the available cool water to a

| Data
The model was calibrated with SRWRC salmon spawning dates and locations from carcass surveys conducted by the California Cooperative Anadromous Fish and Habitat Data Program (Killam, 2021). Survival to RBDD was determined by expanding fry passage estimates at RBDD divided by estimates of female spawner carcass counts (Poytress, 2016;Rea, 2019Rea, , 2021   For a fixed oxygen flux of 25 μgO 2 /hr a change in temperature from 11 to 13 C increases the hypoxia duration by 2 days, while a drop in oxygen flux from 25 to 20 μgO 2 /hr over the same temperature range increases the duration by 3-4 days. At lower fluxes, temperature has no effect on hypoxia duration.

| Egg-to-fry model
A main goal of the calibration was to characterize the model fit over a range of CWs defined by the window boundaries ATU 0 and ATU Y .
Functionally, these boundaries were determined by fitting ATU Y for selected δ and calculating ATU 0 with Equation (5). We focused on two CW types covering different phases of embryo development (Velsen, 1980). The O-type CW depicts late phase of organogenesis prior to hatching (ATU range 400 to 520 C Á d) and the E-type CW depicts spawning to alevin emergence (ATU range 0-950 C Á d). Table 3 gives the representative parameters for these window types.  Table 3.
To explore the effect of model fits across the range of CW types we first derived P sets in Figure 3 by limiting the ATU Y search to the range 100 to 947 C Á d (Figure 3). This encompassed the hatch stage, although the search algorithm was not constrained to that stage. The fitted ATU Y versus δ has a hockey stick shape with a break at δ* = 33 days (Figure 3a). The lower boundary, ATU 0 , for the O-type CW clusters at approximately 400 C Á d for δ < δ* and thereafter declines to zero corresponding to the E-type CW (Figure 3b). The ATU M peaks at δ* and then declines resulting in the O-and E-type CWs locating below the hatching range (Figure 3c gold region), ATU H = 515 ± 22 C Á d @ 10-12 C (Alderdice & Velsen, 1978).
Between δ 1 and δ 2 the ATU M is within or above ATU H . Note, the δ* break in ATU 0 and ATU Y graphs (Figure 3a,b) corresponds to the peak ATU M , which is also above the range of ATU H (Figure 3c). The critical temperature for the onset of TDM, T crit , has a hockey stick shape with a breakpoint at δ 1 (Figure 3d). The density-dependent carrying capacity D varies between 70 and 90 redds/km for these model runs

| Reservoir operation simulation
Because all CW types reasonably fit the data, we simulated the effects of alternative reservoir operations paired with their associated factual and counterfactual CEIP assumptions ( Figure 5). Factual simulations for Operation I, which targets cold water over the entire incubation, assume an E-type CEIP and simulations for Operation II, which targets cold water between the hatch of the first and last redds, assume an O-type CEIP. Counterfactual simulations reversed the CEIP assumptions. We compared results by normalizing Tier 2-4 survivals to the T A B L E 3 Egg-to-Fry model results for two critical window (CW) types

| DISCUSSION
In this study, we modeled the effect of temperature during CWs on embryo development and survival. While a CW may be of any duration and centered at any phase of development, the model fit to the two-decade time series of SRWRC salmon incubation resulted in two CW types fitting the data best. The E-type CW covers all incubation, and the O-type covers the pre-hatch period of organogenesis.
Notably, the selection of the CW types has little effect on the predicted TDM. To explore a possible underlying reason for this insensitivity note that b δ is related to δ by Equation (3) such that Equation (2) simplifies to V ¼ exp ÀαΔ ð Þwhere the mean temperature differential over the CW, Δ, is essentially independent of δ (details in Appendix S1B). Thus, the TDM largely is a function of α, which depends on the oxygen flux to the redd, and Δ, which depends on mean temperature over the critical window. Furthermore, we postulate that variability in Δ increases when ATU M is outside the range of the actual TDM processes. Exceeding the range is exemplified by the interval between δ 1 and δ 2 in which the highest opt (Figure 3d) corresponds to the region in which ATU M ≥ ATU H (Figure 3c). In this F I G U R E 3 Plots of parameters P (Table 1) for 111 fitting delimited by pairs of critical window widths (δ) and accumulated temperature units (ATU Y ). In panels a-h critical window (CW) types are shown as clusters O-type ( ) and E-type ( ) with other points not in clusters ( ). Circled points depict Table 3 parameters. Vertical line δ Ã demarks CW width with largest opt (poor model fit) and slope change in ATU versus δ plots (panels a-c). Vertical lines δ 1 and δ 2 demark CW widths where middle of the CW (ATU M ) is equal to or larger than the ATU range for hatching depicted by gold band (panel c  (Rombough, 1994). Also, the model fit improves with declining δ (Figure 3a). In contrast, no direct mechanism is readily identified for mortality processes across an E-type CW extending from fertilization to fry emergence. However, an E-type CW could involve post hatching hypoxic stress associated with microbial degradation of the redd environment (Sear et al., 2016)

| SUMMARY
The paper had three goals; fit the pattern of SRWRC salmon survival over the two-decades of data, explain the pattern in terms of fish embryogenesis and use the framework to guide design and evaluation of Shasta Reservoir operations.
Using the CW framework, we found that both short and long durations CWs better fit the survival data when the window middle was located several days before hatching. Notably, CWs centered in the hatching stage itself yielded lesser fits.
We sought to explain the effect of the CW on model fits using principles of embryogenesis. Literature and the results of this work together suggest the duration of hypoxia and mortality may be a few days to a week prior to hatching. Surprisingly, the models suggest that hypoxia persisting greater than a week would be largely independent of temperature and this insensitivity may account for the low variability in T crit with CW type. Additionally, the CW framework arises in physiological and ecological studies (Mueller et al., 2015;Uchida, Uesaka, Yamamoto, Takeda, & Irie, 2018) and is implicit in studies that consider cross life-stage effects (Del Rio et al., 2021;Gosselin et al., 2021;Stoks & C ordoba-Aguilar, 2012). We suggest that further studies in a CW framework will provide new insights into embryogenesis and population regulation in a warming climate.
The third goal has been to move beyond journal publications and open-source software distribution as the endpoint of model development. We advocate those models actively used in resource management are most useful if implemented as real-time web-based systems available to managers and the public. As described in Appendix S1A, the model presented in this paper is part of the Sacramento River Fish model that is linked to a database (SacPAS: Central Valley Prediction & Assessment of Salmon) containing information from multiple historical and real-time river and fish data sources (e.g., CDWR, 2021) as well as real-time information on reservoir operations and forecasts (CVTEMP, 2022).
In a brief analysis, we explored benefits and risks of applying alternative CW hypotheses with different reservoir operations. The