Physicochemical and urban land-use characteristics associated with resistance to precipitation in estuaries vary across scales

Estuaries are subject to frequent stressors, including elevated nutrient loading and extreme hydrologic events, which impact water quality and disrupt ecosystem stability and function. The capacity of an estuary to resist changes in function in response to precipitation events is an important descriptor of ecosystem response dynamics, especially when interpreted within the context of baseline system condition and ecological state. However, the factors related to estuarine responses to extreme precipitation remain poorly constrained. This knowledge gap complicates our ability to identify estuaries that are more likely to undergo major shifts in ecosystem services and predict the effects of urban and precipitation disturbances on estuarine water quality. We investigate which physicochemical and land-use characteristics are associated with ecological resistance to precipitation – defined as the magnitude of ecosystem change induced by an event – in five disparate estuaries distributed across the continental United States. Using long-term meteorological and water quality data from the National Estuarine Research Reserve System along with land use/land cover and population data, we examine relationships between the resistance index – a proxy for ecosystem stability calculated using dissolved oxygen – and physicochemical and urban land use characteristics on local-to-continental scales. Contrary to our initial hypothesis, we found that more urbanized estuaries tended to be more resistant to precipitation events in this dataset, possibly due to persistent disturbances to their baseline dissolved oxygen levels, and that water temperature, water column depth, turbidity, nitrogen, and chlorophyll a showed significant but variable associations with resistance at the continental scale. These continental-scale patterns were modulated by estuarine salinity and varied across individual estuaries, where additional relationships between resistance and salinity, phosphate concentrations, N : P, tree cover, and cropland emerged. Our findings suggest that the relationships between urbanization and estuarine stability are complex and context-dependent, and thus management strategies need to consider both broad generalizations and local conditions. Considering emerging stressors from new environmental scenarios and from urbanization, these results may help improve representation of the impacts of disturbances in large-scale models and inform management decisions regarding estuarine water quality.

1 Introduction

Estuaries are highly dynamic environments that often connect freshwater and saltwater systems, cycle organic matter and nutrients from land to oceans, and provide essential ecosystem services (Bianchi, 2007; He and Silliman, 2019). The function of estuarine ecosystems as unique sites for carbon and nutrient cycling, and habitats for macro- and micro-flora and fauna, relies on stability of a predictable range of dynamic processes like temperature fluctuations, hydrology and nutrient mixing. However, anthropogenic activities and extreme precipitation events threaten estuarine ecosystem stability and function (Kemp et al., 2009; Zhang et al., 2010). Predicted increases in urban population size and in the frequency and intensity of precipitation events highlight the urgency to understand the response of estuaries to urbanization and new precipitation patterns (Kyzar et al., 2021; Li et al., 2019; Martínez et al., 2007; Pickett et al., 2011). Yet, the factors that impact urban estuaries' response to precipitation are not fully understood.

When combined with intense precipitation, watershed urbanization and associated changes in land use/land cover (LULC) often result in increased stream hydrological flashiness (Gannon et al., 2022; Grimm et al., 2008; Reisinger et al., 2017). Hydrological flashiness induces high flow rates that cause changes in channel morphology (Booth and Jackson, 1997; Gregory, 2011; Leopold, 1968; Vietz et al., 2016), habitat destruction (Walsh et al., 2005), and disruption of microbial metabolic processes (Reisinger et al., 2017; Uehlinger, 2000). Flashiness can also drastically affect primary production – a regulatory component of dissolved oxygen (DO) dynamics in aquatic environments – through increases in flow velocity, transport of phytoplankton, sediment migration, and light limitation through increased turbidity (Bernot et al., 2010; Fisher et al., 1982; McSweeney et al., 2017; Reisinger et al., 2017; Uehlinger, 2000).

While DO is dynamic and depends on myriad biological, chemical, and physical processes, it is essential for many estuarine processes and has been widely used as an indicator for overall ecosystem function (Abdul-Aziz et al., 2007; Abdul-Aziz and Gebreslase, 2023; Chapra, 2008; Cox, 2003; Kannel et al., 2007; Zhi et al., 2021). Dissolved oxygen is critical to maintaining the life cycles of macro- and micro-fauna, and supports the biogeochemical cycling of carbon and nutrients by serving as a terminal electron acceptor (Bernhardt et al., 2018; Chapra, 2008; Zarnetske et al., 2012). Urban areas in particular have been linked to altered aquatic DO concentrations and other water quality parameters (Bernhardt et al., 2008; Chang, 2005; Freeman et al., 2019; Vietz et al., 2016). Given the tight link between aquatic DO, land use, and diverse ecosystem functions, DO is an important indicator and holistic measure of estuarine condition.

In addition to impacts on hydrology, extreme precipitation events can increase nutrient delivery and salinity change, both of which influence DO dynamics, particularly in waterways adjacent to urban-type LULC (Walsh et al., 2005). Elevated nitrogen (N) influx can overstimulate primary production leading to short-term DO increases followed by spikes in microbial respiration and DO demand. Freshwater influx (i.e., runoff) and storm surges can alter salinity, strengthen stratification, reduce DO replenishment to bottom waters, and induce hypoxia (Rabalais et al., 2010; Wetz and Yoskowitz, 2013; Zhang et al., 2010). However, the relative importance of these processes is expected to vary among estuaries and across spatial scales. For example, in well-mixed estuaries, precipitation-driven mixing may temporarily increase DO through reaeration, whereas in stratified or poorly flushed estuaries, enhanced nutrient loading combined with reduced vertical mixing may promote sustained DO depletion.

Identifying common responses to major precipitation events across disparate estuaries can help project long-term estuarine function under increasing urbanization and more extreme precipitation events. However, such commonalities have been difficult to discern and decipher, in part, due to the complex and interconnected dynamics within individual estuaries and their varied responses to precipitation. While Ombadi and Varadharajan (2022) report contrasting effects of urbanization on salinity during flood events when regional dynamics are considered, a continental-scale study by Kaushal et al. (2018) suggests that anthropogenic activity is associated with increasing salinity in waterways. However, Kaushal et al. (2018) recognize that regional, weather-related, LULC, and geologic variabilities also influence salinization patterns. Similarly, continental-scale evaluations showed that small watersheds appear consistently less flashy than large watersheds and that there is a substantial amount of variability in these relationships at regional scale (Baker et al., 2004; Gannon et al., 2022; Hopkins et al., 2015; Poff et al., 2006). Such variation in relationships across scales may be particularly prevalent in ecosystems influenced by anthropogenic activities (Hopkins et al., 2015; Poff et al., 2006), which demonstrates the importance of considering multiple spatial scales in understanding estuarine responses to changes in precipitation patterns and watershed land use.

