Modeling the Influence of Outflow and Community Structure on an Endangered Fish Population in the Upper San Francisco Estuary

The aim of this community modeling study was to evaluate potential mechanisms by which freshwater outflow in the upper San Francisco Estuary, CA, controls the fall habitat and abundance of subadult delta smelt Hypomesus transpacificus and its community. Through analyses of the community matrix, community stability and the direction of change of community variables were qualitatively and quantitatively modeled under four outflow–input scenarios. Three subsystems were modeled in the low salinity zone (1–6 psu), each overlapping the location corresponding to the distance from the mouth of the estuary to upstream positions where the near-bottom 2 psu isohaline (X2) is at 74, 81, and 85 km (corresponding to high-, mid-, and low-outflows). Results suggested communities were qualitatively stable at each X2 position, but simulations showed the percent of stable models decreased from lowto high-X2 positions. Under all outflow–input scenarios, the predicted qualitative population responses of delta smelt were: (1) consistently positive for the low X2 position, and (2) uncertain under both midand high-X2 positions. Qualitative predictions were generally consistent with quantitative simulations and with the relations between relative abundance of delta smelt and X2. Thus, high outflow seems beneficial to subadult delta smelt when X2 reaches 74 km during fall.


Introduction
The complexity of estuarine ecosystems is influenced by the broad range and extent of anthropogenic effects in these transitional waters [1,2], including morphological and hydrological alteration [3,4], entrainment of plankton and nekton to water pumps [5,6], contaminants [3,7], introduced species [8,9], and climate change [10,11]. Changes in freshwater flow, natural or otherwise, can perturb the estuarine salinity gradient [12,13], and influence the abundance and distribution of estuarine organisms through a variety of flow-mediated processes (e.g., [9,14,15]). Moreover, estuarine food webs can be sensitive to small variations in the freshwater input [16,17]. Thus, understanding how freshwater flow controls species abundance, their distribution and interactions requires explicit consideration of spatial and temporal scales and abiotic and biotic habitat components. Impaired freshwater flows can diminish the amount of habitat and the habitat quality for estuarine-dependent species [18,19]. Effects of impaired flows on fish populations include reductions in growth [20,21], survival [21], abundance [10,21], and biomass [22]. This highlights the variety of ecological roles freshwater flows may have on the population dynamics and sustainability of estuarine-dependent species.
The identification of specific ecological mechanisms linking species responses to freshwater flow in estuaries remains challenging due to the highly complex nature of such ecosystems [17,23]. In particular, there is a pressing need for analytical tools to evaluate ecological hypotheses in such ecosystems, and to inform adaptive management of freshwater flows to conserve fish species at risk of extinction (e.g., [9,24]). Because species interactions are notoriously difficult to quantify [25], qualitative analysis of the community matrix is particularly suited to model ecological interactions since essential community properties can be determined by ecosystem topology and sign structure of the interactions, rather than by the values of variables and parameters themselves [26][27][28]. Analysis of qualitative interactions provides a mathematically rigorous foundation for understanding ecosystem behavior [29][30][31] and has been used in a wide range of aquatic ecosystems to evaluate how community stability and community variables respond to abiotic and biotic perturbations, including, nutrient input in freshwater food webs (e.g., [26,32]); introduced species in freshwater communities (e.g., [30,33]) and environmental disturbance in kelp-urchin communities [34]. Qualitatively models in transitional waters have included evaluations of community responses to dumping of mine wastes in shallow-coastal communities (e.g., [35,36]); and increased ocean acidification in an estuarine community [37]. However, there is a need for further metrics to interpret community responses and to determine whether the predicted influence of a perturbation on a given community variable or species of interest is consistent with field observations.
The upper San Francisco Estuary (hereafter upper SF Estuary; California; Figure 1) is an imperiled ecosystem and the central hub of California's water supply infrastructure [38]. Outflow is the most ecologically important source of freshwater to the SF Estuary and represents the net flow at Chipps Island in the Suisun Region (river km 74; Figure 1) after subtracting water diversions and depletions in the upper SF Estuary and its watershed [6,9]. An ecologically meaningful metric of outflow forcing on the upstream distribution of the salinity field in the SF Estuary is X2, which is defined as the distance (km) from the mouth of the estuary (Golden Gate Bridge), to upstream locations where the near-bottom tidally averaged salinity is 2 psu [39]. The upper SF Estuary is home to several threatened and endangered species including the endemic delta smelt Hypomesus transpacificus, an osmerid particularly sensitive to human impacts due to its low fecundity, predominant annual life cycle, and limited spatial distribution [38,40]. Spawning of delta smelt occurs mostly in freshwater during spring and their planktonic larvae tend to gradually move downstream into the low salinity zone (hereafter LSZ, 1-6 psu), [24], where juveniles and subadults generally rear during summer and fall. Despite major long-term changes in pelagic food webs in the upper SF Estuary, calanoid copepods remain one of their most important prey [41]. Delta smelt was listed as threatened (federal and state) in 1993 and up-listed to endangered (state) in 2010 [42].
Understanding the ecological processes underlying the response of delta smelt to fall outflows is relevant to adaptive management in the upper SF Estuary [24,43], and such management actions are intended to improve fall habitat of delta smelt when preceding precipitation and runoff conditions result in high outflows [43]. Outflow controls the location and range of the LSZ (Figure 1), the extent of the abiotic habitat available for subadult delta smelt [10,44], and the position of X2 which controls the upstream distribution of subadult delta smelt during the fall [45]. Nevertheless, because X2 has not been found to be a simple predictor of delta smelt abundance, the need to investigate more complex mechanisms or additional variables has been suggested [39,44]. The high outflows observed in the upper SF Estuary during 2011 prompted monitoring of abiotic and biotic factors, which along with long-term monitoring in the SF Estuary, enabled the performance of interannual comparisons of abiotic and community responses to changes in the position of fall X2 [24]. Although several of the predicted abiotic and biotic responses in that study could not be confirmed from available field data, the abundance of delta smelt in 2011 and early 2012 was one the highest reported since the early 2000s. Yet, delta smelt became very rare during the 2012-2016 drought and failed to increase its abundance during the wet year 2017, indicating the need to evaluate new management options (e.g., [7,9,46]).
Given the complex abiotic and biotic responses to fall outflow in the upper SF estuary [24,38], there is a pressing need for further analytical approaches to evaluate how outflow controls the essential community structure and functional group interactions influencing delta smelt in the LSZ under particular outflow scenarios (hereafter delta smelt subsystems), and to investigate the likely mechanisms controlling the population response of delta smelt in different subsystems ( Figure 1).
Water 2019, 11, x FOR PEER REVIEW 3 of 27 essential community structure and functional group interactions influencing delta smelt in the LSZ under particular outflow scenarios (hereafter delta smelt subsystems), and to investigate the likely mechanisms controlling the population response of delta smelt in different subsystems ( Figure 1). The objectives of this study were to evaluate: (1) the stability of the delta smelt subsystems in the LSZ when fall outflows are alternatively maintained at low-, mid-, and high-X2 positions corresponding to 74, 81, and 85 km, (2) the predicted direction of change and its uncertainty for the abundance of delta smelt and other species and trophic levels in each subsystem under four outflow input scenarios, and (3) whether field patterns of relative abundance for delta smelt and X2 are consistent with the predictions of community models. Three community metrics were adapted to model the response of delta smelt and other community variables to outflows, first, the weighted prediction matrix [47] was modified to denote both the direction of change and its uncertainty in a signed weighted prediction matrix (Ws), second, a modified community matrix specifying press inputs ( ) was used to estimate the combined effect of outflow press perturbations on two or more community variables, and third, the proportional net change of adjoints in the negative community matrix (Δp) was used as a measure of the determinacy in the direction of change of variables in quantitative simulations. Finally, statistical relations between the field relative abundance of subadult The objectives of this study were to evaluate: (1) the stability of the delta smelt subsystems in the LSZ when fall outflows are alternatively maintained at low-, mid-, and high-X2 positions corresponding to 74, 81, and 85 km, (2) the predicted direction of change and its uncertainty for the abundance of delta smelt and other species and trophic levels in each subsystem under four outflow input scenarios, and (3) whether field patterns of relative abundance for delta smelt and X2 are consistent with the predictions of community models. Three community metrics were adapted to model the response of delta smelt and other community variables to outflows, first, the weighted prediction matrix [47] was modified to denote both the direction of change and its uncertainty in a signed weighted prediction matrix (Ws), second, a modified community matrix specifying press inputs (A P ) was used to estimate the combined effect of outflow press perturbations on two or more community variables, and third, the proportional net change of adjoints in the negative community matrix (∆p) was used as a measure of the determinacy in the direction of change of variables in quantitative simulations. Finally, statistical relations between the field relative abundance of subadult delta smelt and X2 enabled to evaluate the predicted population response of delta smelt in each of the modeled subsystems.

