Introduction

The abundance of organisms in a closed population depends on birth and mortality rates. In open systems, abundance will also be influenced by immigration and emigration. Subsidies of nutrients and prey populations from one location of abundance to another of relative scarcity, which is a form of immigration, are common in marine and estuarine ecosystems where exchange rates can be high (Polis et al. 1997). For example, subsidies can be found in marine-derived nutrients transported within spawning salmon migrating to freshwater tributaries (Gende et al. 2002), in larval recruitment (Cowen et al. 2006; Morgan et al. 2017), and in subsidies of zooplankton to inshore regions (Morgan et al. 2017).

In the San Francisco Estuary (SFE), productive shoals provide spatial subsidies of phytoplankton to light-limited channels (Lucas et al. 1999a) and areas of low benthic grazing with higher phytoplankton densities can subsidize areas of high benthic grazing with lower densities (Lucas et al. 1999b, 2009; Kimmerer and Thompson 2014). Notably, the introduced brackish-water clam Potamocorbula amurensis has grazed down phytoplankton biomass in the low-salinity zone (LSZ, salinity 0.5–5, Practical Salinity Scale) during most springs and every summer since 1986 (Alpine and Cloern 1992; Cloern and Jassby 2008; Lucas et al. 2016). This grazing contributes to a persistently low phytoplankton biomass that is partially offset through spatial subsidies driven by hydrodynamic transport (Kimmerer and Thompson 2014).

Populations of calanoid copepods in the LSZ in summer–autumn have also been reduced through a combination of the lower phytoplankton densities and high consumption of copepod nauplii by clams (Kimmerer et al. 1994; Kimmerer and Lougee 2015). The most abundant calanoid copepod in the LSZ in summer is Pseudodiaptomus forbesi, which maintains a population in the LSZ only through a subsidy from its freshwater population center in the Sacramento-San Joaquin Delta (Delta; Fig. 1; Kimmerer et al. 2018, 2019). P. forbesi is a small (~1 mm as adult) copepod introduced from subtropical regions of Asia (Orsi and Walter 1991), and is abundant from ~May to November in freshwater. The high-salinity limit of its distribution is probably not controlled by salinity tolerance, but rather through predation by clams and predatory copepods (Kayfetz and Kimmerer 2017).

Fig. 1
figure 1

Map of the northern San Francisco Estuary with boundaries of the low-salinity box for X2 scenarios of 67, 74, 81, and 85 km. Black triangles represent sampling stations, red dots indicate the water diversion facilities of the State Water Project and federal Central Valley Project, and the grey border line is the boundary of the legal Delta. Blue shading indicates depth in meters. The red square of the inset map indicates the position of the upper estuary in California

Full size image

The LSZ provides habitat for a plethora of fish species (Kimmerer et al. 2013), but many of them have declined in abundance or shifted their distributions toward habitats with presumably greater food availability (Sommer et al. 2007). Most of these fishes consume zooplankton during larval to juvenile stages, and although other copepods can be sporadically important, P. forbesi is the most commonly consumed prey, particularly for critically endangered delta smelt, Hypomesus transpacificus (Slater and Baxter 2014; Slater et al. 2019). The majority of juvenile delta smelt rear in the LSZ (Brown et al. 2014; Bush 2017), and the decline in delta smelt has been attributed in part to poor feeding conditions there, particularly during summer and autumn (Mac Nally et al. 2010; Slater and Baxter 2014; Hammock et al. 2015). The decline of this fish has continued to the point where it is seldom caught in routine monitoring surveys that historically captured delta smelt commonly and in large numbers.

In September–October of recent wet years, net Delta outflow (calculated seaward flow of freshwater at the western margin of the Delta, hereafter “outflow,” CNRA 2021), has been augmented by increasing flows from reservoirs or by reducing diversions (“exports”) from the tidal freshwater reach of the Delta. Augmentation of outflow is one of several measures instituted by state and federal agencies to improve habitat in response to a decline in delta smelt abundance. Regulatory permits for state and federal water diversions now include a requirement for flow augmentations intended to increase food availability for delta smelt in the LSZ during summer and autumn (USFWS 2019a; CDFW 2020).

