Identifying minimum freshwater habitat conditions for an endangered fish using life cycle analysis

Identifying the most important factors affecting population growth in animal life cycles is an important tool of species conservation. Delta Smelt (Hypomesus transpacificus), an annual fish endemic to the San Francisco Estuary in California (USA), has been provided legal protection since 1993 but 30 years later exists in a conservation‐reliant state on the brink of extinction. Despite considerable controversies about what factors are most responsible for the species' decline, no population growth rate sensitivity comparisons between the most important factors regulating growth have been done. Nor has anyone attempted to quantitatively identify habitat conditions needed to support positive population growth. We developed a set of stage‐structured population models to link habitat indices regulating recruitment of new generations of fish as they metamorphosed into juveniles and the subsequent survival of those fish over several seasons until they reached adulthood. These models are used to quantify drivers of growth rate variation over 30 years. Several complimentary sensitivity analyses indicated freshwater outflow to the estuary during summer had the largest potential to change population growth. Multiple habitat metrics (e.g., food availability, temperature) influencing recruitment and life stage specific survival rates across different seasons interacted in nonlinear ways to determine habitat conditions and water management targets associated with positive population growth. We discuss the implications for freshwater management, Delta Smelt conservation, and the challenges climate change poses for co‐implementation of these two societal priorities.


| INTRODUCTION
Population growth rates are important components of modern species status assessments aimed at describing the status and trends of populations and providing summaries of current and future factors determining abundances (Smith et al., 2018).Estimates of growth rates and models describing the sources of their variability provide, among other things, quantitative assessment about causes of species abundance changes and support identification of minimum necessary conservation actions intended to reverse declines (Morris & Doak, 2002).For fish residing within watersheds whose freshwater flows are diverted and otherwise altered to support large human populations, growth rate models have the potential to contribute to the oft contentious debate of just how much water they need to persist.The challenges of native fish conservation in arid areas like the North American west are only expected to grow in light of ongoing climate changes (Gido et al., 2023).
One such species whose freshwater habitat has and continues to be altered is Delta Smelt (Hypomesus transpacificus), a small fish (Figure 1) afforded protections under United States federal (USFWS, 1993) and California state (CFGC, 2009) laws.The Delta Smelt is an annual species endemic to the upper portion of the San Francisco Estuary in California, United States.Their spatial distribution and life cycle are somewhat correlated (Moyle et al., 2016), with mature adults found mainly in the Sacramento-San Joaquin Delta (hereafter referred to as the "Delta") and Suisun Marsh and Bay habitats (Polansky et al., 2018).While different rearing strategies have been documented (Hobbs et al., 2019;Lewis et al., 2021), most of the population spends the majority of its life cycle within oligohaline habitats, generally spawning in freshwater areas during the winter and rearing in brackish water areas from spring and early summer through the subsequent winter (Hobbs et al., 2019;Lewis et al., 2021;Moyle et al., 2016).Generally, no individuals are found in the southern Delta as temperatures increase during the summer and fall, and occasionally when flows are high some individuals move west into the eastern areas of the San Pablo Bay (Figure 1).
The San Francisco Estuary's upper portion where Delta Smelt occur broadly consists of the eastern San Pablo Bay, Suisun Marsh and Bay, and the Delta (Figure 1).Substantial changes to both biotic, for example, predator and prey communities (Nobriga & Smith, 2020;Sommer et al., 2007;Winder & Jassby, 2011), and abiotic, for example, hydrologic and physical structure (Gross et al., 2018;Herbold et al., 2022;Hutton et al., 2017;Whipple et al., 2012), features of this portion of the estuary have occurred over recent decades.Most of these changes have been implicated to various extents in the decline of Delta Smelt (Moyle et al., 2016), but particularly germane from a habitat management perspective is the role of freshwater in shaping features of their environment and subsequent survival.While the upper estuary hydrology is complex (detailed in Section 2.1), two simplifications shown in Figure 1 illustrate two primary freshwater hydrodynamic processes contributing to Delta Smelt population dynamics.One pertains to Delta outflow, the amount of freshwater moving from the Delta into Suisun Bay and eventually to the Pacific Ocean, which contributes to the position and function of low salinity habitat important for rearing during the summer and fall (Feyrer et al., 2011;Jassby et al., 1995;MacWilliams et al., 2015;Nobriga et al., 2008).Another relates to the largest point source diversions of freshwater out of the southern Delta by two large pumps, which can create negative net movement of water (i.e., movement in the opposite direction of the Pacific Ocean) that entrains Delta Smelt during the winter and spring into areas they cannot survive (Grimaldo et al., 2009;Kimmerer, 2008).As freshwater passing through the Delta directly supports annual economic activities measured in the billions of dollars and the livelihoods of tens of millions of people, proposed protections involving restrictions to anthropogenic manipulation of the freshwater processes within Delta Smelt habitat has generated considerable legal and political controversy (Scoville, 2020).One contribution of the analysis provided here is to move beyond quantifying how habitat conditions affect recruitment and survival to understanding their consequences for growth rate variability.
Recent population modeling by Kimmerer and Rose (2018) and Smith et al. (2021) has provided some direct quantitative insight into Delta Smelt population growth rates by comparing changes in growth rates between a limited set of scenarios reflecting different habitat conditions.The former study focused on the relative influence of food availability versus entrainment throughout the affected parts of the life cycle.The latter focused mainly on entrainment.However, not done was population growth rate sensitivity analyses to changes in the covariates determining recruitment and survival across all life stages, or how growth rate expected values and variability change over a continuous change in these covariates.Such analyses are needed to compare the relative importance of each of the modeled vital rates (recruitment and survival) along with the covariates that are the best-supported drivers of those vital rates, and to identify habitat condition thresholds delineating negative and positive population growth.These more explicit sensitivity analyses can also help address management questions about relative effects of timing and quantity of freshwater diversions on Delta Smelt's population growth rate.
Here we use a life cycle analysis to quantify how Delta Smelt population growth rates depended on habitat conditions using 30 years of monitoring data on both abundances and environmental conditions.Despite abundance declines at the decadal time scale, not every cohort has had negative growth (Hamilton & Murphy, 2022;Rose et al., 2013).Following the hierarchical population model framework outlined in Newman et al. (2014), we build a set of quantitative population models to evaluate how early life stage recruitment rates and subsequent life stage specific survival rates vary between cohorts with negative and positive population growth, and how changes in each of the underlying covariates controlling these vital rate distributions in turn regulate the population growth rate.The effects of changing multiple covariates simultaneously on the population growth rate are studied to gain more refined insight into the way different factors interact to control the success or failure of each generation, and to identify different combinations of conditions that are predicted to lead to positive population growth.These are used to quantify how much freshwater or other habitat modifications would be needed to stabilize the Delta Smelt population with a given probability.

