*4.1. 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].

### *4.2. 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-<sup>α</sup>, 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 e ffect 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.

### *4.3. Delta Smelt Abundance under Di*ff*erent 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 di fferent 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 sugges<sup>t</sup> 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 di fference 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 di fficult to distinguish specific outflow input mechanisms among years.

### *4.4. 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 di fferences in qualitatively derived stability among X2 positions, the use of Δ*p*<sup>ˆ</sup> 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<sup>P</sup>*). 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<sup>P</sup>* can provide independent verification of detailed computations involving inputs to several community variables. For quantitatively specified communities, the use of matrix *A<sup>P</sup>* in combination with Δ*p*<sup>ˆ</sup> 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*<sup>ˆ</sup> accounted for the percent of non-stable models in simulation analyses, the use of Δ*p*<sup>ˆ</sup> can complement qualitative analyses of the community matrix and help to evaluate hypotheses about community stability and their response to perturbations.

Despite the di fferences 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–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 managemen<sup>t</sup> in these ecosystems.