The objective of this study was to estimate the magnitude of the subsidy in P. forbesi that resulted from the outflow augmentation. We estimated subsidies of P. forbesi into the LSZ during 2017, a year of relatively high outflow. Two previous studies provide background for this study. The first showed that abundance of P. forbesi in the LSZ was higher during summers with high freshwater outflow (Fig. 4 in Kimmerer et al. 2018). The second analyzed abundance data with a box-modeling approach to estimate spatial variability in mortality of P. forbesi (Kimmerer et al. 2019). That study found that mortality of the nauplius larvae was exceptionally high in the LSZ, suggesting that abundance in the LSZ can be maintained only through spatial subsidies from other locations. This study builds on Kimmerer et al. (2019) in two ways. First, the focus of Kimmerer et al. (2019) was mortality and its potential causes. Here, we are more concerned with estimating the magnitudes and effects of the subsidies, and thereby evaluating the effectiveness of flow augmentation for this purpose. Second, the previous study used mainly fixed sampling stations to populate the spatial boxes with abundance data, whereas we sampled for copepods using a random sampling station assignment to better characterize spatial patterns of abundance and thereby the magnitude of the subsidy.

Methods

Our study uses salinity and geographic position to characterize the distribution of the copepods. Estuarine plankton have long been known to sort along the salinity gradient so that each species is most abundant within some range of salinity (Jeffries 1962; Collins and Williams 1981; Miller 1983). This can reflect salinity tolerance or the effects of predators on population distributions (Rippingale and Hodgkin 1975), including that of P. forbesi (Kayfetz and Kimmerer 2017). Regardless of the causes of these patterns, salinity is a more reliable guide than location to the distributions of planktonic species in estuaries with highly variable river flow.

Our box model approach followed Kimmerer et al. (2019) by dividing the upper estuary into a series of boxes with boundaries determined by salinity and examining population dynamics and movement between boxes under steady-state conditions within a box. Since the abundance of P. forbesi is closely tied to salinity, having spatial boxes with similar salinity patterns among flow scenarios made the model calculations more straightforward and interpretable than they would have been with spatially fixed boxes. However, we used a generalized random tesselation stratified spatially balanced survey design (Stevens and Olsen 2004) to collect zooplankton samples and assign each sample to a box based on the salinity at the sample site; thus, the assigned boxes could differ among flow scenarios. Three-dimensional hydrodynamic and particle-tracking models were then used to estimate the tidal and net movement of water and copepods between boxes, and a simple population-dynamics model was used to assess mortality and the magnitude of the subsidy.

We modeled movements and subsidies under four values of outflow as indexed by X2, the distance up the main channel of the estuary to where the tidally and daily averaged near-bottom salinity was 2 (Jassby et al. 1995; MacWilliams et al. 2015). The four flow scenarios were run at X2 values that covered the typical range for autumn: 67, 74, 81, and 85 km, corresponding to steady outflows of 536, 348, 235, and 190 m3 s−1, respectively (based on Eq. (8) with parameters in Table 5 of MacWilliams et al. 2015).

Assumptions of this method include the following: (1) salinity is a suitable predictor of the distributions of copepods; (2) a steady-state approach is adequate; (3) mortality for a given sample is constant over a series of copepod life stages, e.g., all nauplii or all copepodites; and (4) making the spatial (salinity) distributions discrete does not overly distort the dynamics. The first assumption is discussed above. The steady-state assumption is appropriate to eliminate effects related to transient adjustments to flow and minor effects such as spring-neap tidal variability that would be associated with a specific point in time (MacWilliams et al. 2015) and, instead, provide a representative estimate of daily exchange between boxes for each X2 value. The age (life stage) distribution of copepods within each sample is also assumed to be stable (Kimmerer 2015). An assumption of constant mortality across several life stages is necessary to be able to model mortality based on differences in abundance across life stages (Kimmerer 2015). Box modeling for examining spatial variability in estuaries has a long history of characterizing distributions and transformations of properties (Officer 1980; Officer and Nichols 1980; Testa and Kemp 2008). The alternative to a discrete (box) model is a full individual-based model, which is beyond the scope of our study.