| Study area and freshwater management
Here we expand on details about Delta Smelt habitat with attention to several key features about the upper San Francisco Estuary, Delta outflow and water diversions from the Delta, to provide management relevant context.The Delta's flow regime has been substantially modified by the storage and diversion of water throughout the watershed resulting in lower flows much of the year.During the summer and early fall however, base flows are elevated beyond what naturally occurred sustained by reservoir releases that repel estuarine salinity to support the export of fresh water from the Delta (Andrews et al., 2017;Gross et al., 2018).Absent water diversions, river inflows into the Delta would become outflows to Suisun Bay.Delta outflow is a net flow (i.e., tidally filtered quantity of water) that is dispersed through a strongly tidal estuarine environment in which the tidal flow is often 1-2 orders of magnitude higher than the net flow (Jassby et al., 1995).Except during very high river flows, this results in the Suisun Bay region experiencing higher salinity than the Delta; Suisun Bay is a frequent location for the estuary's largest low-salinity zone, water with a salt concentration of about 2-5 psu.Delta outflow is regulated in part using a metric called X2, which demarcates the location where near bottom salinity is 2 psu; X2 is reported as a km distance from the Golden Gate Bridge (Figure 1).The statistical relationship between Delta outflow and X2 is inverse because higher river flows push X2 closer to the Golden Gate Bridge, lessening the km distance between the two.
A second important hydrodynamic metric considered here is the sum of Old and Middle river (OMR) flows (Figure 1).The contemporary Delta bears little resemblance to the large floodplain wetland system that existed until shortly after the California Gold Rush began in 1849 (Whipple et al., 2012), and is presently an extensive system of leveed and interconnected canals and shipping channels at the confluence of the Sacramento and San Joaquin Rivers.Most of California's precipitation falls in the north state and consequently inflows to the Delta from the Sacramento River basin average 4-5 times higher than those from the San Joaquin River basin (Kimmerer, 2004).Water is diverted from the waterways of the Delta by numerous public and private entities; however, the export of water by the Bureau of Reclamation's Central Valley Project (CVP) and the California State Water Project (SWP) has the largest hydrodynamic influence on how water moves into and out of the Delta (Hutton et al., 2019).Diversions by these projects often lead to net 'upstream' flow toward the water export diversions (Kimmerer, 2008).These flow reversals are believed to contribute to entrainment of fishes in the diverted water; for Delta Smelt this happens in the winter and spring during which spawning behavior and temperatures contribute to a shift in their population distribution into the Delta region.Thus, combined net OMR flows has been used as a regulatory compliance metric.