We aim to uncover generalizable responses to large precipitation events across five estuaries that span a gradient of urbanization and physicochemical characteristics. Using DO as an integrative indicator of ecosystem function (Abdul-Aziz and Gebreslase, 2023), we evaluate estuarine resistance to precipitation – defined here as the ability of estuaries to maintain DO stability (Isbell et al., 2015; Lake, 2013; McCluney et al., 2014; Pimm, 1984; Utz et al., 2016; Van Meerbeek et al., 2021). In the context of physicochemical and land-use factors, we evaluate estuarine resistance at: (1) the continental scale (i.e., across all estuaries); (2) across estuaries grouped by salinity; and (3) within each estuary. We hypothesize that urbanization decreases estuarine resistance to precipitation, and that relationships between resistance and physicochemical and land-use factors will vary with spatial scale and ambient salinity. This study helps advance understanding of how ongoing changes in the Earth system influence estuarine ecosystem function.

2 Methods

2.1 Estuaries and monitoring locations

We used long-term water quality monitoring data from five estuaries in the National Estuarine Research Reserve System (NOAA NERRS, 2019) to understand factors associated with their resistance to precipitation events. Lake Superior (LKS), NERR WI; Chesapeake Bay Maryland (CBM), NERR MD (Jug Bay only); Guana Tolomato Matanzas (GTM), NERR FL; Weeks Bay (WKB), NERR AL; and San Francisco Bay (SFB), NERR CA span ecoregions, land uses, and salinity (0.1–35 ppt) (Fig. 1, Table 1). Across all estuaries, there were a total of 19 monitoring locations.

Figure 1 Selected National Estuarine Research Reserve (NERR) stations. (a) Map of monitoring locations and land use/land cover within associated watersheds at Lake Superior (LKS) NERR, Chesapeake Bay, Maryland (CBM – Jug Bay) NERR, Guana Tolomato Matanzas (GTM) NERR, Weeks Bay (WKB) NERR, and San Francisco Bay (SFB) NERR. (b) Salinity from 2012 to 2022 for each monitoring location at each NERR (all n>150 000). Boxes indicate interquartile range. Means are shown in white circles. Black lines inside the boxes indicate medians.

Table 1 Land use/land cover (LULC) and population density in five National Estuarine Research Reserve System (NERRS) estuaries within 10 km of water quality and nutrient monitoring locations.

LULC definitions per ESRI (Karra et al., 2021): Water – areas where water is present throughout the year, excluding human-made structures like docks. Tree cover – vegetation with a closed or dense canopy ≥ 15 m in height. Flooded vegetation – areas where water and vegetation intermix, which are flooded seasonally or predominantly throughout the year. Crops – human-planted vegetation (e.g., cereals, grasses, and crops) that does not reach tree height. Built area – human-made structures, including roads, railroad networks, parking spaces, and industrial and residential buildings. Bare ground – areas dominated by rock, soil, sand (e.g., deserts) with sparse-to-no vegetation throughout the year. Rangeland – homogeneous grasses, mixes of vegetation below tree-height with rock and soil, or forest clearings.

Briefly, LKS NERR is a freshwater estuary situated at the confluence of the St. Louis River and Lake Superior, and includes four continuous water quality monitoring locations influenced by lake seiche: Blatnik Bridge, Oliver Bridge, Pokegama Bay, and Barker's Island. CBM NERR at Jug Bay is a semi-diurnal minimally tidal estuary. Jug Bay is located in the upper tidal reaches of the Patuxent River. It is dominated by tidal freshwater marsh, and it comprises Railroad Bridge, Iron Pot Landing, and Mataponi Creek monitoring locations. WKB NERR is a partially-mixed (i.e., lagoonal/ semi-isolated), diurnal micro-tidal estuary located on the eastern shore of the Mobile Bay. WKB NERR comprises Fish River, Weeks Bay, Middle Bay, and Magnolia River monitoring locations. SFB NERR estuary is a mixed, semi-diurnal tidal system comprising four monitoring locations across two embayments. The locations include China Camp and Gallinas Creek in the San Pablo Bay embayment within the China Camp State Park, and First Mallard and Second Mallard in the Suisun Bay embayment within the Rush Ranch Open Space Preserve. GTM NERR is a mixed, semi-diurnal tidal estuary associated with the Tolomato and Guana River estuarine systems to the north and the Matanzas River estuary to the south. Monitoring locations include Pine Island, Fort Matanzas, San Sebastian, and Pellicer Creek, with the latter site noted to experience stratification (i.e., salt-wedge).

Upstream reservoirs are reported for several monitoring locations across the set of estuaries (Table S1 in the Supplement).

2.2 Land use and land cover

To assess the relationship between resistance and urbanization, we used LULC data at 10 m resolution from ESRI (2017–2020) (Karra et al., 2021) and population density data at 100 m resolution from World Population Hub (Bondarenko et al., 2020). The NERR watersheds shapefiles were obtained from NERRS (NOAA NERRS, 2019). We analyzed LULC and population density using QGIS 3.30.3 (QGIS Geographic Information System, 2024) equipped with a semi-automatic classification plug-in, which generates percentages of each class within a given area based on the classification raster. The resulting breakdown of LULC and population density within watersheds is shown in Table S2.

Further, because the resistance index was calculated for precipitation events across short time-scales (days), large watershed areas could overrepresent hydrologic contributions to water quality at monitoring locations (watershed area range: 114–11 704 km², Table S2). To better match the temporal scale of precipitation responses, we restricted watershed characterization to a 10 km proximity zone surrounding each monitoring location. This approach captures watershed areas most likely to influence water quality over short time scales, and improves comparability among estuaries (Table 1). LULC and population density were quantified within these proximity zones. For the SFB NERR watershed, separate 10 km proximity zones were delineated for San Pablo and Suisun embayments to account for different hydrologic dynamics.