Zooplankton Sampling

We sampled zooplankton in the northern estuary during four surveys in the autumn of 2017: Oct. 11 to 16, Oct. 23 to 26, Nov. 6 to 10, and Nov. 20 to 22 (Fig. 1). Mean X2 during these events was narrowly constrained to 75, 76, 78, and 77, respectively, corresponding to a narrow range of steady outflows from 277 to 328 m3 s−1, calculated as above. Up to 10 sites were sampled per day following a generalized random tessellation stratified sampling (GRTS) design (Stevens and Olsen 2004), in which the northern estuary was divided into five polygons. Half of the sites were in shallow water (< 3 m) and the other half were in deeper water (> 3 m) based a digital elevation model created for the Delta using bathymetry and topography data relative to the North American Vertical Datum of 1988 (Fregoso et al. 2017). We chose 3 m as a depth threshold because it is close to the median depth of the estuary. Salinity was measured at the time of each sample to assign zooplankton data to spatial boxes (below). At each site, a vertical tow was taken with a 50 cm diameter, 50 µm mesh net pulled from the bottom, or at most from a depth of 10 m, to the surface. The volume of water filtered was calculated as the product of the area of the opening, the distance towed vertically through the water column, and a net efficiency of 70% (Kimmerer and Slaughter 2016). Upon retrieval, the net was washed down and the sample concentrated into labeled jars with formaldehyde at a final concentration of ~4% by volume.

In the laboratory, subsamples of the zooplankton samples were taken using an adjustable pipette, and in each subsample the life stages of P. forbesi were identified and counted as in Kimmerer et al. (2018). Adults were sexed and counted, copepodites were staged and counted, and nauplii and eggs in both loose and attached egg sacs were counted. To provide high enough counts of each stage for mortality estimates, we took additional subsamples to obtain at least 30 of each life stage. Abundance (m−3) was determined by dividing the count of each life stage by the product of the subsample fraction for that life stage and the volume filtered in the original sample.

Box Modeling

Our box model approach closely followed methods in Kimmerer et al. (2019; Fig. 2). The three seaward boxes represented the following: (i) very low-salinity zone (VLSZ, salinity 0.2–0.5); (ii) low-salinity zone (LSZ, 0.5–5.0); and (iii) high salinity zone (SAL, > 5.0). The box boundaries varied in location with X2, which indexes the approximate center of the LSZ box. The two most landward boxes (salinity < 0.2) were the Sacramento River box (SAC) and the San Joaquin River box (SJ). Both adjoined the VLSZ box at its landward margin at salinity 0.2.

Fig. 2
figure 2

Flow diagram of calculations. Shaded shapes represent calculations made for this paper, whereas copepod development time is unshaded because it involves parameter estimates from Kimmerer et al. (2018). Hexagons are three-dimensional simulation models and rectangles are the results of calculations based on data and 3-D model output

Full size image

Hydrodynamics were simulated using the three-dimensional RMA UnTRIM San Francisco Estuary Model (Gross et al. 2010; Andrews et al. 2017; Kimmerer et al. 2019). UnTRIM solves the discretized Reynolds-averaged shallow-water equations on an unstructured grid (Casulli and Walters 2000). It allows for the wetting and drying of computation cells, and a sub-grid scale representation of bathymetry (Casulli and Stelling 2011). Vertical turbulent mixing in the model is parameterized using a k-ε closure, which solves one equation for turbulent kinetic energy (k) and another for turbulent dissipation (ε) using published parameter values (Warner et al. 2005). The model domain extends from the coastal ocean through the San Francisco Estuary including the Delta. The model was applied to simulate hydrodynamics from October 18, 2016 to December 1, 2017, using idealized tides with the M2 tidal period modified to 12 h, so that tides repeated daily and therefore did not vary on a spring-neap cycle. The wind also repeated a daily pattern, which was representative of average summer and autumn wind conditions. This approach allows transport to be estimated for a single day of representative tidal and wind conditions (Kimmerer et al. 2019).