| Population model
An overview of the population model is shown in Figure 2a, which tracks four life stages within each annual cohort corresponding to spring, summer, fall, and winter abundances.Denote the abundances by N i,t for stage i and cohort t, t ¼ 1, 2, …, 29 which are linked through time by cohort specific stochastic vital rates.Starting with the adult abundance, the first vital rate is a per-capita recruitment rate r t connecting the number of mature adults in cohort t À 1 to their early life stage offspring in generation t, followed by a sequence of seasonal and stage specific survival rates s i,t , i ¼ 1, 2, 3.Here recruitment is used to refer to the number of early life stage 'post-larval' fish per spawning adult, where post-larval fish are those less than 1 year old in spring through early summer that have swim bladders and ray fins.The population dynamics are N 1,t ¼ N 4,tÀ1 r t and N iþ1,t ¼ N i,t s i,t .
Recruitment is modeled using covariates indexing habitat metrics during late winter and spring time periods; these covariates can be associated with environmental conditions most important for either the spawning adults of cohort t À 1 or their progeny of cohort t, so we have attempted to be clear about which life stage is thought to be influenced by best-supported covariates.The survival rates are modeled with covariates that all apply to the maturing cohort approximately corresponding to habitat metrics during the summer, fall, and wintertime periods.The vital rates are random variables with r t following a log-normal den- , and s i,t following logitnormal density functions f S i,t s i,t ,μ S i ,t x S i ð Þ;σ S i À Á .The mean parameter values depend on covariates in an additive way, where k is the total number of covariates used to model recruitment, but to simplify notation we will describe the general equations assuming a single covariate per vital rate.For these vital rate distributions both the expected values and variances are cohort specific by dependence on the time specific covariate values.For a given set of covariates and realized random effects ϵ Á,t , the realized vital rates are . The state process is thus a multivariate normal random variable after transformation with cohort specific mean (log , and the diagonal entries of Σ are the process variance parameters and the off-diagonal entries are zero.
Two modifications to this population dynamics model were then considered.First, the covariate dependence by r t was replaced with one of two density-dependent functions, the Ricker model In these models, r 0 is an intrinsic growth rate at low abundances, K is an upper threshold at which growth rates become less than one, and a < K is the value below which decreases in N 4,tÀ1 lead to decreases in the per-capita growth rate (depensatory density dependence) conceptualized here as arising because of reduced mate-finding success.Second, in turning our attention to analysis of vital rates, correlation between the s 2,t and s 3,t random effects were evident even though μ S 2 ,t and μ S 3 ,t were not correlated.Thus, we allowed for correlation between s 2,t and s 3,t by specifying where ρ is the correlation between the random effects of s 2,t and s 3,t .(Examination of the correlation between residuals from the top two models indicated that other types of covariance structures were not supported as indi- The realized population growth rate of cohort t is λ t ¼ r t s 1,t s 2,t s 3,t , which is a product of random variables.Because the vital rates depend on covariates which are different for each cohort, the distribution function of λ t is non-stationary with cohort specific mean and variance. Subsequent focus is on the realized, cohort specific, growth rates and distributions conditional on specific environmental covariates rather than an asymptotic distribution for λ t , consistent with acknowledging the importance of non-stationary environmental variation in population growth rate analyses (Koons et al., 2016).A useful expression about λ t is P λ t > λ T ð Þ¼p, where λ T is a target growth rate and p is the associated probability that this target is met.

| Vital rate covariates
The covariates used to model recruitment and survival are summarized in Table 1.These were selected from a more extensive set of ecological and abiotic indicators considered in Polansky et al. (2021) by selecting those T A B L E 1 Model covariate summary grouped by vital rate membership.Data are generally averages over the indicated time period.

Covariate
Predominant implied mechanism and literature support

Recruitment (March-May)
Water temperature a Lower mean temperatures reflect extended spawning season and in turn increasing egg supply (Bennett, 2005).Summarized as average degrees Celsius (deg C).
Adult prey index a Biomass index of adult copepods and mysids.Higher prey densities for adult spawners increase the likelihood of multiple spawns per female, which in turn increases egg supply (Bennett, 2005;Damon et al., 2016).Summarized as biomass per unit volume (BPUV).
Summer survival (June-August) Delta outflow a Higher outflow increases the quality and quantity of low-salinity habitat and relocates it seaward where summer water temperatures are often cooler and water is often more turbid (Bever et al., 2016;Feyrer et al., 2011;Nobriga et al., 2008).Summarized as millions of cubic meters per day (MM m 3 /day).
Secchi depth Higher water transparency (low turbidity) is associated with lower occupancy and higher mortality (Ferrari et al., 2014;Peterson & Barajas, 2018).
Juvenile prey index Biomass index of larval, juvenile, and adult copepods and mysids.Low prey densities during periods of peak water temperature are believed to reduce feeding success and increase mortality (Maunder & Deriso, 2011;Nobriga & Smith, 2020).Summarized in BPUV.
Fall survival (September-November) Delta outflow Similar to summer Delta outflow.
Secchi depth a Similar to summer Secchi depth.

Juvenile striped bass (Morone saxatilis)
Abundance index of competition and predation risk to juvenile Delta Smelt.Summarized as an abundance index.Index calculated the same way as Delta Smelt abundance indices described in Polansky et al. (2019).
Winter survival (December-March) Old and Middle rivers (OMR) flows a Net flow direction in two channels that deliver water to major water diversion systems.OMR flows that are increasingly net negative are associated with higher fractions of the population being entrained into major water diversions (Smith et al., 2021).
South delta Secchi depth a Measure of water clarity in the vicinity of the greatest point sources of water diversion.Delta Smelt are most prevalent in turbid water resulting in greater entrainment risk when turbid water is being exported (Grimaldo et al., 2009(Grimaldo et al., , 2021)).This factor affects entrainment risk synergistically with OMR flow.

Adult striped bass
Abundance index measuring predation risk to adult Delta Smelt constructed using the method described in Polansky et al. (2019).
a Included in reduced model covariate set.
POLANSKY ET AL.
with the highest relative importance and biological reasonableness (i.e., the sign of the slope parameter associated with a covariate was consistent with the expected effect in determining vital rate values).Covariates were also selected to avoid grouping highly collinear pairs within a season (e.g., Delta outflow and the X2 location), although the choice of which to retain is somewhat arbitrary and any particular application may benefit from choosing the covariate of most interest.Data, both predictor variables and abundance indices, were updated from this earlier study to include information from 1991 through 2020, with 2019 being the last cohort for which non-zero abundance indices were available.

| Model types, fitting, and selection
State-space models (Newman et al., 2014) were used for model fitting to address measurement errors and relative biases in the abundance estimates of different life stages.State-space models consist of a state process model describing the population dynamics, and an observation model (described below) linking the estimated abundances to the abundances predicted by the state process model.The state process model set consisted of a total of 10 different models structurally described in Section 2.2, but which differed based on how recruitment was modeled, which covariates were included, and whether or not correlation between fall and winter survival random effects was estimated or fixed at zero (Table 2).With an eye toward parsimony, we compared models with the full set of covariates to those based on a reduced set (see Table 1).The reduced covariate set was selected from the full covariate set based on choosing the top one or two covariates with the largest (in absolute value) slope coefficients for each vital rate process taking all models with the full covariate set into consideration (i.e., different models having in common the full covariate set did not have very different effect size estimates for each of the covariate making a common reduction in this way feasible).Having at least one covariate per vital rate model facilitates estimability of the state-space models constructed here (Polansky et al., 2021(Polansky et al., , 2023)).This was intended to constrain the model dimension explosion resulting from increases in the number of covariates and focus on the most important factors controlling growth rate variation.
The observation model linking abundance estimates b N i,t to state process model predictions N i,t is given by b , where ψ i is a survey specific parameter to account for life stage specific estimate biases, and b . Data from multiple long-term monitoring programs conducted throughout the Delta was used to obtain estimates b N i,t and c CV i,t using the methods described in Polansky et al. (2019).Since 2002 the surveys providing data used for estimating abundance has been different for each life stage, while prior to 2002 the estimates of N 3,t and N 4,t were based on the same survey method carried out in different times of year.The surveys providing data to estimate abundances for stages 2, 3 and stage 4 prior to 2002, were not designed for Delta Smelt, so ψ i was estimated for i ¼ 2, 3 and ψ 4 ¼ ψ 3 for the observation model of mature adult abundances prior to 2002: T A B L E 2 Model features where estimation of ρ refers to inclusion of random effect covariance between the stage 2 and 3 survival rates, number of parameters (k), minimum negative log-likelihood (nll), and sample size corrected Akaike information criteria differences ΔAICc) relative to model m 7 .Models were fit to standardized covariate values using maximum likelihood implemented with Template Model Builder (Kristensen et al., 2016) within the R software environment (R Core Team, 2023).Properties of inference for this type of model without covariance between random effects have been reported on in some detail elsewhere (Polansky et al., 2021(Polansky et al., , 2023)).Model comparisons and selection was based on sample size corrected Akaike information criteria (AICc), one-step ahead prediction residuals (Thygesen et al., 2017), estimated vital rates, and parameter point estimates and standard errors.

| Realized population growth rate sensitivities
Given estimated model parameters from a fitted model, the coefficients β R,1 and β S i ,1 reflect the effects of covariates on vital rates.The relative ability of the past (cohort specific) covariates to change past realized growth rates λ t were compared using the magnitudes of the partial derivatives with respect to each covariate, ∂λ t =∂x Á,t j j .The derivatives were evaluated at the realized covariate values x Á,t and the maximum likelihood estimated parameter values and random effects.Without loss of generality, assume the first covariate of a vital rate is of interest.If the covariate models the recruitment process, this magnitude is β R,1 λ t , which is an exponential curve as a function of x R,1,t .If the covariate models survival process S i , the magnitude is ∂λ t =∂x S i ,1,t j j¼ which is a symmetric function in x S i ,1,t whose center is determined by the sum of the intercept parameter β S i ,0 and realized random effect ϵ S i ,t and which approaches zero asymptotically as the covariate goes to large or small values.One implication is that sensitivity values computed in this way will not necessarily reflect whether a vital rate is limiting population growth in the sense that a small sensitivity may correspond to a survival rate near 0 or 1.The covariate with the largest sensitivity among those modeling a common vital rate will be the one with the largest slope coefficient.

| Population growth rate variation
Analysis of growth rate changes across a continuous range of covariate values was carried out to identify what habitat conditions will likely lead to positive growth.As a product of random vital rates given covariate values, the growth rate is a random variable whose density and cumulative distribution functions, f λ x ð Þ and F λ x ð Þ, respectively, are product distributions that depend on the covariates and vital rate distribution parameters.Given a set of input covariates, vital rates were simulated using the model parameterized with the maximum likelihood parameter estimates.The Delta is an ecosystem with considerable intra-and interannual variation in environmental conditions.As a consequence, process variation overwhelmingly dominates parameter estimate uncertainty in determining predicted vital rate variability (Polansky et al., 2021) so parameter estimate uncertainty was not included.
Two approaches were used to set covariate values from which to predict recruitment and survival.The first simulation method used covariate sets that accounted for correlations and non-linear relationships between those used in model fitting (Figure A.2). Multivariate kernel smoothing (Duong, 2022) was used to estimate their joint probability distribution function f X x ð Þ with which to simulate two million values from the joint distribution of X, after discarding any jointly simulated x that had one or more of its elements outside the range of observed values for the corresponding covariate.This allowed us to sidestep the need to specify phenomenological or mechanistic relationships among the covariates while accounting for any associations.The estimated f X x ð Þ was used to obtain two million simulated values from the covariate joint distribution after discarding any sets that had one or more of its elements outside the range of observed values for the corresponding covariate.A shortcoming of this first approach is that how f λ x ð Þ changes as a function of a one or more covariates of interest will be confounded with the changes to the remaining non-focal covariates according to f X x ð Þ, potentially distorting understanding about the effects of changing the focal covariate of interest alone.A practical consequence of this is that a covariate which has a positive (negative) effect on a vital rate may not show a monotonic positive (negative) relationship with λ t .Additionally, the results will be somewhat sensitive to the choice of using sequential adult life stage abundances to calculate growth rates.For example, if growth rates are based on ratios of sequential fall abundances λ t ¼ N 3,t =N 3,tÀ1 then the covariate set associated with a realized growth rate will be x S 3 ,t , x R,tþ1 , x S 1 ,tþ1 , x S 2 ,tþ1 f g instead of x R,t , x S 1 ,t , x S 2 ,t ,x S 3 ,t f g and the resulting joint density function estimated from these sets f X x ð Þ will be different.The second approach for setting covariate values independently varied each over a grid of values to isolate their effects on growth rate variation.This level of abstraction also removes the dependence of results on the choice of which life stage is used to compute growth rates.We focused on a grid of values emphasizing changes in the hydrodynamic variables summer Delta outflow and winter OMR flows, both of which are particularly important from biological, statistical, climatical, and management perspectives.The grid consisted of 100 equally spaced points ranging between their minimum and maximum values.The remaining four covariates were firstly set at their mean values, and secondly were individually set at either their minimum or maximum values.For each covariate set 10,000 simulated sets of vital rates were generated to capture process variation in these and the growth rates.To condense the basic patterns into a manageable set of results, we first report on growth rate variability over changes in summer Delta outflow and winter OMR with all remaining covariates set to their mean.This is compared with a bookend set of results representing either an overall worst-or overall best-case scenario whereby each of the remaining covariates are chosen so that the corresponding vital rate expected values are at their minimum or maximum with respect to the covariate.We note caution is needed with these predictions because some sets of values on the grid may not be realistic or contained with the range of values used in model fitting (e.g., the maximum spring temperature value paired with the maximum summer Delta outflow value is not near any prior observations).
Growth rate variation across changes in covariates was summarized in two ways.First, the geometric mean of all simulated growth rates associated with a given covariate vector (or covariate vectors lying within a cuboid) was calculated as a measure of the average growth rate.Second, the probability of positive growth P λ t > λ T ¼ 1 ð Þ¼p, calculated as the proportion of simulated growth rates larger than 1, was also examined.We paid special attention to the covariate values where the geometric mean became larger than 1, and where the probability of positive growth is 0.5.While p ¼ 0:5 is not necessarily protective from a conservation perspective (see Section 4), it can be viewed as one interpretation of California's coequal goals in managing the Delta to provide a reliable water supply and for wildlife conservation.We note that the figures provide information about conditions leading to different values of p given λ T ¼ 1.

| RESULTS
Models generally agreed in the qualitative patterns of their residuals and estimated vital rates (Appendix B,Figures B.2 and B.3). Differences were apparent in AICc values, with models using the reduced covariate set and including a correlation parameter ρ between ϵ S 2 ,t and ϵ S 3 ,t random effects having smaller AICc values, and m 7 with density-dependent recruitment without an Allee effect had the smallest AICc overall (Table 2).Parameter estimates changed most noticeably according to whether or not ρ was included in the model, the inclusion of which tended to reduce the winter covariate effect sizes (Table B.1).Early juvenile (N 1 ) carrying capacity K estimates ranged between 16.08 and 17.76 million for density-dependent recruitment models without an Allee effect, and between 3.25 and 3.92 million for densitydependent recruitment models with an Allee effect (Table B.1).The Allee effect threshold a was not well estimated in both models including this parameter (Table B.1, note very large standard errors for the estimates of a).Because of the similarities in models m 7 and m 3 in terms of the different model comparison metrics, and because current adult abundance estimates are much less than any of the estimated carrying capacities, we applied model m 3 to report on relationships between covariates and recruitment, survival, and population growth rates.
Covariates used by model m 3 include spring temperature and the adult prey abundance index for r, summer Delta outflow for s 1 , fall Secchi depth for s 2 , and winter OMR flows and south Delta Secchi depth for s 3 .The distribution of recruitment and survival rates, grouped according to whether the growth rate was negative or positive, indicates that cohorts with positive realized growth have distinctly better recruitment and summer survival rates (Figure 2b).Changes in recruitment and survival over the covariates modeling them are shown in the top panels of Figure 3.The effects of spring temperature and adult food abundance on recruitment were opposite but similar in magnitude, while comparisons of survival rate functions showed that summer Delta outflow has the strongest effect and least uncertainty (Figure 3).We interpret these results to indicate a temperature-driven effect on the duration of the spawning season and that well-fed adults produce more or healthier offspring, that increased habitat quality stemming from increased outflows and higher turbidity positively influences survival over the summer and fall, and that entrainment in water diversions can decrease population level survival during the transition to the adult life stage.
The realized population growth rate mean, median, and geometric mean for model m 3 are 1.41, 0.69, and 0.83, respectively.Of the 29 growth rates, 18 had λ t < 1 and 11 had λ t > 1.The proportion of cohorts with λ t < 1 generally increases through time (Appendix C, Figure C.1).The maximum absolute growth rate sensitivity was associated with spring adult prey three times and summer outflow 26 times (Appendix C, Figure C.1).The finding that summer Delta outflow is the covariate that most frequently had the potential to change past realized growth rates holds across all the other models (Table B.2) and thus appears to be robustly supported.
Growth rate variation in response to changes in covariates when they are sampled from f X x ð Þ are shown in the middle and bottom panel rows of Figure 3.The qualitative direction, but not the exact shape, of vital rate changes as functions of their covariates is generally passed through to the changes in the growth rate distribution (compare panels within columns of Figure 3).The approximate values at which P λ t > 1 ð Þ¼0:5 are 15.8 C for spring temperature, 18 biomass per unit volume (BPUV) of spring adult copepod and mysid (ACM), 25 million (MM) m 3 /day for summer Delta outflow, 68 cm for fall Secchi depth, À190 m 3 /sec for winter OMR flows, and 60 cm for winter south Delta Secchi depth (Figure 3).Realized growth rate changes with increases in spring adult prey and winter OMR flows are clearly not monotonic, partly a result of their larger values being associated with smaller values of summer Delta outflow (Figure A.2).
Growth rate changes in response to changes in the recruitment and survival model covariates using the second method that independently varies each covariate are shown in Figure 4 (also see Appendix C,Figures C.3 and C.4).When the remaining covariates represent the worst-case scenario, summer Delta outflow and winter OMR flows interact to determine the boundary of expected positive growth and its probability of occurring.With remaining covariates at their average values, the importance of winter OMR flows declines (contour lines become approximately vertical), and when they are at a best-case scenario positive growth is always expected with a probability near 1 across any combination of summer Delta outflow and winter OMR flows.Note that the vertical contour lines at about 25 MM m 3 /day for summer Delta outflow, indicating where P λ t > 1 ð Þ¼0:5, matches the findings from simulations based on f X x ð Þ shown in Figure 3.
Detailed examination of covariates and individual realized vital rates indicates vital rates other than summer survival can still be important in regulating the population.The cohorts with the lowest four growth rates also had the lowest four recruitment values indicative of poor per-capita reproductive output or poor survival through the unmonitored larval life stage (Table C.1).Of the 18 cohorts with λ t < 1, the lowest survival rate was during summer, fall, and winter 4, 11, and 3 times, respectively.In each of these cases the product of all vital rates except the minimum survival term was larger than one, so that a sufficient increase in the minimum survival rate could have in theory resulted in λ t > 1. Example cohorts with λ t < 1 where the later life stage survival rates are particularly small include 1994,1996,2005,2013

| DISCUSSION
We used state-space models for fitting population dynamics models to combine imperfect abundance estimates, and this enabled seasonal resolution modeling of covariate impacts to Delta Smelt population growth rates.For the derivative-based realized growth rate sensitivity analysis, we temporarily treated the realized growth rate as a deterministic function of vital rate covariates, akin to the perturbation analyses of periodic models within the matrix modeling literature (Caswell, 2001).Incorporating both covariate uncertainty and process variation makes the population growth rate a random variable with a product distribution whose properties (e.g., means, quantiles) depend on the values of the underlying covariates and random effects.Analysis in this case is made possible using simulations with two perspectives that can be taken with respect to the treatment of the covariate randomness illustrated.These approaches highlighted the ability to extract enough information about the growth rate distribution to identify non-linearities in its response to covariates and what subsets of covariates will provide optimal conservation solutions given competing objectives on resources.
The realized growth rate sensitivity analysis indicated that spring recruitment and summer survival, and the underlying covariates controlling their deterministic values, were the most influential factors affecting Delta Smelt population growth rates.The importance of summertime processes, particularly Delta Smelt body condition as affected by food, agrees with a model sensitivity analysis by Zhang et al. (2022) of an individual based life cycle model developed by Rose et al. (2013).Similarly, greater changes to the growth rate due to changes in summer outflows in comparison to changes in winter south Delta OMR flows were found by Smith et al. (2021).
The growth rate variation analyses broadly supported the sensitivity analysis findings that changes in growth rates and probability of positive growth strongly responded to spring temperature and summer outflow.The strengths of the analytical approach based on derivatives are that it is a classic method with deep understanding about its use in stage structured population models (Caswell, 2001), allows a relatively straightforward comparison of the different rates of change in λ t due to changes in underlying covariates, and provides a "simple" answer for high level communication that is unencumbered by the necessary contextualization involved when discussing random variables.However, the importance of process stochasticity is one of the defining features of the science of population ecology (Bjørnstad & Grenfell, 2001;Newman et al., 2014).Due to the difficulty in applying analysis approaches for stochastic stage structured population models outlined by Tuljapurkar (1990), in part because of recognition of the non-stationary vital rate factors involved here, the simulation approaches developed here allowed practical progress.More generally, the growth rate variation analyses about how λ t and its probability P λ t > λ T ð Þ¼p varied across covariate space provided insight into habitat metrics needed to support positive growth.These also identified values for which further changes result in reduced benefits or detriments to the growth as a result of nonlinearities in the relationships.For example, increases in spring temperatures have relatively little impact until about 15 C, beyond which further increases are associated with marked declines in expected growth rates, while the benefits of summer flow increases begin to taper off once it exceeds about 50 MM m 3 /day.Similar examinations about how λ t and p (for a given choice of λ T ) vary across changes in values of multiple covariates simultaneously allowed us to quantify how covariates interact within and across seasons in determining growth rate variability.The changes are generally complex in that they are nonlinear in any given covariate dimension, are influenced by interactions between covariates determining vital rate distributions of different life stages, and therefore ultimately depend on all covariates across all life stages.
The analyses here can be used to inform minimal target values for habitat metrics, or in support of decisions involving multiple competing objectives (Kennedy et al., 2008).For instance, a goal could be phrased in terms of a minimum desired target growth rate λ T,min or a minimum acceptable probability of positive growth p min .The points or contours (more generally level surfaces) in the covariate space delineating where λ t ¼ λ T,min or p ¼ p min can be used to set these kinds of targets when monotonic relationships exist between the covariate and the measure of growth rate.When there are competing objectives, these values become the best choice in meeting both objectives.Perhaps the most obvious set of competing objectives here is the use of freshwater during the summer to provide outflow in support of Delta Smelt vs. the diversion of that water in support of municipal or other economic uses.Given a choice about a minimum target growth rate and/or a minimum probability of positive growth, the points (contours) satisfying this target are the non-dominated solutions (Pareto frontier) such that any further increases in water used for one objective is detrimental to the other.We note these solutions should be interpreted as approximate because of the large amounts of stochasticity in the system.

| Conclusions and management implications
Over the 20th and into the 21st century the Delta has generally experienced increased spring temperatures and decreased flows (Herbold et al., 2022).But the flow declines are not uniform compared to what would be realized under unimpaired flows (Gross et al., 2018).Flow volumes are generally lower from December-June and higher from July-October.It is possible the lost June outflow that was part of the historical snowmelt hydrograph experienced in California's Central Valley is important to Delta Smelt survival over the summer.If so, increasing outflow in June might increase the frequency of years that the Delta Smelt population grows.However, a recent modeling experiment coupling Delta Smelt bioenergetics with other habitat suitability metrics suggests that foraging conditions deteriorate as the summer progresses and that higher outflows help better align the Delta Smelt's spatial distribution with the best remaining foraging habitat conditions (Smith & Nobriga, 2023; see also Kimmerer & Rose, 2018;Zhang et al., 2022).This implies that outflow after June may be the mechanistic driver of Delta Smelt survival, or that extended periods of very elevated outflow in the summer are necessary for robust population growth given the very altered ecological context in which Delta Smelt now live.
Beginning in December 2021 the United States Fish and Wildlife Service and partners initiated an ambitious experiment on population supplementation by releasing hatchery raised adult Delta Smelt into the Delta.We found no evidence that low spawning numbers are leading to low per-capita recruitment rates; however, we note that detecting Allee effects in marine systems is notoriously difficult (Perälä & Kuparinen, 2017).Similar difficulties likely apply to estuarine fishes as well.Further, the ability to measure the effects of supplementation on the population growth rate is a subtly complicated problem at the intersection of supplementation effort, ensuing environmental conditions, and subsequent sampling program effort and efficiency.Nevertheless, the results here suggest that for any supplementation program to result in a self-sustaining population, some habitat restoration will likely be needed in concert, and some amount of freshwater flow is needed to maintain estuarine habitat.
The findings here suggest summer, not fall or winterspring, is the most important season for freshwater flow augmentation to assist Delta Smelt population growth rate.Our example management target of P λ t > 1 ð Þ¼0:5, is not especially protective; nonetheless, it translates into a daily summer outflow of 25 MM m 3 /day, and historically this has been associated with X2 averaging around 68 km and no greater than 82 km.Of the subset of years from 1991 to 2019 with summer outflow averaging greater than 25 MM m 3 /day during June-August, daily X2 values generally also had distributions between approximately 60 to 80 km (mean = 68 km, min = 52 km, max = 82 km) (Figure 5a).If outflow of 25 MM m 3 /day were used for fish conservation for very long it would quickly meet substantial political resistance and in critically dry years it is not even theoretically achievable (Figure 5b).This emphasizes the importance of multiyear planning to provide stronger protections in wet years when targets can be met to offset declines in dry years.The Delta is in an area of North America that from 2000 to 2021, experienced its driest conditions on record since 800 CE (Williams et al., 2022).Incorporating climate trends and patterns of reoccurring climate scale events (e.g., El Niño) associated with precipitation patterns in California into multi-year projections of population growth rate distributions will be a complex endeavor but should help inform annual targets for λ T and p. Adaptive management within a year can also be applied, starting with records of spring temperature with which to condition subsequent outflow management targets, and similarly fall and winter targets could be set based on early conditions.As climate change alters the portfolio of possible actions in any given year and their likelihood of success, both intra-and interannual planning and forecasting would be important elements of a successfully implemented adaptive strategy.

ACKNOWLEDGMENTS
The viewpoints expressed are those of the authors and do not necessarily reflect the opinions of the U.S.

F
I G U R E 1 Geographic location of the study system along with a photograph of a mature adult Delta Smelt.The major rivers of the Sacramento and San Joaquin river basins that drain into the upper San Francisco Estuary where Delta Smelt occur are shown in the California boundary inset.The bottom detail map includes approximate distances to the Golden Gate Bridge (GGB, vertical lines) and Delta smelt habitat range (outer black polygon).Two important hydrodynamic flow directions, Delta outflow and Old and Middle rivers (OMR) flows, are shown by the arrows.Photo credit: Matthew Dekar (USFWS).

F
I G U R E 2 (a) Life cycle model diagram.(b) Box plots of model m 3 vital rate point estimates across cohorts given negative (λ t < 1) or positive (λ t > 1) realized growth rate.

F
I G U R E 3 Demographic parameters as functions of covariates used to model them based on model m 3 .Top row-recruitment and survival rate medians (solid lines), means (dashed lines), and interquartile ranges (shading), based on 10,000 simulations for a given covariate value.Points are at the realized vital rate and realized covariate values.Middle row-log-population growth rate ( log λ t ð Þ) distributions based on 2 million simulations from the joint distribution f X x ð Þ of realized x values, colored according to density from low (blue) to high (yellow); points are at realized covariate and log λ t ð Þ value pairs.Bottom row-probability that the simulated growth rates in the middle row are larger than 1.

F
I G U R E 4 Model m 3 predicted growth rates (top row) and probability that growth rate is larger than 1 (bottom row) based on 10,000 simulated sets of vital rates given a set of covariate values.Contours are over a range of summer Delta outflow and winter OMR flow values given the remaining model covariate values representing three different scenarios indicated in the panel title.Worst-the remaining covariates are set at the minimum or maximum value so that recruitment or survival is minimized.Average-the remaining covariates are set at their average values.Best-the remaining covariates are set at the minimum or maximum values so that recruitment or survival is maximized.
Department of the Interior, the U.S. Fish and Wildlife Service, or the other member agencies of the Interagency Ecological Program.Editors Toni Lyn Morelli and Brian D. Healy, along with Katherine Osborn and several anonymous reviewers provided helpful comments on earlier versions of this manuscript.The strategy of using kernel density methods for covariate simulation as applied here was suggested by Ken B. Newman.Funding was provided by the California Department of Water Resources and the U.S. Fish and Wildlife Service endangered species program.DATA AVAILABILITY STATEMENTData used in this analysis is available online from the Environmental Data Initiative with package ID=edi.1433.2,Polansky, L., L. Mitchell, and K. Newman.2024.Delta smelt (Hypomesus transpacificus) life cycle model input data.ver 2. Environmental Data Initiative.https://doi.org/10.6073/pasta/356c59a91513a8ff355b3b09b2393c05.ORCID Leo Polansky https://orcid.org/0000-0003-2738-8079Lara Mitchell https://orcid.org/0000-0002-8439-8508Matthew L. Nobriga https://orcid.org/0000-0002-6758-9257FI G U R E 5 (a) Box plots of daily X2 locations during June-August for each year used in the model.Years with an asterisk had average outflow greater than 25 MM m 3 /day.Fill colors are based on whether or not the population growth rate is positive (see legend).Horizontal line is drawn at 82 km.(b) Percent of the total Delta inflow during June-August of 50 years (1970-2019) that would be needed to equal an average outflow of 25 MM m 3 /day over this same time interval, grouped according to the California Sacramento Valley water year index.
Time trends of different covariates as well as correlations between pairs of covariates are shown in Appendix A, Figures A.1 and A.2.