All data were collected using NERRS standard operating procedures. Briefly, water column DO (mg L⁻¹, calculated from % air saturation, temperature, and salinity at the time of measurement; thereby accounting for the impacts of temperature and salinity fluctuations on DO concentration), temperature, turbidity, salinity, and depth were measured at 15 min intervals using synchronized YSI EXO2 multiparameter sondes. Meteorological conditions, including precipitation, were also measured at 15 min intervals using NERRS standard weather station instrumentation. PO${}_{\mathrm{4}}^{\mathrm{3}-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{2}}^{-}$, and chlorophyll a (Chl a) were measured monthly from grab samples and analyzed in the lab (NOAA NERRS, 2019). Samples for nutrients were analyzed following the U.S. Environmental Protection Agency (EPA) methods (O'Dell, 1996b, a; U.S. EPA., 1993a). For Chl a analysis, the samples were collected as whole water, then filtered onto 0.45 µm pore size glass-fiber filters (0.7 µm at GTM NERR) and processed following APHA (2001) and U.S. EPA. (1993b) methods. While Chl a is also measured as chlorophyll fluorescence using an optical sensor with YSI sondes, these data were not included in our analysis due to the confounding effects of temperature, turbidity, fluorescent dissolved organic matter and other factors on sensor-based Chl a measurements (and sensor and lab-based measurements often do not correspond to each other; Dix et al., 2026). We omitted all data flagged as “suspect” or “out of range”.

Additionally, we calculated water column depth as the sum of measured water depth plus the distance between the sonde and sediment bed (Table S1). We also calculated the sum of NO${}_{\mathrm{3}}^{-}$, NO${}_{\mathrm{2}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ to assess dissolved inorganic nitrogen (DIN) concentrations and the ratio of DIN to PO${}_{\mathrm{4}}^{\mathrm{3}-}$ (hereafter, N : P).

2.3 Determination of major precipitation events

Because hydrologic dynamics and related estuarine functions can vary dramatically with annual weather conditions, we first selected one “wet” and one “dry” year for each estuary using precipitation records from nearby airports. The purpose of selecting years with disparate rainfall patterns was to encompass the maximum range of variability in expected estuarine resistance. Following Murrell et al. (2018), we calculated the long-term interquartile range (IQR, 1990–2020) of total monthly precipitation for each estuary and then selected relatively wet/dry years based on the number of months plotting above/below IQR and total annual precipitation. Possible wet and dry years were further filtered based on the completeness of NERRS data available for each estuary (Fig. S1). Long-term precipitation records included: Duluth International, Washington Reagan International, Jacksonville International, Birmingham, and San Francisco International airports available from the National Centers for Environmental Information.

Further, we selected major precipitation events within each wet and dry year by plotting daily precipitation using data from NERR meteorological stations (Fig. S2). Specifically, because the definition of “major” precipitation events is hard to quantify and it varies across estuaries, we first considered hurricanes, tropical storms, Nor'easters, atmospheric rivers, and other major storm events noted within NERR metadata sheets when selecting precipitation events. For example, at GTM NERR we focused on tropical storms Colin, Julia, and Hermine, hurricanes Matthew and Irma, and two Nor'easters (Table S3). Data availability was a second consideration – we removed possible events for which there was a substantial amount of missing data. Lastly, because metabolic and hydrologic processes vary across seasons, we chose to focus on warm season events, with the exception of SFB NERR where most precipitation occurs in the cool season but seasonal temperature fluctuations are generally lower than in other systems (Figs. S2 and S3, Table S3).

2.4 Calculation of resistance

To investigate physicochemical and urban land-use characteristics associated with estuarine responses to precipitation events, we calculated the resistance index described in Orwin and Wardle (2004). The resistance index is a normalized parameter (−1 to +1) describing the magnitude of shift in a response variable from an initial condition. It has been used across a wide variety of ecosystems and response variables, including aquatic ecosystems (Thayne et al., 2022, 2023; Tsai et al., 2011). The resistance index is calculated as:

$$\begin{array}{}\text{(1)}& \mathrm{Resistance}=\mathrm{1}-{\displaystyle \frac{\mathrm{2}\left|D_{\mathrm{0}}\right|}{\left(C_{\mathrm{0}}+\left|D_{\mathrm{0}}\right|\right)}}\end{array}$$

where C₀= concentration of the response variable pre-disturbance and D₀= difference between the concentration of the response variable pre- and post-disturbance (P₀) (i.e., $D_{\mathrm{0}}=C_{\mathrm{0}}-$ P₀). An index value of +1 indicates the highest possible resistance. Index values between 0 and 1 show that the magnitude of response variable shift is less than the magnitude of the baseline (i.e., $\left|D_{\mathrm{0}}\right|\le C_{\mathrm{0}}$). A resistance index of 0 indicates that the shift in the response variable is equivalent to the magnitude of the baseline (i.e., $\left|D_{\mathrm{0}}\right|=C_{\mathrm{0}}$), whereas index values between 0 and −1 reflect that change in the response variable is greater than the magnitude of the baseline (i.e., $\left|D_{\mathrm{0}}\right|>C_{\mathrm{0}}$). Overall, index values closer to 1 indicate more resistant systems (Orwin and Wardle, 2004).

While the resistance index is indicative of the ability of a system to maintain its pre-disturbance functions and ecological state (i.e., ecosystem stability), we emphasize that it is a normalized value that does not in itself convey information about overall estuarine function. Similarly, the resistance index alone cannot infer if a shift has a positive or a negative effect on ecosystem function. The resistance index enables the comparison of the amount of change induced by disturbance across vastly different estuaries.

In parallel, it is also useful to consider the absolute value of the response variable (in this case DO) pre- and post-disturbance, which conveys information on the ambient state of an estuary and the directionality of its response to the disturbance. We therefore present C₀ and P₀ (Figs. 2a–e and S4, Table S4) to define differences in DO within and across estuaries. We also present the resistance index, which pairs C₀ and P₀ values for the same event (Fig. 2f–j), to understand the ecosystem stability within and across estuaries after precipitation.

Figure 2 Variation in pre- and post-disturbance distribution of dissolved oxygen and resistance within individual estuaries. Boxes show the quartiles of the dataset and the whiskers show the rest of the distribution. Means are shown in white circles, and medians are shown in black solid lines. (a–e) Distribution of dissolved oxygen concentrations prior to (C₀, yellow boxes) and post- (P₀, pink boxes) precipitation. (f–j) Resistance index across monitoring locations at: Lake Superior (LKS) NERR (all n=7), Chesapeake Bay, Maryland (CBM – Jug Bay) NERR (all n=10), Guana Tolomato-Matanzas (GTM) NERR (all n=7), Weeks Bay (WKB) NERR (all n=15), San Francisco Bay (SFB) NERR (n=8 except Second Mallard n=7).

Because DO is critical to myriad functions that regulate aquatic ecosystems (Abdul-Aziz et al., 2007; Abdul-Aziz and Gebreslase, 2023; Caffrey, 2004; Mulholland et al., 2001; Murrell et al., 2018; Odum, 1956) and because bi-directional flow makes estuarine metabolism virtually impossible to model (Loken et al., 2021), we calculated the resistance index using temperature- and salinity-adjusted DO (mg L⁻¹). Additionally, the resistance index is highly sensitive to the researcher-defined baseline and post-disturbance time periods that are used to calculate C₀ and D₀. Therefore, we calculated C₀ as the average DO during a manually-curated timespan preceding each precipitation event (∼24 h to 6 d, without precipitation) where DO dynamics looked unremarkable (Fig. S5, Table S3). The timespan for estimating P₀ was also manually selected in the context of each event, defined here as the maximum displacement from C₀ during and after the event (Fig. S5, Table S3). We also verified that the resistance index calculated using DO (mg L⁻¹) versus DO (% sat.) showed no notable differences (Fig. S6).

2.5 Statistical analysis

To compare resistance, nutrient concentrations, and concentrations of DO pre- and post-precipitation within and among estuaries, we used ANOVA with post-hoc Tukey HSD or Kruskal-Wallis test, as appropriate based on the Shapiro-Wilk normality test.

To test specific physicochemical and land-use factor associations with estuarine resistance, we used linear regressions at continental and local scales independently (i.e., all estuaries combined versus within each individual estuary) and within salinity-based groups (i.e., using a threshold of average annual salinity <10 or ≥10 ppt calculated for wet and dry years separately). Specifically, we tested the relationships between resistance and turbidity, salinity, water temperature, water column depth, DIN, PO${}_{\mathrm{4}}^{\mathrm{3}-}$, N : P, Chl a, percent trees, percent crops, percent built area, and population density. For linear regressions that involved land-use factors, we used LULC and population density from the 10 km proximity zone (Table 1). Because of discrepancies in the temporal intervals at which measurements were collected, we conducted separate regressions involving turbidity, salinity, water temperature, water column depth (collected at 15 min intervals) versus DIN, PO${}_{\mathrm{4}}^{\mathrm{3}-}$, N : P, and Chl a (collected at monthly intervals) predictors. The low-salinity group included all LKS NERR and CBM NERR locations, First Mallard and Second Mallard locations at SFB NERR during wet and dry years, as well as Fish River, Middle Bay, Magnolia River, and Weeks Bay (Weeks Bay during the dry year only) at WKB NERR estuary. The high-salinity group included all monitoring locations at GTM NERR, China Camp and Gallinas Creek at SFB NERR estuary for wet and dry years, and Weeks Bay at WKB NERR during the wet year.

We used annual mean values for continental-scale, local-scale, and salinity-based groups regressions involving nutrients and Chl a because monthly sampling intervals of lab-based measurements did not always correspond with selected precipitation events. This resulted in two data points per monitoring location for nutrients and Chl a regressions. While not directly associated with any particular precipitation event, relationships of nutrients and Chl a with resistance values analyzed on an annual basis carry valuable information about how the ambient conditions of the system can impact an estuary's response to precipitation. This knowledge is essential for deriving and testing hypotheses that describe why a certain estuary may respond to a storm event in a particular way.

For regressions involving sensor-based measurements, we attempted to provide as much resolution as possible into specific events. We therefore used values averaged over the event-specific baseline periods (Table S3) for sensor measurements across scales and salinity-based groupings. Values averaged over event-specific baseline periods were matched with corresponding event-specific resistance index resulting in one data point per precipitation event for each monitoring location. We did not include LULC and population density as predictors of resistance at individual estuaries because of the close proximity of some monitoring locations to one another (<800 m).

To ensure the robustness of linear regression analysis, we assessed the dataset for multicollinearity (Fig. S7) and conducted principal component analysis (PCA) with Varimax rotation on loadings (Fig. S8). For PCA, we used sensor-based, 15 min interval measurements (temperature, salinity, turbidity, water column depth) and omitted the nutrient data collected at monthly intervals to keep most data points in our analysis while remaining consistent with continental-scale linear regressions. Most variables were largely uncorrelated and the PCA results agreed with continental-scale linear regressions.

Statistical analyses and visualizations were performed in Python 3.10.11 using scipy, statsmodels, seaborn, scikit-learn, and matplotlib libraries.

3 Results

3.1 Land use/land cover and nutrient concentrations across estuaries

The five study estuaries spanned strong gradients in LULC and nutrient concentrations. Both embayments of the SFB NERR had higher percentages of urban-type land (e.g., built area) and population density than any other estuary, followed by CBM NERR (Table 1). Agricultural land was more prevalent at WKB NERR (30.36 %) compared to other estuaries. GTM NERR and LKS NERR had more mixed LULC, with high proportions of tree cover. SFB and CBM NERRs had high DIN concentrations compared to LKS and GTM NERRs (p<0.01). Mean DIN values at all SFB and CBM NERRs monitoring locations were >0.50 mg-N L⁻¹ versus <0.16 mg-N L⁻¹ at LKS and GTM NERRs (Fig. S9). Phosphate concentrations were the highest at SFB NERR (mean across all monitoring locations =0.135 mg-P L⁻¹, SD =0.08; means at all other estuaries <0.03 mg-P L⁻¹, p<0.001, Fig. S9). Overall, mean N : P across LKS, CBM, GTM, WKB, and SFB NERRs estuaries was 28.04 (SD =28.08), 26.56 (SD =10.75), 2.84 (SD =0.91), 83.82 (SD =55.27), and 5.10 (SD =1.4), respectively (Fig. S10), with GTM and SFB NERRs indicating N-limiting conditions based on a mass-based 7.2N : 1P Redfield ratio (Redfield, 1934).

3.2 Changes in dissolved oxygen and resistance to precipitation across estuaries and monitoring locations

Dissolved oxygen concentrations pre- and post-precipitation differed across all estuaries, when evaluating the overall directionality (pre- [C₀]: F=45.6, p<0.0001; post- [P₀] : F=17.1, p<0.0001; Fig. 2a–e). Across all precipitation events (i.e., non-specific to an event), SFB NERR had the highest pre-precipitation DO concentration of all estuaries. Generally, at SFB and CBM NERRs (more urban) DO declined following precipitation (p<0.01). At LKS NERR, precipitation events significantly increased DO concentration (F=5.1, p=0.03). There was no significant difference between the overall pre- and post-precipitation DO concentrations at GTM and at WKB NERRs (more agricultural) (GTM NERR: F=0.5, p=0.48; WKB NERR: F=1.23, p=0.27).

Resistance also differed across estuaries ($F(\mathrm{4},\mathrm{172})=\mathrm{26.3}$, p<0.0001, Fig. 2). SFB NERR monitoring locations had the highest mean resistance (mean =0.66), while monitoring locations at WKB NERR were the least resistant (mean =0.23). Within individual estuaries, resistance varied across monitoring locations at LKS NERR ($F(\mathrm{3},\mathrm{24})=\mathrm{13.9}$, p<0.0001), CBM NERR ($F(\mathrm{2},\mathrm{27})=\mathrm{7.9}$, p<0.01), and SFB NERR ($F(\mathrm{3},\mathrm{27})=\mathrm{10.2}$, p<0.001) but was not significantly different between monitoring locations within GTM NERR (H(3)=3.5, p=0.32) and WKB NERR ($F(\mathrm{3},\mathrm{56})=\mathrm{0.5}$, p=0.65) (Fig. 2f–j). Resistance was most variable across monitoring locations at LKS NERR (−0.26 to 0.89) and least variable across locations at SFB NERR (0.44 to 0.86).

3.3 Continental, salinity-based, and local relationships between resistance and physicochemical factors, land use/land cover, and population density

When data from the five estuaries were considered together (i.e., continental scale), we found generally weak but significant positive relationships between resistance and water column depth (p<0.0001, R²=0.11), turbidity (p=0.0014, R²=0.06), log(DIN) (p=0.037, R²=0.12), percent built area (p=0.014, R²=0.16), and population density (p<0.001, R²=0.29), as well as significant negative relationships with water temperature (p<0.0001, R²=0.42) and Chl a (p<0.0001, R²=0.43) (Fig. 3, Table S5).

Figure 3 Continental-scale and salinity-based relationships of resistance with physicochemical factors, land use/land cover, and population density. Continental-scale regressions considered all monitoring locations across all estuaries. Significant relationships (p<0.05) are shown in dotted green, solid red, and solid blue lines for continental scale, high-salinity and low-salinity estuaries, respectively. Standard errors of the mean are shown in vertical and horizontal lines. The LULC parameters and population density were used from within the 10 km proximity zone adjoined to the monitoring locations. Please refer to Table S5 for resulting statistics and to Fig. S7 for information on covariance among predictor variables.

When grouped by salinity, estuarine resistance was more strongly related to physicochemical and land-use factors (Fig. 3, Table S5). Within low-salinity estuaries, resistance was positively related to depth of the water column (p<0.0001, R²=0.20), which is consistent with continental-scale results; and negatively related to salinity and percent cropland (p<0.0001, R²=0.23, and p=0.031, R²=0.19, respectively), relationships not found on continental scale. Within high-salinity estuaries, mean resistance showed positive relationships with annual mean log(DIN) (p=0.004, R²=0.54) and percent built area (p<0.001, R²=0.73), consistent with continental-scale results. Also, resistance in high-salinity estuaries was positively related to turbidity (p<0.0001, R²=0.28). Observations present in high-salinity estuaries but absent from continental-scale evaluations included a negative relationship between mean resistance and tree cover (p=0.012, R²=0.45), as well as a positive relationship with N : P (p<0.0001, R²=0.80). Additionally, mean resistance in both low- and high-salinity groups was positively related to population density (p=0.033, R²=0.18, and p<0.001, R²=0.73, respectively), and negatively related to water temperature (p<0.0001, R²=0.36, and p<0.0001, R²=0.58, respectively) and Chl a (p<0.001, R²=0.49, and p=0.042, R²=0.32, respectively). The temperature and Chl a relationships were consistent with continental-scale results. Generally, the low-salinity estuaries had fewer significant relationships between resistance and physicochemical factors, LULC, and population density compared to high-salinity estuaries. The strength of the relationships between resistance, physicochemical and LULC factors generally increased with estuarine salinity. In the high-salinity group, the coefficient of determination (R²) was above 0.5 in five of eight significant relationships, whereas all significant relationships identified in the low-salinity group and in the continental-scale analysis had R² values below 0.5.

At local scales (i.e., within each estuary), resistance was related to some physicochemical factors not observed in continental-scale or salinity-based groups. The number, strength, and direction of relationships between local-scale resistance and physicochemical factors varied substantially across estuaries (Figs. 4, S11, S12). GTM NERR had the highest number of significant relationships between resistance and physicochemical factors. Resistance was negatively related to water temperature (p<0.001, R²=0.52), PO${}_{\mathrm{4}}^{\mathrm{3}-}$ (p=0.013, R²=0.67) and Chl a concentrations (p=0.009, R²=0.71), as well as positively related to salinity and N : P (p=0.02, R²=0.19 and p=0.009, R²=0.71, respectively). In contrast, at CBM NERR, resistance was related only to water column depth (positive, p=0.013, R²=0.21). At LKS NERR, resistance was positively related to water column depth (p<0.0001, R²=0.50) and negatively related to turbidity (p=0.004, R²=0.28). At WKB NERR, relationships between resistance and water temperature, salinity, turbidity, and Chl-aconcentrations were all negative (p<0.001, R²=0.41; p=0.047, R²=0.07; p=0.002, R²=0.16; and p=0.03, R²=0.58, respectively). At SFB NERR, resistance was negatively related to water column depth (p<0.0001, R²=0.48), and positively related to salinity (p=0.00016, R²=0.42). There was no overarching relationship between resistance and total precipitation amount of each event, except for significant but weak negative relationships at GTM and WKB NERRs estuaries (at GTM NERR: p=0.026, R²=0.18; at WKB NERR: p=0.002, R²=0.15; Fig. S13).

Figure 4 Relationships between resistance and physicochemical factors for each estuary. Significant correlations (p<0.05) are shown with black lines. Standard errors of the mean are shown in vertical and horizontal black lines for relationships using annual means for chlorophyll a (Chl a), PO${}_{\mathrm{4}}^{\mathrm{3}-}$, and N : P. For additional results see Figs. S11, S12, and S14.

In summary, some relationships between estuarine resistance and physicochemical factors and urban land use appeared to be more generalizable while others varied across scales (Fig. 5). For instance, resistance was related to water temperature, water column depth, turbidity, and Chl a across continental and local scales, and within salinity-based estuary groups. Dissolved inorganic N, percent built area, and population density were all related to resistance at the continental scale and in salinity-based groups. Both local and salinity-based evaluations revealed that N : P, salinity, and turbidity were related to resistance. Unique relationships included PO${}_{\mathrm{4}}^{\mathrm{3}-}$ at the local scale (GTM NERR only), tree cover in high-salinity and percent cropland in low-salinity estuaries. The strength of the relationships varied among groups, with higher coefficients of determination present in the high-salinity estuarine group (i.e., most R²>0.5, Fig. 3, Table S5).

Figure 5 Cross-scale significant relationships between estuarine resistance and physicochemical and land-use factors. (a) Venn diagram of resistance relationships with physicochemical and land-use factors in high- vs. low-salinity estuaries. Positive or negative relationships are indicated with “+” and “−”, respectively. Italicized text identifies relationships with coefficient of determination (R²)<0.5 across both high- and low-salinity groups. Factors in the overlap zone with regular text have contrasting R² across groups. (b) Venn diagram of relationships between resistance and physicochemical and land-use factors in continental, local, and salinity-based groups. Estuarine resistance with land use/land cover and population density marked with asterisks (^*) were not evaluated at local scale due to overlap in the spatial domains of some monitoring locations. Italicized text identifies continental-scale relationships with R²<0.5.

4 Discussion

Discerning estuarine response patterns following precipitation is important for predicting future impacts of urbanization and changing global weather patterns on estuarine ecosystems. Previous studies have shown that patterns identified at large scales may not be applicable across different ecoregions, geologies, and/or other ecosystem factors (Baker et al., 2004; Gannon et al., 2022; Hopkins et al., 2015; Kaushal et al., 2018; Ombadi and Varadharajan, 2022; Poff et al., 2006). Our results underscore the importance of cross-scale evaluations that can elucidate commonalities in estuarine response to precipitation, as well as variability in the factors associated with resistance across individual estuaries.

We show that while relationships between estuarine resistance, urbanization, and DIN may exist at the continental scale, they do not consistently apply within individual estuaries. This is because an estuary's resistance depends on myriad specific factors in addition to many factors identified at larger scales. In contrast to our overarching hypothesis, we found that urbanized estuaries tend to have higher resistance than pristine estuaries. These results suggest that the effects of watershed urbanization may impact estuarine stability either by providing a mechanism that allows estuaries to dampen large shifts in DO in response to major precipitation events; or by disturbing the baseline DO to an extent where even a major precipitation event would not produce a significant shift in DO, making the system appear highly resistant.

4.1 Dissolved oxygen dynamics differ between estuaries influenced by urban or agricultural land use and those with minimal urbanization

There were vast differences in geometry, circulation, and hydrologic conditions between urbanized estuaries and estuaries surrounded by agricultural land in this study. Yet, despite these differences, several estuaries exhibited similar DO responses to precipitation. Precipitation generally reduced DO concentrations at SFB and CBM NERRs (Fig. 2a–e), whereas LKS NERR, the estuary with the lowest proportion of urbanized land cover, experienced an overall increase in DO concentration following precipitation.

While a wide range of physical factors can impact DO in estuaries including channel geometry, river discharge, circulation, wind patterns, residence time, and upstream reservoirs (Kemp and Boynton, 1980; Raimonet and Cloern, 2017; Raymond et al., 2012; Raymond and Cole, 2001; Scully, 2010; Zheng et al., 2024), urban estuaries in particular are often impacted by a combination of these processes. Urban estuaries often serve as basins for wastewater treatment outflows, which can supply continued freshwater discharge and nutrients during baseflow. Both urban and agriculturally influenced estuaries are prone to increased nutrient loading during and shortly after precipitation events (Chapin et al., 2004; Costanzo et al., 2003; Mallin et al., 2009), which impacts primary production and microbial metabolism and may lead to declines in DO concentrations (e.g., algal blooms). Similarly, power-generating dams and other large reservoirs often present in urban estuaries can strongly influence DO dynamics through their regulation of flow conditions, salinity, temperature, microbial processes, and solute residence time (Abbott et al., 2022; Ferencz et al., 2021b, a).

Therefore, closely examining the directionality of the change in DO response to precipitation is important for informing the environmental significance of storm events for a specific estuary. An increase in DO following precipitation could reduce hypoxia through water column reaeration (Bianucci et al., 2018; Bohórquez-Bedoya et al., 2024) and have positive effects on nutrient cycling and growth conditions for macro- and micro-fauna (Harris et al., 2015). A decrease in DO following precipitation may reflect the opposite effects (e.g., hypoxia), leading to changes in microbial community structure, disruptions in nutrient cycles, and fish die-off (Diaz et al., 1992; Llansó, 1992). Further, the extent to which ecosystem function is linked to DO changes is intertwined with its ambient ecological state and the magnitude of the DO change. An increase in DO following a storm may be less consequential for estuaries with DO-replete ambient conditions as compared to estuaries with chronically depleted DO. This implies that a comprehensive assessment of storm impacts on estuarine ecosystem stability and function should consider resistance index alongside historical DO conditions and the magnitude and directionality of DO change.

For example, resistance index values for the more pristine and tree-cover dominated Pokegama Bay at LKS NERR were comparable to most monitoring locations at WKB NERR, which were associated with agriculture. However, Pokegama Bay and WKB NERR experienced different patterns in DO change following precipitation (Fig. 2). While the full nature of the causes and impacts of DO change in each estuary is beyond the scope of this study, the differences in absolute DO change suggest that different underlying processes influence stability in estuarine dynamics in response to precipitation across different estuaries. Our goal is to identify common predictors of estuarine resistance at the continental scale despite substantial variability in processes, while also identifying salinity-dependent and estuary-specific predictors that help elucidate dynamics that may be particularly influential in individual systems.

4.2 Urbanization and dissolved inorganic nitrogen correspond with elevated resistance to precipitation at the continental scale

Higher resistance to precipitation in the most urbanized estuaries and overarching relationships between urban LULC and resistance across all estuaries reflect that estuaries within urban watersheds typically experience relatively small shifts in DO following precipitation (Figs. 2, 3). These results contradict our hypothesis that urban estuaries should show low resistance to precipitation because of greater physical and chemical disturbances like flashiness, streambed scouring, N loading, and turbidity that impact DO dynamics (Bernhardt et al., 2008; Groffman et al., 2004; Hession et al., 2003; Hopkinson and Vallino, 1995; Walsh et al., 2005).

It is possible that watershed urbanization could equip estuaries with adaptations that help dampen precipitation impacts on shifts in DO. For instance, an increase in flashiness in waterways would increase flow velocity and water column reaeration (Raymond et al., 2012; Raymond and Cole, 2001) and contribute to phytoplankton removal via transport or turbidity-driven light attenuation (Caffrey, 2004; Pennock and Sharp, 1986). This impact may be particularly important for estuaries whose baseline DO conditions are influenced by algal blooms, which induce large diel DO fluctuations by overproducing DO during the day and severely depleting it at night (Chapin et al., 2004; Ni et al., 2020). Managing phytoplankton overgrowth could help maintain DO closer to baseline conditions. Supporting this explanation, turbidity-driven limitations on phytoplankton were previously reported for the SFB estuary (Cloern, 1987). We also found that turbidity was positively related to built area, Chl a was negatively related to population density, and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ concentration was negatively related to tree-cover (Figs. S7, S15). Tree-populated riparian zones, often reduced or absent in urban environments, can help regulate phytoplankton overgrowth and DO dynamics by significantly reducing PO${}_{\mathrm{4}}^{\mathrm{3}-}$ concentrations (Keller and Fox, 2019).

Relatively high levels of ambient DIN concentrations and altered baseline DO dynamics in urban estuaries may also contribute to smaller DO responses to precipitation. Here, estuarine DIN concentrations were positively associated with urbanization and were significantly, but weakly, related to resistance at the continental scale (Fig. 3, S15). While moderate levels of N support biological metabolisms (Howarth, 1988; Howarth and Marino, 2006; Vitousek and Howarth, 1991; Zhang et al., 2021), nutrient loading, including N, is a significant problem for urban aquatic environments (Beman et al., 2005; Bernot et al., 2010; Bettez et al., 2015; Black et al., 2011; Hopkinson and Vallino, 1995; Mulholland et al., 2008; Reisinger et al., 2016). For urban estuaries, where N loading is associated with chronically low DO availability, this could mean that precipitation-driven increases in DO under simultaneous high N input would be small, which would result in higher resistance. Nitrogen limitation, in contrast, is prevalent in coastal marine systems and can negatively impact overall estuarine function by constraining biological activity (Elser et al., 2007; Guildford and Hecky, 2000; Paerl, 2018; Paerl and Piehler, 2008). Relationships of resistance with N and urbanization were particularly evident in more coastal (high-salinity) estuaries where we found positive relationships between resistance, N : P, and built area (Fig. 3). Thus, altered N dynamics in urban estuaries may be key to understanding their comparatively small changes in response to precipitation.

While the mechanisms behind the positive relationship between resistance and N are not clear, changes in concentrations of NH${}_{\mathrm{4}}^{+}$ were previously suggested to regulate algal blooms in the SFB estuary (Dugdale et al., 2007; Parker et al., 2012). Elevated NH${}_{\mathrm{4}}^{+}$ concentrations resulting from changes in water-treatment practices in San Francisco were shown to inhibit NO${}_{\mathrm{3}}^{-}$ uptake and suppress algal blooms, thereby altering a major control on DO dynamics in SFB. Therefore, the prevalent N species and shifts in DIN composition following precipitation events should be considered in attempts to explain the underlying mechanisms behind the positive relationship between DIN and resistance.

Alternatively, we also consider whether baseline DO concentrations in urban estuaries are already disturbed to such an extent that precipitation events cause minimal further disruption. Urbanization itself can promote large diel DO fluctuations, as suggested by Gold et al. (2020), and if a precipitation event induces fluctuations that are similar in magnitude, the post-precipitation DO concentrations will not deviate significantly from baseline conditions. Under such a scenario, urban estuaries will appear more resistant to precipitation events compared to more pristine systems.

4.3 Water column depth, temperature, turbidity, and Chl a are related to resistance across all scales

We found generalizable patterns (across all scales), in which resistance was positively but weakly related to water column depth and turbidity and negatively related to water temperature and Chl a (Figs. 3–5). As discussed above, phytoplankton dynamics may be an important factor in regulating DO in response to disturbance due to their tight linkage with N loading, which is consistent with results previously reported by others (Thayne et al., 2023). This is further underscored by the existence of a relationship between Chl a and resistance across continental, salinity-based, and local analyses. In parallel, the relationship between resistance and turbidity was positive at the continental scale and in the high-salinity group (Fig. 3), while at LKS resistance and turbidity correlated negatively (Figs. 3, 4). The positive relationship between resistance and turbidity and simultaneous negative relationship with Chl a can be attributed to turbidity-driven light attenuating conditions that restrict phytoplankton growth, which is further supported by a negative relationship between turbidity and Chl a (Fig. S7).

Water column depth can also influence estuarine resistance to precipitation through its effects on DO dynamics. Through dilution, deeper estuaries may show a more buffered response to short-term hydrologic changes following precipitation by attenuating changes associated with freshwater inflow, nutrient loading, turbulence, and gas exchange. Similarly, deeper estuaries generally have longer equilibration time with environmental conditions and less diel variability in parameters like temperature (Caissie, 2006; Macan, 1958). These characteristics can contribute to more stable baseline DO conditions and more moderate responses to precipitation. However, greater depth can also enhance vertical stratification and reduce oxygen replenishment to bottom waters, potentially increasing susceptibility to hypoxia (Cloern, 2001; Cloern and Jassby, 2010). Deeper estuarine channels are frequently associated with urbanization, which exhibited high resistance in this study, because channel incision, dredging, and altered sediment transport are more common in urban and agricultural watersheds (O'Driscoll et al., 2010; Simon and Rinaldi, 2006; Walsh et al., 2005). Accordingly, water column depth appears to be an important factor influencing estuarine resistance to precipitation, although its effects on DO dynamics depend on interactions with stratification, hydrodynamics, and watershed characteristics.

Lastly, because temperature controls various chemical and biological processes that impact DO availability (e.g., microbial growth, oxygen solubility), many studies have focused on the impact of rising temperature on DO dynamics in estuaries (Apple et al., 2006; Caffrey, 2003; Caffrey et al., 2014). Previous studies have linked changes in weather patterns to thermal pollution of aquatic systems following rain events (Zahn et al., 2021) and link elevated global temperatures to decreased primary production (Song et al., 2018). As such, estuaries with elevated ambient temperatures may have a decreased capacity to resist disturbances to DO dynamics relative to estuaries with more moderate temperatures.

4.4 High variability in factors associated with resistance at local scales

The substantial differences in physicochemical factors linked to resistance across individual estuaries highlight the need to consider both generalized and local estuarine patterns when predicting and managing the effects of major precipitation events and watershed urbanization on water quality (Figs. 3–5). Although several factors were significantly associated with resistance, these relationships explained relatively little of the observed variance and may have been influenced by non-independence among reported observations from the same monitoring sites and estuaries. Individual estuaries may require consideration of factors beyond water temperature, water column depth, turbidity, and Chl a – identified here as generalizable at the continental scale – and be mindful of differences in the directionality of relationships between physicochemical factors and resistance among estuaries. For example, we found that resistance at CBM NERR is related to one factor (water column depth), while at GTM NERR, resistance is related to five factors (salinity, water temperature, Chl a, PO${}_{\mathrm{4}}^{\mathrm{3}-}$, and N : P), of which three were not identified at the continental scale. Similarly, although resistance and salinity are positively related at GTM and SFB NERRs, at WKB NERR this relationship is negative. Likewise, while the relationship between resistance and water column depth is negative at SFB NERR, at CBM and LKS NERRs this relationship is positive.

Variability in responses to stressors, at both the system and individual factor levels, is common. Understanding how individual physicochemical factors respond to precipitation broadly and within each system provides insight into the observed relationships between those factors and resistance. For example, Chl a and dissolved inorganic phosphorus have been shown to respond more strongly to storms in some estuaries than in others due to differences in light availability, grazing, and nutrient concentrations (Chen et al., 2015; Cloern, 2001; Cloern and Jassby, 2010; Dix et al., 2013; Liao et al., 2021; Zhang et al., 2022). Likewise, water temperature and salinity also have variable responses to precipitation (Buelo et al., 2023; Chen et al., 2015; Dix et al., 2008), leading to differences in biological processes and phytoplankton activity across estuaries (Apple et al., 2008). System dependencies on groundwater discharge – as a driver of nutrient inputs and a regulator of DO dynamics (Brookfield et al., 2021; Kornelsen and Coulibaly, 2014) – or on microbial community structure (Cheng et al., 2021) should also be appraised when evaluating factors associated with urban estuaries' resistance to precipitation.

4.5 Conceptual model for resistance to precipitation based on estuary characteristics

Through observed associations between resistance and physicochemical and land-use factors in this dataset, we propose the following conceptual model for expected resistance based on estuarine characteristics (Fig. 6, Table S1). This conceptual model is intended as a hypothesis-generating framework for future work rather than a deterministic classification of estuarine behavior. We suggest that the highest resistance may occur in tidal/marine-dominated or urban-influenced estuaries that are well-mixed and have relatively short residence time, similar to SFB NERR. In such estuaries, the disturbances to DO dynamics driven by precipitation events, riverine freshwater inputs and/or terrestrial-aquatic exchange across mixed land uses (e.g., changes in salinity, turbidity, and/or nutrients), are mediated by tidal mixing. In contrast, the lowest resistance index could be expected in enclosed, shallow, poorly mixed (i.e., stratified) estuaries with long residence time, strong influence from river systems, and/or surrounded by agricultural LULC, similar to WKB NERR, Pellicer Creek (GTM NERR), or Pokegama Bay (LKS NERR). In such estuaries, the abrupt disruptions to ambient salinity, turbidity, and/or nutrient concentrations driven by precipitation events are expected to be associated with longer-lasting changes and larger fluctuations in DO, which in turn would result in low resistance.

Figure 6 Conceptual model for expected estuarine resistance to precipitation events based on estuarine characteristics.

Because estuaries are characterized by a continuum of interacting variables, we expect that many estuaries will show moderate resistance that varies with precipitation events of different sizes. Based on the associations revealed in this dataset, medium-high resistance may occur in well-mixed estuarine systems with engineered flow, relatively short residence time, and periodic nutrient loading conditions (e.g., CBM NERR, remaining LKS and GTM NERRs sites). Conversely, medium-low resistance may be more common in systems with partial stratification, periodic nutrient loading and moderate residence time. However, these patterns are context-dependent and likely influenced by additional interacting factors not explicitly resolved here.

Although numerous estuary-specific interactions can produce different responses to stressors, our results provide a framework for understanding generalizable responses to precipitation based on measurable characteristics, while also highlighting the importance of system variability in determining estuarine resistance. We recognize that additional scales of investigation (e.g., regional) and other ecological, hydrological and environmental factors could further refine predictions of estuarine resistance. Incorporating these elements in future studies could strengthen our ability to anticipate estuarine responses to disturbances.

5 Conclusions

In light of increasing urbanization and changing rainfall patterns, cross-scale evaluations of estuarine responses to precipitation events are imperative for developing effective management strategies. We observed that urban estuaries tended to show higher resistance to changes in DO in response to precipitation than more pristine estuaries; although the relationships varied across systems. Cross-scale analyses identified that water column depth, water temperature, turbidity, and Chl a variables broadly associate with estuarine resistance. However, across different scales, we found that system variability results in additional factors that are important to consider when managing the responses to major precipitation events of individual estuaries. Based on our results, we suggest that investigations targeting N, nutrient delivery mechanisms, and microbial activity could help understand what makes urban estuaries more resistant to large shifts in DO following precipitation. We also highlight that high-resolution water quality and nutrient data surrounding precipitation events, along with careful consideration of local variability and models for system responses to precipitation, are needed to help elucidate the underlying mechanisms for high resistance of urban estuaries. This study provides a basis for improving guidelines and predictive capabilities addressing system response to future weather and urbanization scenarios.