The Flexible Integration and Staggered-grid Hydrodynamics Particle Tracking Model (FISH-PTM; Ketefian et al. 2016) used calculated velocities and eddy diffusivities from the UnTRIM model at a half-hour interval to simulate individual particle trajectories in the study area through time. These particles represented P. forbesi under two alternative swimming behaviors: tidal migration for copepodites and adults, and passive behavior for nauplii. Tidal migration was simulated using an upward swimming speed of 0.25 mm s−1 during flood tide and a downward swimming speed of 0.75 mm s−1 during ebb tide, which resulted in the closest match between simulated particles and vertical distributions of copepods observed in the estuary (Kimmerer et al. 2014). The modeled tidal migration behavior is effective at retaining copepods in regions of salinity stratification, but less so in freshwater regions. Nauplii (< 0.2 mm long) were assumed to be too small to swim continuously at speeds required to overcome vertical turbulence and were therefore treated as passive (neutrally buoyant) particles.

The FISH-PTM model was initialized by seeding each box with a spatially uniform density of one particle per 5000 m3, resulting in approximately 400,000 particles for each of the eight combinations of simulated behavior and X2. The FISH-PTM then simulated forward over 24 h (two semi-diurnal tidal cycles). Particles were transported between boxes by tidal and net currents. After 24 h, the number of particles present in a destination box from each source box was divided by the original number of particles seeded in the source box. This fraction represented the proportional exchange from a source box to a destination box. The calculation was repeated for all source and destination boxes to populate a “Gross Exchange Matrix” (Kremer et al. 2010). Losses of particles at water diversion locations (sinks; Fig. 1) were also calculated and these particles were removed from the simulation.

The two swimming behaviors were applied for each of the four X2 scenarios, resulting in eight Gross Exchange Matrices. Omitting X2 and life-stage subscripts for simplicity, the proportional loss was calculated from

$${E}_{bjk}=\frac{{P}_{bk, t=1}}{{P}_{bj,t=0}}$$
(1)

where the gross exchange matrix element Ebjk is the proportion of particles with swimming behavior b that originated from source box j and moved to destination box k within one day. Pbj,t=0 is the number of particles seeded in box j during initiation at time zero, and Pbk,t=1 is the number of particles in box (or sink) k one day after the release time in box j.

Subsidy Calculations

The proportional subsidy and net gain of copepods to each box is described below as a proportion of the steady-state density in the box. The gain to each box, i.e., the ratio of the subsidy (copepods day−1) to the extant population size in the box, depends on exchange rates and differences in copepod densities and volumes among boxes. Following the notation in Eq. (1), the gains (or subsidies) were calculated from

$${g}_{bkj}=\frac{{n}_{k}{E}_{bkj}{V}_{k}}{{n}_{j}{V}_{j}}$$
(2)

where gbkj is the proportional gain in numbers in box j due to transport from box k, nj and nk are the mean densities of P. forbesi (m−3) across all sites and surveys in a box, and Vj and Vk are the volumes (m3). The exchange proportions and box volumes in Eq. (2) vary with X2, but the mean population densities within each box do not because they were calculated from the zooplankton samples. Density in the freshwater boxes was historically invariant with freshwater flow, but density in the LSZ was elevated at flows above about 300 m3 s−1, corresponding to our two higher-flow scenarios (see Fig. 2b of Kimmerer et al. 2018). This means that under the 67 km X2 scenario and to a lesser extent the 74 km scenario, the magnitude of the subsidy was somewhat overestimated by our method. The proportional net gain into box j is the sum of the subsidies (Eq. (2)) from all other boxes combined, less the sum of the proportional losses (Eq. (1)) from box j to the other boxes:

$${G}_{bj}=\sum\nolimits_{j\ne k}{g}_{bkj}- \sum\nolimits_{j\ne k}{E}_{bjk}$$
(3)

Mortality Calculations