Study Area
The SF Estuary is the largest estuary in the west coast of the United States (Figure 1), and has been profoundly altered physically, chemically, and biologically [6,8]. Multiple stressors in the upper SF Estuary have been implicated in the pelagic organism decline occurring in the early 2000s [48,49]. Reservoirs upstream of the SF Estuary, and substantial water diversions within an extensive network of interconnected channels surrounding leveed agricultural islands (the Delta) have greatly altered the hydrodynamics of this estuary [13,50]. Two dataset are central to the evaluation of hydrologic conditions in the SF Estuary: (1) Dayflow, a 2-D mass balance model of freshwater outflows from the tributaries of the Upper SF Estuary [51], and (2) the California Data Exchange Center, a query-based portal which maintains, and operates an extensive hydrologic data collection network [52]. In addition, the Interagency Ecological Program (IEP) provides long-term monitoring data of pelagic, planktonic, and benthic organisms in the SF Estuary [53]. The field data generated by the IEP have been collected consistently over time at multiple locations representing the study areas in the upper SF Estuary (e.g., [24,38,54]). The temporal and spatial coverage of these field surveys has enabled to evaluate the role of hydrodynamic forcing on the distribution and relative abundance of planktonic, pelagic, and benthic species (e.g., [24,38,[55][56][57]). Additional monitoring of aquatic macrophytes by the University of California at Davis, has been conducted extensively using a combination of field sampling and remote sensing in the upper SF Estuary (e.g., [58,59]).
The LSZ area varies greatly with outflow, and is about 90 km 2 when X2 is positioned at 74 km, 50 km 2 when X2 is positioned at 81 km, and 40 km 2 when X2 is positioned at 85 km [24]. Altered fresh water outflow in the upper SF Estuary has greatly changed the abiotic habitat for pelagic fishes [9,54], and facilitated the introduction of species [60,61]. A high proportion of pelagic and planktonic organisms in the upper SF Estuary are also lost due to prevailing upstream flows in the south Delta caused by massive water pumps [5,6]. As is the case of the lower SF Estuary, most communities in the upper SF Estuary are dominated by introduced species [61,62]. The introduced macrophyte Egeria densa was first reported in the Delta in 1946 [63] and became the most abundant submersed plant in the Delta, where it has greatly altered community structure and functions [58,64,65]. The introduction of the clam Potamocorbula amurensis in the late 1980s caused major trophic changes in the benthic and pelagic communities of the SF Estuary [60,66], greatly adding to the filter feeding activity of the clam Corbicula fluminea, an earlier introduction in the upper SF Estuary [67]. In addition, the cyanobacterium Microcystis aeruginosa was first observed in the Delta in 1999 [68] and produces toxic algal blooms that can be harmful to upper trophic levels [69,70].

Community Matrix Models
The community structure and interactions in this study were based on qualitative analysis of the community matrix [29,71]. Given a classical dynamical ecosystem of n interacting variables (N i ), the change in the equilibrium level of species n i * is dn * i dt = f i(n 1 , n 2 , . . . , n n ; c 1 , c 2 , . . . , c n ) (i = 1, n).
where the n i represent species abundances or abiotic factors. The c i are parameters that govern biological rates (e.g., birth, mortality, growth), and depend on the species and environment [72,73]. Modeled subsystems were represented by a signed digraphs, where element a ij of a qualitatively defined community matrix (A o ) of order n denotes the net effect to variable i from variable j (for i and j from 1 to n). For example, a subsystem composed of three trophic levels can be represented by three variables and their interactions ( Figure 2). Three possible qualitative effects to variable i from variable j were assigned to each matrix element: a i,j = 1, a i,j = −1 and a i,j = 0, which in the neighborhood of any equilibrium point respectively denote positive, negative, and no effect on the instantaneous growth of variable i due to increase in the level of variable j. The stability of the community matrix was evaluated using the characteristic polynomial: p(λ) = |A o − λI| = a 0 λ n + a 1 λ n − 1 + a 2 λ n − 2 + . . . a n = 0 (2) where, A o is as defined earlier, I is the n-dimensional identity matrix, and λ are the roots or eigenvalues [71]. Following a perturbation, negative real roots means the community can return to equilibrium. Conversely, positive real roots indicate the ecosystem moves away from equilibrium, as in unstable ecosystems. Multiple zero roots can result in unstable equilibrium [73], and zero overall feedback indicates neutral stability [74].
Water 2019, 11, x FOR PEER REVIEW 5 of 27 of variable i due to increase in the level of variable j. The stability of the community matrix was evaluated using the characteristic polynomial: where, A o is as defined earlier, I is the n-dimensional identity matrix, and  are the roots or eigenvalues [71]. Following a perturbation, negative real roots means the community can return to equilibrium. Conversely, positive real roots indicate the ecosystem moves away from equilibrium, as in unstable ecosystems. Multiple zero roots can result in unstable equilibrium [73], and zero overall feedback indicates neutral stability [74]. showing qualitative interaction terms for 3 variables such as species or trophic groups, with diagonal matrix elements −a11, −a22, and −a33 representing negative self-effects of each variable and zero denoting no effect, and (B) its corresponding signed-digraph. Open circles denote variables (Ni), with lines between variables ending in arrows and bubbles respectively denoting positive and negative effects, and curved lines originating and ending in the same variable denoting self-effects, after [29].

Community Structure and Outflow Scenarios
Qualitative analysis was used to model three flow-dependent fall community structures in the upper SF Estuary under high-, mid-, and low-net Delta outflows corresponding respectively to X2 positions of 74, 81, and 85 km [24]. The three essential components of qualitative models (community variables, interactions or links, and driving force or perturbation) were based on the following five steps: (1) Compilation of long-term data and studies in the upper SF Estuary, including: monitoring data for delta smelt and other species [53], outflow and X2 data [51], interpretation of the spatiotemporal distribution of the salinity field in the LSZ based on the three-dimensional UnTRIM hydrodynamic model for the upper SF Estuary [24,75], and baseline knowledge of the modeled ecosystem (e.g., [6,10,39,40,54,55,63,76,77]). (2) Consideration of dynamic and stationary factors for developing conceptual models of estuarine communities [78], where dynamic factors include physico-chemical and biological characteristics of the low salinity habitat corresponding to the LSZ each X2 position, and stationary factors include geographically fixed habitat features at each X2 position such as substrate, erodible sediment, and bathymetry in the low salinity habitat [24,38]. (3) Refinement of conceptual models into different subsystems describing the essential community variables influencing subadult delta smelt based on ecological syntheses of long-term field data and studies [24,38]. Community variables selected for each of the three modeled subsystems were based on functional groups (e.g., [74,79]), and their predominant spatio-temporal overlap with dynamic and stationary abiotic and biotic factors. Except for delta smelt and species with significant ecological impacts (the cyanobacterium M. aeruginosa [80], the clams P. amurensis and Figure 2. Example of (A) a 3 × 3 community matrix (A o ) showing qualitative interaction terms for 3 variables such as species or trophic groups, with diagonal matrix elements −a 11 , −a 22 , and −a 33 representing negative self-effects of each variable and zero denoting no effect, and (B) its corresponding signed-digraph. Open circles denote variables (N i ), with lines between variables ending in arrows and bubbles respectively denoting positive and negative effects, and curved lines originating and ending in the same variable denoting self-effects, after [29].

Community Structure and Outflow Scenarios
Qualitative analysis was used to model three flow-dependent fall community structures in the upper SF Estuary under high-, mid-, and low-net Delta outflows corresponding respectively to X2 positions of 74, 81, and 85 km [24]. The three essential components of qualitative models (community variables, interactions or links, and driving force or perturbation) were based on the following five steps: (1) Compilation of long-term data and studies in the upper SF Estuary, including: monitoring data for delta smelt and other species [53], outflow and X2 data [51], interpretation of the spatio-temporal distribution of the salinity field in the LSZ based on the three-dimensional UnTRIM hydrodynamic model for the upper SF Estuary [24,75], and baseline knowledge of the modeled ecosystem (e.g., [6,10,39,40,54,55,63,76,77]). (2) Consideration of dynamic and stationary factors for developing conceptual models of estuarine communities [78], where dynamic factors include physico-chemical and biological characteristics of the low salinity habitat corresponding to the LSZ each X2 position, and stationary factors include geographically fixed habitat features at each X2 position such as substrate, erodible sediment, and bathymetry in the low salinity habitat [24,38]. (3) Refinement of conceptual models into different subsystems describing the essential community variables influencing subadult delta smelt based on ecological syntheses of long-term field data and studies [24,38]. Community variables selected for each of the three modeled subsystems were based on functional groups (e.g., [74,79]), and their predominant spatio-temporal overlap with dynamic and stationary abiotic and biotic factors. Except for delta smelt and species with significant ecological impacts (the cyanobacterium M. aeruginosa [80], the clams P. amurensis and C. fluminea [56], and the macrophyte E. densa [58], Table 1), other functional groups included trophic levels to minimize redundant species interactions (e.g., [74,81]), (Table 1). 323 [24] 227 [24] 142 [24] Major abiotic forcing factor controlling the overlap between dynamic and stationary factors in the upper San Francisco Estuary, the position and size of the low salinity zone and the low salinity habitat of delta smelt and community interactions [24,38] (4) Reformulation of conceptual models into signed digraphs based on qualitative model guidelines (e.g., [71,79,93]). Negative self-effect (self-damping) was assumed to arise for each variable from density-dependent growth rate or a limited source, as in the case of nutrients [71]. Each community variable was then implicitly connected to other variables or abiotic factors through negative feedback [71,74]. Reported community interactions considered for the modeled subsystem included: predation (+, −); interference competition (−, −); and amensalism (0, −), (Appendix A). (5) Estimation of the direction of change of community variables (+, 0, −) in response to increased outflow (Table 1). Four outflow input scenarios were modeled, with the first scenario accounting for the effect of outflow on the previously referred species having significant ecological impacts. The outflow inputs for the three subsequent scenarios were used to evaluate whether cumulative outflow inputs in each subsystem could reinforce or reverse potential responses on delta smelt and other community variables. These scenarios included: scenario 1 + phytoplankton (scenario 2), scenario 2 + zooplankton (scenario 3), and scenario 3 + delta smelt (scenario 4).

Qualitative Analyses
The modeled subsystem at each X2 position was considered stable if the two Routh-Hurwitz (R-H) criteria [94] were met. The first R-H criterion requires all roots of the characteristic polynomial be different from zero and have the same sign. The second R-H criterion requires all Hurwitz determinants be positive, resulting in stronger feedback at lower feedback levels than at higher levels. The effect of sustained input (C h ), on the rate of change of a hth variable corresponding to the jth variable in the community matrix was predicted assuming an ecosystem in moving equilibrium [29,71]. Such input could result in positive, negative, or no effect on the equilibrium magnitude or carrying capacity (N i * ) of a species. The inverse of the negative community matrix (−A ij −1 ) provides an estimate of change in the equilibrium level of variable N i * carrying capacity due to change in parameter C h [95].
where −C ji T is the adjoint of the negative community matrix (Adj − A ij ), and |−A ij | is the determinant of the negative community matrix. Since |−A ij | is constant and positive for each element of −A −1 in stable ecosystems, the predicted response under local equilibrium was inferred from the Adj − A: which corresponds to the transpose of the negative cofactor matrix [33]. The absolute number of positive and negative complementary feedback cycles (Tc) which contributes to each element of the Adj − A was used to derive model uncertainty. Because Tc differed from zero in all present models, it was deemed practical to modify the weighted prediction matrix [47] to express both the direction of change in response to input and its uncertainty by defining a signed weighted prediction matrix (Ws): where Adj − A is as defined in Equation (4), and "→" is a vectorized matrix operator denoting element-by-element division for each element of the Ws matrix [47], and Tc is as defined earlier and 0. Hence, both the lowest Ws value (−1, negative response) and highest Ws value (+1, positive response) indicate unconditional sign determinacy and Ws = 0 denote total uncertainty in the direction of change, while absolute values of Ws ≥ 0.5 and < 1 are indicative of high level of sign determinacy, after [27].
To estimate the overall influence of a press input (P k ) due to a disturbance k acting directly upon two or more variables, the combined effect of P k on the direction of change for community variable i and their corresponding uncertainty level was computed from the corresponding matrix element Ws ik (Equation (5)) and these were derived from a modified community matrix specifying press inputs (A P ): where a 11 to a nn represent the interaction elements a ij of community matrix A o , and P ik denotes the press input k on community variable i, which in a qualitatively defined ecosystem is specified as having positive (1), negative (−1) or no effect (0). The last row of matrix A P includes zero elements to denote no effect of community variables on P k , and a negative element (−1) to preserve the stability of matrix A P provided the community matrix A o is stable. Qualitative analyses were conducted in Maple 18 [47].

Quantitative Simulations
Simulations using A P further enabled to quantitatively assess the stability of the community and the direction of change of its variables relative to the qualitatively specified matrix with same structure and sign of interactions. Simulations were performed by assigning interaction strengths between 0.01 and 1.0 to non-zero elements using a pseudorandom number generator from a uniform distribution [31]. The overall direction of change and determinacy in a quantitatively-specified ecosystem was derived from the proportional net change for each community variable (∆p): where f (+Adj − A ik )* and f (−Adj − A ik )* respectively denote the absolute frequencies of positive and negative Adj − A ik in stable matrices, and N is the total number of simulations for each perturbation scenario, including stable and non-stable matrices (N = 10,000). As in the case of Ws, the value of ∆p provided both the overall direction of change and its uncertainty for each community variable, with extreme potential values (∆p = −1, negative response and ∆p = 1, positive response) both indicating no uncertainty in the direction of change among simulations, ∆p ≥ |0.5| and < |1|indicating low uncertainty, and ∆p = 0 denoting total uncertainty. Quantitative simulations were performed in a program written in C# version 3.0 [96].

Statistical Analyses
To determine the extent to which the predicted response to outflow input by each variable in qualitative models represented the overall direction of change in quantitative simulations, ∆p was plotted vs. Adj − A for all scenarios at each X2 position. Values of ∆p were also plotted vs. Ws to compare the direction of change and uncertainty between qualitative and quantitative predictions. These variables were fitted at each X2 position using Lowess smoothing. The predicted population response for subadult delta smelt by community models at each X2 position was compared with the relative abundance for subadult delta smelt as a function of X2 to evaluate model predictions. Relative abundance data for subadult delta smelt, monitored since 1967, was obtained from the California Department of Fish and Wildlife's Fall midwater trawl index for the period 1967-2017 [97]. The average position of X2 during September-October was derived from the California Department of Water Resources Dayflow program [51]. Linear regression analysis was used to evaluate the slope of the relative abundance for delta smelt during fall (Y) versus the average X2 position in September-October (X). Regressions were conducted for three ranges of X2, which included the three X2 positions in community models (X2 ≥ 74, ≥81, and ≥85 km). The entire historical record of subadult delta smelt abundance as a function of the average September-October X2, was further considered as a baseline relative to the three X2 positions evaluated in the present community models and a previous fall outflow evaluation [24]. Data were analyzed and plotted using Curve Expert Professional (version 2.6; D.G. Hyams).

Community Stability
Compared to the models at the low X2 position, the modeled communities at the mid-and high-X2 positions included a higher number of introduced species and more non-trophic interactions (Figure 3, Appendix A). Although all qualitative community models were stable, the percent of stable simulations decreased greatly from downstream (low X2) to upstream (high X2) positions in the LSZ, and such decrease in stability seemed primarily due to the increasing percent of models with stronger high-level feedbacks than low-level feedbacks at the mid-and high-X2 positions ( Table 2).

Community Response to Outflow
The predicted change in abundance for delta smelt and other community variables to outflow input varied greatly among the three X2 positions (Figure 4), both in the case of qualitative models ( Figure 5) and simulations ( Figure 6). In response to outflow input for the delta smelt subsystem at the low X2 position (Figure 4), the predicted direction of change in the abundance of all community variables was positive, except for the clam P. amurensis, with predicted responses of species and trophic levels generally showing high sign determinacy, particularly for delta smelt, its predators, and zooplankton across most outflow input scenarios (Figures 5 and 6). Outflow input at the mid X2 position (Figure 4), revealed the cyanobacterium M. aeruginosa was the only community variable with a consistently negative response and high sign determinacy in both qualitative models and simulations. Yet, the clam C. fluminea and predators of delta smelt showed high sign determinacy under most quantitative simulation scenarios at the mid X2 position (Figures 5 and 6).
The cyanobacterium M. aeruginosa suggested a negative response to outflow at the high X2 position (Figure 4), but to a lower extent compared to the mid X2 position (Figures 5 and 6). The Brazilian waterweed, E. densa was the only other variable at the high X2 position indicating high degree of qualitative sign determinacy but only at one outflow scenario ( Figure 5). In contrast, simulations at the high X2 position indicated all the predicted responses of community variables had low sign determinacy ( Figure 6).
The relation between ∆p and Adj − A indicated the predicted direction of change for species and trophic groups in response to outflow input was generally consistent between qualitatively and quantitatively specified subsystems at each X2 position; this despite a plateau in ∆p was towards its highest values ( Figure 7A). As inferred from ∆p and Ws, the direction of change and degree of certainty in the direction of change in quantitatively specified simulations were associated with the corresponding qualitative model, with such relations appearing progressively more sigmoidal toward higher-X2 positions ( Figure 7B). The cyanobacterium M. aeruginosa suggested a negative response to outflow at the high X2 position (Figure 4), but to a lower extent compared to the mid X2 position (Figures 5 and 6). The Brazilian waterweed, E. densa was the only other variable at the high X2 position indicating high degree of qualitative sign determinacy but only at one outflow scenario ( Figure 5). In contrast, simulations at the high X2 position indicated all the predicted responses of community variables had low sign determinacy ( Figure 6).     Figure 3) and legend denotes the four outflow input scenarios (Figure 4). Values of │Δp│ can range between 1 (unconditional sign determinacy) and 0 (complete uncertainty), with │Δp│ between 0.5 (dotted lines) and <1 implying low uncertainty in the predicted direction of change.
The relation between Δp̂ and Adj − A indicated the predicted direction of change for species and trophic groups in response to outflow input was generally consistent between qualitatively and quantitatively specified subsystems at each X2 position; this despite a plateau in Δp̂ was towards its highest values ( Figure 7A). As inferred from Δp̂ and Ws, the direction of change and degree of certainty in the direction of change in quantitatively specified simulations were associated with the corresponding qualitative model, with such relations appearing progressively more sigmoidal toward higher-X2 positions ( Figure 7B).

Delta Smelt Abundance and X2
When only considering years in which the average position of X2 in September-October was located at 81 km and further upstream, the relative abundance of subadult delta smelt was not associated to the position of X2 (Figure 8). However, when all years in which X2 was located at 74 km or further upstream were included, the relative abundance of subadult delta smelt was inversely associated to the average position of X2 in September and October (Figure 8). Moreover, when considering all the years with available relative abundance for subadult delta smelt over the period 1967-2017, an inverse association between the relative abundance of subadult delta smelt and the average position of X2 during September-October was also detected, and such relation included X2 positions as low as 65 km (Figure 8).

Delta Smelt Abundance and X2
When only considering years in which the average position of X2 in September-October was located at 81 km and further upstream, the relative abundance of subadult delta smelt was not associated to the position of X2 (Figure 8). However, when all years in which X2 was located at 74 km or further upstream were included, the relative abundance of subadult delta smelt was inversely associated to the average position of X2 in September and October (Figure 8). Moreover, when considering all the years with available relative abundance for subadult delta smelt over the period 1967-2017, an inverse association between the relative abundance of subadult delta smelt and the average position of X2 during September-October was also detected, and such relation included X2 positions as low as 65 km (Figure 8).

Community Stability
While the modeled subsystems in the upper SF Estuary at the low-, mid-, and high-X2 positions were qualitative stable based on the R-H criteria, quantitative simulations using such stability criteria revealed a significantly higher percent of stable models at the low X2 position, which could promote enhanced ecosystem resilience when the LSZ is predominantly located in the Suisun region ( Figure  1). In contrast, a decreasing number of simulations met both R-H stability criteria at the mid-and high-X2 positions, respectively (Table 2). Such spatial stability pattern supports the hypothesis that fall outflows maintaining the X2 position at 74 km are more conducive to sustaining a community structure supporting the delta smelt population compared to upstream X2 positions. The stability patterns inferred from simulations are consistent with marked reductions in upstream fall habitat quality for subadult delta smelt and with the low abundance for this species since the early 2000s [24,54]. These changes have implicated long-term stressors in the upper SF Estuary, including water diversions, introduced species and contaminants [6,13,50]. As in the case of other ecosystems where native estuarine fish species evolved in turbid waters [98], the decreased turbidity in the Delta is partly due to the spread of E. densa and other introduced macrophytes, which have favored the population expansion of introduced predators [38,99]. Although herbicides have been used in an attempt to control macrophytes in the Delta, the extent to which such efforts have been effective is unclear. Nonetheless, some transitional waters formerly dominated by macrophytes seem to have shifted to pelagic ecosystems through the long-term use of herbicides [100].

Community Stability
While the modeled subsystems in the upper SF Estuary at the low-, mid-, and high-X2 positions were qualitative stable based on the R-H criteria, quantitative simulations using such stability criteria revealed a significantly higher percent of stable models at the low X2 position, which could promote enhanced ecosystem resilience when the LSZ is predominantly located in the Suisun region ( Figure 1). In contrast, a decreasing number of simulations met both R-H stability criteria at the mid-and high-X2 positions, respectively (Table 2). Such spatial stability pattern supports the hypothesis that fall outflows maintaining the X2 position at 74 km are more conducive to sustaining a community structure supporting the delta smelt population compared to upstream X2 positions. The stability patterns inferred from simulations are consistent with marked reductions in upstream fall habitat quality for subadult delta smelt and with the low abundance for this species since the early 2000s [24,54]. These changes have implicated long-term stressors in the upper SF Estuary, including water diversions, introduced species and contaminants [6,13,50]. As in the case of other ecosystems where native estuarine fish species evolved in turbid waters [98], the decreased turbidity in the Delta is partly due to the spread of E. densa and other introduced macrophytes, which have favored the population expansion of introduced predators [38,99]. Although herbicides have been used in an attempt to control macrophytes in the Delta, the extent to which such efforts have been effective is unclear.
Nonetheless, some transitional waters formerly dominated by macrophytes seem to have shifted to pelagic ecosystems through the long-term use of herbicides [100].

Community Model Predictions
Comparisons between predicted population responses of community variables considered in this modeling study and field data [24] were possible for phytoplankton, zooplankton, P. amurensis, and delta smelt. Consistent with the direction of change and sign determinacy for delta smelt predictions at the low X2 position (Figures 4-6), the relative abundance of subadult delta smelt in fall 2011 when X2 was at 74 km was a ten-fold higher compared to years with higher fall X2 positions since the decline of pelagic species in early 2000s [38]. As in the case of the stability for quantitatively simulated subsystems, the response of delta smelt to sustained outflow was influenced by the location of the LSZ and hence, the underlying community structure. When outflow caused X2 to reach 74 km in the Suisun region, a positive population response was suggested for delta smelt and such response had high sign determinacy despite the ambiguous response of P. amurensis under most outflow input scenarios.
Field data for phytoplankton in fall 2011 showed either similar or greater average phytoplankton biomass, as measured by chlorophyll-α, at the low-X2 position relative to the mid-and high-X2 positions [24]. Such field pattern could result from variability in the mass balance of phytoplankton production, total grazing, and transport [89,101]. A positive overall phytoplankton response to outflow input was only predicted for the models at the low-X2 position, and such prediction had low uncertainty or marginally low uncertainty for scenarios 2 and 4 (Figures 4-6). Hence, suggesting a conditional phytoplankton response to outflow.
Despite the high variability in zooplankton biomass contributed by calanoids-a key prey item of delta smelt-such biomass seemed slightly higher at the low X2 position in fall 2011 compared to years in which the LSZ was found at higher X2 positions [24]. A positive response of zooplankton was predicted in the LSZ at the three X2 positions but such predictions only had low uncertainty at the low-X2 position under most outflow scenarios (Figures 4-6). Importantly, considering the larger area of the LSZ when X2 is located at 74 km relative to higher X2 positions (Figure 1), the potential carrying capacity for species and trophic levels could be significantly higher when the LSZ overlaps Suisun Bay in the fall.
The field abundance of P. amurensis in the Suisun Bay region was lower during October 2011 when X2 was at 75 km compared to the same location in the lower-outflow year 2010, when fall X2 was at 85 km, which is consistent with a negative local response of this species to outflow ( [56], their Table 1), and with the low uncertainty in the direction of change under scenario 1 for the low-X2 position (Figures 4-6). In contrast, the uncertain direction of change of P. amurensis under other outflow input scenarios at the low X2 position implies the negative direct effect of outflow on P. amurensis would be offset by potential outflow-induced increases in its food supply. Because the approximate center of distribution and abundance of the clam P. amurensis is Suisun Bay and the downstream embayment, San Pablo Bay ( [56,67], their Figure 3) its abundance is likely higher when the LSZ overlaps Suisun Bay relative to LSZ locations further upstream. Moreover, changes in its benthic abundance would not be as responsive to changes in the salinity field compared to pelagic species.

Delta Smelt Abundance under Different Ranges of X2
The predicted response for delta smelt in community models was generally consistent with the patterns of relative abundance for subadult delta smelt across different ranges of X2 (Figure 8). With the exception of 2011 and 2017, no years since 2000 have shown the average of X2 in September-October located in the Suisun region. However, as in the case of year 2017, low values of X2 in fall since 1967 have not always been associated with high relative abundance for subadult delta smelt (Figure 8). This suggest that additional abiotic factors can be important in controlling the abundance of delta smelt in some years. Unlike 2011, the summer 2017 was warmer than usual, which could have reduced the survival for juvenile delta smelt (e.g., [11,38]). Importantly, since year 2000 there has been a substantially lower occurrence of years where the average position of X2 in September-October was ≤75 km in the upper SF Estuary (11.1%) compared to the period 1967-1999 (39.4%). A contributing factor to such difference is the long-term increase in freshwater diversions, both upstream of the SF Estuary and in the south Delta. For example, total water diversions from 1986 to 2005 held X2 upstream of 71 km nearly 80% of the time. In contrast, assuming no upstream water diversions, X2 would have been maintained upstream of 71 km only 50% of the time [13].
The predicted response for delta smelt in community models was generally consistent with the reported decline of abiotic habitat quality indices of at higher positions of X2 (e.g., [10,54]). The present study suggests a key role of community structure in the LSZ given its influence on community stability and response to outflow. Hence, supporting the hypothesis that fall outflow exerts a positive influence on the abundance of subadult delta smelt when the X2 is positioned at 74 km [24]. When considering the four outflow scenarios for each of the three X2 positions, a generally similar response of community variables was apparent across outflow input scenarios, and this was suggested despite the overall ambiguity of subsystems at the mid-and high-X2 positions (Figures 5 and 6). This implies the propagation of outflow inputs through multiple feedbacks in the community could make it difficult to distinguish specific outflow input mechanisms among years.

Community Models and Prediction Metrics
Because qualitative models sacrifice precision for generality and realism [71], the qualitative models used in this study were intended to provide insights on the community stability and the complex mechanisms mediating the response of delta smelt and its community to fall outflows, rather than quantifying changes in the abundance of community variables in the upper SF Estuary. Given the large differences in qualitatively derived stability among X2 positions, the use of ∆p in quantitatively specified subsystems seemed a useful community metric to evaluate both the overall direction of change of community variables and the determinacy in the overall direction of change across subsystems with the same topology but varying interaction strengths. Although estimating the direction of change and the degree of uncertainty in the response of community variables to multiple inputs entails several steps (e.g., determining qualitative direction of change of community variables, the absolute feedback matrix and the weighted predictions from each element of the Adj − A, [47]), the present study showed that the response of community variables to multiple inputs can be alternatively obtained by defining the press inputs for each variable in a modified community matrix (A P ). This enabled readily estimation of the qualitative response of each community variable, and their corresponding weighted predictions. In the case of communities without zero absolute feedback matrix elements, it is also feasible to define signed weighted predictions to facilitate comparisons between qualitative models and quantitative simulations. Moreover, the use of matrix A P can provide independent verification of detailed computations involving inputs to several community variables. For quantitatively specified communities, the use of matrix A P in combination with ∆p seemed particularly useful to estimate the determinacy in the direction of change for each community variable under simulated press inputs and interaction strengths among community variables. Although the signed weighted feedback seemed a useful measure of uncertainty for qualitatively stable models, and was significantly related to the estimated uncertainty of quantitatively specified models ( Figure 7B), quantitative simulations showed that interaction strength can greatly influence the stability of community models with the same sign structure of stable qualitative models. This in turn could influence the determinacy of the predicted direction of change following a perturbation. However, because the determinacy level indicated by ∆p accounted for the percent of non-stable models in simulation analyses, the use of ∆p can complement qualitative analyses of the community matrix and help to evaluate hypotheses about community stability and their response to perturbations.
Despite the differences in objectives and methods in past modeling studies related to pelagic fishes and their communities in Upper SF Estuary, independent analyses of flow-related metrics (e.g., X2, salinity, specific conductance, water diversions) have revealed that hydrological conditions play a major ecological role through their interactions with multiple natural and anthropogenic forces (e.g., [6,9,102]). For example, generalized additive models (GAM) have shown association between spatial and temporal patterns in specific conductance and turbidity and long-term declines in habitat quality for delta smelt and other pelagic fishes (e.g., [54,99]). Additional GAM-based habitat suitability scenarios for delta smelt suggested changes in outflow due to future development and climate change could lead to increased habitat declines for this species across a broad range of hydrological conditions [10]. A Bayesian-based multivariate autoregressive model (MAR) revealed increased spring and fall X2 could be linked to the decline of several pelagic fishes and their zooplankton prey since the early 2000s [55]. That study further suggested negative influence of spring and winter south Delta diversions and high summer water temperature on delta smelt. Another MAR model evaluating mechanisms of species invasions and their impacts in the Upper SF Estuary suggested salinity intrusion was a primary pressure facilitating species invasions and emphasized a central role introduced species in altering multi-trophic interactions [57]. Individual-based models for delta smelt have shown a combination of abiotic and biotic factors could control its population, including hydrodynamics, water temperature, zooplankton density, long-term changes in the prey-base, entrainment losses due to water diversions in the south Delta, and salinity (e.g., [103][104][105]). Moreover, hydrodynamic and statistical modeling showed the decline in phytoplankton in the upper SF Estuary could be largely explained by the combined impacts of increased freshwater diversions from the south Delta since the late 1960s, and the filter-feeding activity of the clam P. amurensis [106]. Hydrodynamic and statistical modeling also revealed that winter-spring diversions altered the salinity habitat for several species of pelagic fishes, and these changes were significantly associated to their corresponding interannual population trends [9]. Thus, these modeling studies lend further support to the connection between hydrological alteration and multiple processes contributing to the decline of pelagic fishes in the upper SF Estuary [13,38,102].
Given the physical variability of estuarine ecosystems (e.g., [6,13,15]), consideration of the variability in the salinity field seems particularly relevant when evaluating hydrodynamic forcing on estuarine communities. This study is the first qualitative analysis of the community matrix evaluating the role of freshwater flows on aquatic communities, and the first qualitative modeling study in transitional waters based on a hydrological model [75]. Qualitative community models in other estuaries have examined other types of disturbances, for example an intertidal community composed of native and introduced species (Yaquina Bay, Oregon, USA) suggested weak or absent overall feedback, which, to some extent could have buffered such community from impacts of introduced species [74]. In another estuarine community (Willapa Bay, Washington, USA), qualitative networks suggested ocean acidification could stimulate production of phytoplankton and eelgrasses [37]. Yet, the opposite predicted responses for two shellfish species were sensitive to the assumed perturbation on primary producers. This highlights the value of qualitative modeling to identify likely community outcomes while recognizing the uncertainties of complex ecological interactions. As in the case of the SF Estuary, other estuarine ecosystems have been primarily modeled quantitatively. For example, an oxygen-balance model coupled to fish growth models in a range of estuary types in South-Western Australia indicated that climate change could induce declines in river discharge, exacerbating the impacts of reduced rainfall and water diversions on hypoxia in these estuaries [107]. Such increased hypoxia could then magnify the significant reported differences in habitat quality for the growth and body condition of the teleost Acanthopagrus butcheri in those estuaries. A coastal ecosystem model (Integrated Compartment Model) in southern Louisiana, USA [108] suggested that even doubling current riverine flows would be insufficient to fully counteract the impacts of future sea level rise on forested wetlands, implying that habitat expansion for marine fishes would result in habitat contraction for freshwater fishes. Application of the ecosystem model AQUATOX 3.1 in the Minho Estuary, NW coast of Iberian Peninsula [109] revealed that estuarine production was strongly dependent on hydrodynamics, with temperature and river flow showing antagonistic interactions on macroinvertebrate communities. Given the increasing threats to transitional water ecosystems, the integration of qualitative and quantitative models to assess ecosystems responses to multiple perturbations could further promote ecologically sound management in these ecosystems.

Management Implications
Based on the modeled outflow-mediated community responses and the relations between subadult delta smelt abundance and fall X2 inferred from this study, the sustainability of this species could be compromised in part by the high prevalence of low-and mid-fall-outflows. Given the tendency of the delta smelt population to increase in high outflow years, the population could still be resilient provided the limited outflow in most years is not further reduced by the high percent of upstream water diversions [6,13], particularly in non-wet years [9]. While full understanding of the mechanisms by which outflow influences the abundance of a species like delta smelt does not seem necessary for conservation purposes, there is great value in continued experimental and modeling evaluations on the role of hydrological alteration and other perturbations to mitigate cumulative stressors in the upper SF Estuary, including those anticipated from climate change [10,11], provided such evaluations are implemented using sound ecological and adaptive management guidelines (e.g., [110][111][112]). Considering the long-term population declines for delta smelt and other native species like longfin smelt, Spirinchus thaleichthys, and their record low abundances since the drought beginning in 2012 [7,46], the positive influence of outflow in other seasons for delta smelt and other species (e.g., [7,9,38,102]) suggests the need to integrate flow adaptive management actions across species, seasons and hydrological conditions not restricted to wet years.

Conclusions
This study supports the conclusion that fall outflow in the upper SF Estuary tends to increase the abundance of delta smelt when the position of X2 overlaps Chipps Island (X2 = 74 km), but not when X2 is at 81 km or further upstream. The stability patterns of three modeled subsystems of delta smelt in the upper SF estuary support the hypothesis that fall outflows that maintain the X2 position at 74 km are more likely to result in a stable community structure supporting the delta smelt population compared subsystems located at higher X2 positions. Community models suggested a dual role of outflow in controlling the population dynamics delta smelt in the upper SF Estuary, first by determining the geographical location of the LSZ and the underlying community structure, and second by exerting press perturbations on community variables influencing biological rates, as mediated by trophic and non-trophic interactions. Outflow seemed to control community interactions and stability patterns at each X2 position by forcing the overlap among pelagic and benthic species and trophic levels based on their salinity-dependent and geographically-dependent distributions. Integration of flow adaptive management actions across seasons, species, and hydrological conditions could help to address the declining population trends for delta smelt and other native fish populations in the upper SF Estuary. Qualitative community models can complement ecological syntheses, hydrodynamic models and quantitative ecological models to inform ecologically sound management of aquatic ecosystems.

Conflicts of Interest:
The author declare no conflict of interest. Reference to trade names does not imply endorsement by the U.S. Government.