Mortality rate in a closed population is typically calculated from the rate of change in density from one stage to the next, accounting for the durations of each stage. Copepods develop through six nauplius stages, five copepodite stages, and the terminal adult stage. However, robust estimates of mortality must take into account the uncertainty in both the stage durations and in the density estimates, which are based on counts of individuals and therefore have a variance that scales with the number counted (Kimmerer 2015). In practice this means estimates of mortality must be made across a sequence of life stages, e.g., all copepodites, requiring an assumption that mortality is indistinguishable among all stages within the group, and the life-stage distribution is stable (Kimmerer 2015). Violations of these assumptions manifest as uncertainty in the mortality rate, but mortality in a dynamically changing population can be estimated only through population-dynamic modeling, which is difficult and rarely done.

An additional challenge is that real copepod populations in an estuary are not closed, so mortality calculated as in Kimmerer (2015) is biased by emigration or immigration. Thus, we use the term “apparent mortality” to mean the mortality calculated under the assumption of a closed population, which is subsequently corrected by estimates of the actual movements from the box model (Kimmerer et al. 2019).

We calculated apparent mortality within each box from counts of individuals by life stage, fractions of the sample counted for each life stage, and estimates of the stage durations under the above assumptions (Kimmerer 2015; Kimmerer et al. 2019). We used a Bayesian hierarchical approach to estimate apparent mortality for copepodites (based on abundance by stage and estimated stage durations) and adults (based on the results for copepodites and the abundance of copepodites and adults) following Kimmerer and Lougee (2015). Apparent mortality could not be calculated for nauplii, because counts of nauplii were lower than expected from the counts of the other life stages, likely as a result of predation by Acartiella sinensis (Slaughter et al. 2016), a predatory copepod that is most abundant at salinity of ~1–10.

The actual mortality rate, which we term the in-situ mortality rate, Mbi, was calculated as the sum of apparent mortality mi and the proportional net gain Gi based on:

$${M}_{bi}={m}_{i}+{G}_{i}$$
(4)

where “i” is an index for each unique site-survey sample. The net gain estimates essentially correct the apparent mortality estimate for gains and losses resulting from transport. Gi values were computed by linear interpolation, using the box-specific values from Eq. (3), with X2 values that spanned the X2 value at station ‘i’ at the time of sampling.

Uncertainty Estimates

Uncertainty in estimates of proportional subsidies, proportional net gain, and in-situ mortality were computed by bootstrapping (Efron and Tibshirani 1986). The subsidy and net gain calculations are based on mean P. forbesi densities for each box across sites and surveys (Eq. (2)). However, these means are uncertain owing to variation in estimates of mean density in each box on each survey. We computed the mean and the standard error of the mean for the log density for each box and survey and used these values in the rnorm() function in R (R Core Team 2019) to return random samples of the log of mean densities which were then transformed. These values were used in Eq. (2) to compute uncertainty in the proportional subsidy, and therefore net gain. Uncertainty in in situ mortality estimates depends on both the uncertainty in apparent mortality determined from the Bayesian model, and uncertainty in net gain. We combined these uncertainties by using random samples from the posterior distribution of mi, and the distribution of Gi values from boostrapping (Eq. (4)).

Results

The particle-tracking model predicted very little transport from seaward boxes (LSZ, VLSZ) into the San Joaquin and Sacramento boxes owing to strong net seaward flow (Fig. 3). Exchange from the VLSZ box to the LSZ box was substantial and increased with outflow (i.e., lower X2). Less than 10% of particles in the low-salinity box were transported daily to the SAL box (Fig. 3).

Fig. 3
figure 3

Predicted proportion of particles leaving each source box (panels) and going to destination boxes (colors) under four X2 conditions for particles with tidal migration behavior. The sum of exchanges from each source box to all destination boxes equals one, except in cases where there are losses due to water diversions and other particle sinks. Boxes are lower San Joaquin River (SJ), lower Sacramento River (SAC), very low-salinity zone (VLSZ), low-salinity zone (LSZ), and high salinity (SAL)

Full size image

Population densities of P. forbesi based on the zooplankton samples were highest in the SAC, SJR, and VLSZ boxes, and were much lower in the LSZ and SAL boxes (Fig. 4). Within boxes, densities of copepodites were higher than those of adults, as expected. Densities of nauplii were sometimes lower than those of copepodites, which is impossible in a closed population because the development times for each of those life stages are about the same; therefore this likely reflects the strong subsidy of this life stage from the lower-salinity waters. There was considerably more variation in P. forbesi densities across stations within boxes than there was across surveys.

Fig. 4
figure 4

Box plots showing the variability in P. forbesi densities across boxes by life stage averaged across sampling surveys. The center line in each box plot shows the median and the box shows the interquartile range or IQR. Error bars show Q1-1.5*IQR and Q3+1.5IQR

Full size image

Subsidy calculations depend on box volume (Eq. (2)), which varied among boxes and among X2 levels (Table 1). The volumes of SJR and SAC boxes at the highest X2 level were respectively 30 and 40% lower than at the lowest X2 level due to the landward contraction of the boxes as X2 increased. The volumes of VLSZ box and LSZ box varied less but were not monotonic across X2 levels.

Table 1 The calculated volume of each box (millions of m3) at the four modeled X2 scenarios with corresponding net Delta outflow (NDO). Boxes are lower San Joaquin River (SJ), lower Sacramento River (SAC), very low-salinity zone (VLSZ), low-salinity zone (LSZ), and high salinity (SAL)
Full size table

Proportional subsidies of P. forbesi varied by life stage, box, and X2. Proportional subsidies to the LSZ for nauplii were very high but uncertain (Fig. 5). Proportional subsidies for copepodites and adults were more precise and often large. Subsidies declined with increasing X2 (i.e., less outflow) in the LSZ because of reduced exchange from the VLSZ to the LSZ (Fig. 3).

Fig. 5
figure 5

Proportional subsidies and net gains of P. forbesi (Eq. (2)) to the low-salinity zone (LSZ) by X2 scenario. Results are shown by life stage (columns) averaged across sampling surveys. The box plots are presented as in Fig. 4, with outliers omitted to show X2 patterns more clearly

Full size image

Proportional net gains to the low-salinity box were greatest for nauplii and smallest for adults with similar patterns of uncertainty to those for proportional subsidies (Fig. 5). In the LSZ, the median net gain declined with increasing X2 for copepodites and adults.

In-situ mortality rates for copepodites in the VLSZ and LSZ boxes were often considerably greater than in the SAC and SJ boxes (Fig. 6). This pattern was less evident for adults. Estimates of apparent mortality for P. forbesi nauplii were unreliable because counts were lower than required to support the rate of production estimated for the first copepodite stage. Also, estimates of net gain were very uncertain (Fig. 5), so in situ mortality could not be reliably estimated for this stage.

Fig. 6
figure 6

Station-specific estimates of in-situ mortality for P. forbesi as a function of salinity by life stage and sex. The color of points denotes the box the sample was taken in. Points and error bars show the medians and 80% credible intervals, respectively

Full size image

Discussion

Flow Augmentation Supports Food Subsidy to the Low Salinity Zone

The U.S. Fish and Wildlife Service mandated augmentation of flow in autumn to provide “direct and indirect benefits to delta smelt” under an implicit and untested assumption that it would augment zooplankton prey availability for delta smelt in their low-salinity habitat (USFWS 2008). However, it is not feasible to directly measure zooplankton subsidy in the San Francisco Estuary. An earlier study conducted nine 30 h intensive field sampling events over three years using two boats in two of those years to try to infer movements of copepods from their distributions. They observed tidally phased vertical positioning of copepods and larval fish but could not resolve longitudinal movement (Bennett et al. 2002; Kimmerer et al. 2002). That movement was finally resolved using particle-tracking modeling in which the particles mimicked the observed vertical movement of the copepods (Kimmerer et al. 2014).

Using a box model, this study was designed to test that assumption for P. forbesi, which comprise ~40% of the diet of delta smelt in late spring through summer (Slater and Baxter 2014). The autumn abundance of all life stages of P. forbesi in the LSZ is much lower in brackish than in freshwater habitats (Kayfetz and Kimmerer 2017, Fig. 4 in this study). Our results show that subsidy of copepods to the LSZ increases with outflow augmentation (i.e., lower X2). Higher subsidies to the low-salinity box with greater outflow are partially offset by higher in situ mortality. Net gain of copepods within the low-salinity box increased with outflow; however, increased subsidies to an area may or may not lead to a measurable increase in abundance, which depends on the past history of all population-dynamics processes and can be obscured by high sampling variability.

Proportional Gains and Losses

Subsidies at high flows provide a large contribution in proportion to the low abundance in the receiving environment, but proportional losses from the source environment are low because recruitment of copepods is large, so losses through subsidy are a minor component of population dynamics in freshwater source regions. Using adults as an example, at X2 = 74 km the SJR box had a volume of 6 × 108 m3 and a proportional loss rate to the VLSZ box of ~0.04 day−1, which had a volume of 0.9 × 108 m3 and a proportional subsidy estimated at 0.49 day−1. Two factors come into play: [1] the larger volume and [2] the higher abundance of copepods in the SJR box than in the VLSZ box. Just as augmented flows increase subsidy to the LSZ, they also increase subsidy to the VLSZ because copepod abundance is highest in freshwater regions. Moreover, since abundance increases in the LSZ as flow increases, the most parsimonious conclusions for the VLSZ is that abundance increases slightly there as well. This would be difficult to detect, though, because few monitoring stations were within that small area. This was explored more fully in Kimmerer et al. (2019), in which proportional losses from the SJR box were similarly small.

Using the same zooplankton box model for different years and data from sampling programs with three different designs, Kimmerer et al. (2019) also found that mortality rates of P. forbesi copepodites in the low-salinity box were higher than in other boxes. These findings are consistent with the concept that predation by invasive clams and copepods (A. sinensis) depletes zooplankton densities in low-salinity habitat (Kimmerer et al. 1994, 2018; Slaughter et al. 2016; Kayfetz and Kimmerer 2017), increasing the importance of subsidies from landward sources in repopulating the LSZ with copepods.

Effect of X2

The calculated subsidy for the box model is not influenced by assumptions related to mortality. Rather, it is strongly driven by the large observed density difference between freshwater boxes and the other boxes and subsidy has been observed to be salinity dependent, not geographical (Kimmerer et al. 2018). Our expectation was that the subsidy from zooplankton-rich freshwater sources to the LSZ would be inversely related to X2, and therefore proportional net gain would also be inversely related to X2, i.e., directly related to outflow. Our modeling results support this expectation and indicate that subsidies of P. forbesi declined with increasing X2. The lowest X2 modeled (67 km) has occurred infrequently in the historical record for autumn (CNRA 2021), but was chosen to expand the range of flows to improve our ability to detect effects. Moreover, the subsidies were expected to bolster the population in the LSZ, which is consistent with higher abundance in the LSZ when outflow is high, as was observed in 2011 (Kimmerer et al. 2018).

Utility of Box Models

Box models are a valuable tool for analyzing spatially variable data in a context where the details of the spatial distributions or the underlying dynamics are uncertain (Officer 1980). They sacrifice temporal or spatial resolution for tractability; the alternative would have required a spatially explicit, detailed population dynamic model. Neither rate processes nor movement of zooplankton and other small, rapidly growing organisms can be inferred from abundance data alone. Our mortality estimates used a box modeling approach to distinguish in situ mortality from apparent mortality by accounting for hydrodynamic transport (Kimmerer et al. 2019). We quantified uncertainty in both subsidies and in situ mortality which resulted from large variability in counts among samples. Although increasing the resolution of predictions by decreasing box size may improve representation of transport processes, the resulting tradeoff would increase uncertainty in subsidy and mortality calculations by decreasing the number of observations per box. The approaches used to calculate subsidy and mortality in our study can be applied to other organisms and other ecosystems given sufficient abundance data in each spatial box and an understanding of swimming behavior. However, a calibrated hydrodynamic model of the aquatic system of interest is required to drive a particle-tracking model with an appropriate behavioral representation of the organism of interest.

Some of the inputs to the box model are well characterized while others are uncertain. First, the volumes of modeled boxes are known accurately, and the exchange of particles between boxes can be estimated reliably, using the well-calibrated and validated hydrodynamic and particle-tracking models used in this analysis (Ketefian et al. 2016; Gross et al. 2019). The particles were either passive to represent the small nauplii, or moved up in the water column on the flood and down on the ebb to represent copepodites and adults. This behavior results in a good match with the vertical distributions of the copepods (Kimmerer et al. 2014) and causes retention of particles that reach the deeper, more saline parts of the estuary. Second, the estimates of proportional net gain are largely driven by differences in zooplankton densities among boxes.

Spatial Transport Processes

This study builds on that of Kimmerer et al. (2019) to reinforce the finding that spatial transport processes are an essential element of the population dynamics of P. forbesi, and probably other planktonic taxa. These processes are now seen to be important at a range of spatial scales, from the entire habitat of the species to the scale of individual wetlands and their connections to open water (Kimmerer et al. 2018; Yelton et al. 2022). Transport processes of plankton are heavily dependent on the interaction of tidal flows with spatio-temporal variation in abundance and migratory behavior, presenting a challenge for analysis. Given the availabity of robust, well tested transport modeling described above, and a growing knowledge and understanding of the population ecology and behavior of P. forbesi, a logical next step would be to derive a more mechanistic model of copepod production and subsidy. This would combine transport modeling with a better representation of the vertical movements of copepods, reinforced by data on abundance and vital rates by life stage, to develop spatial mass-balance models. Other species of copepod are known to be important to the diets of small fishes and can be dominant in certain years and seasons including autumn (Slater et al. 2019; Slater and Baxter 2014). Models and analyses incorporating a more complete picture of food items important to delta smelt may improve predictions of how flow augmentations may improve foraging habitat for them.

Management Implications

The implication of these findings for water resource managers is that increasing freshwater outflow to reduce X2 does appear to increase subsidies and net gain of P. forbesi to the LSZ by twofold between an X2 of 85 and 74 and threefold between an X2 of 85 and 67. Average X2 values in wet years since 1995 are 77.2 in September and 79.9 in October; thus, realistic increases in flow and associated subsidies would be more modest. A flow augmentation to achieve an X2 of 67 for 1 month would be impractical, in that it would require about 16% of the capacity of Lake Shasta (the largest reservoir in California), and would cost ~ $100 M based on a steady-flow analysis of MacWilliams et al. (2015) and current water prices, and is likely an underestimate if the next year is dry. However, even more modest augmentations should increase abundance of P. forbesi in the LSZ, but the typically high variability of sampling for zooplankton will limit detectability of an increase to times of very high flow (Kimmerer et al. 2018). In addition, planktivorous fishes and other planktivores may respond to a slight increase in prey abundance by increasing their feeding rate or by shifting their distribution into the LSZ, masking the increase in zooplankton density. Since increasing feeding success is one of the main goals of flow augmentation, the efficacy of the augmentation should be evaluated using a modeling approach like the one presented here or a more complex individual-based model for target fish species such as delta smelt.

Multiple flow-related management actions to benefit delta smelt are specified in the most recent operating permit allowing the federal water diversions to take a limited number of the species (USFWS 2019b), with most of these actions targeting enhancement of food. Based on historical data, delta smelt distribute across the salinity gradient, but occur at higher density in the LSZ than in fresh water. Yet, the concentration of the zooplankton food of delta smelt in summer is much higher in fresh water than in the LSZ (Kimmerer et al. 2018). While studies of delta smelt condition have revealed higher indices of feeding success in fresh water and in Suisun Marsh, and lower indices in Suisun Bay, life in fresh water is more energetically costly (Hammock et al. 2015, 2022). Although most of the delta smelt population migrates to the LSZ to rear as juveniles before moving back to freshwater areas upstream prior to spawning in the late winter and early spring, the life history strategies of delta smelt appear to be more varied than once thought (Hobbs et al. 2019). The approach we applied in this paper may provide a means to evaluate the response of food to flow-related actions in other areas of interest.