Linking watershed nutrient loading to estuary water quality with generalized additive models

Location and data

Lavaca Bay is a 190 km² estuary with the majority of freshwater inflow provided by the Lavaca and Navidad River systems (Fig. 1). Lavaca Bay is relatively shallow with an average depth of 1.2 m and a generally well mixed and turbid water column (Beseres Pollack, Palmer & Montagna, 2011; Montagna et al., 2020). The Garcitas-Arenosa, Placedo Creek, and Cox Bay subwatersheds provide additional freshwater inflows. The proximity to freshwater inputs in the upper Lavaca Bay results in a strong mean salinity and variance gradient (Montagna et al., 2020). Average salinity ranges from 12 psu near the Lavaca River mouth to 23 psu at the mouth of Lavaca Bay. The variance is inversely correlated with mean salinity. The entire watershed area is 8,149 km² and primarily rural. Watershed land cover and land use is 50% grazed pasture and rangeland, 20% cultivated cropland (primarily rows crops such as corn, cotton, and sorghum), and 5% suburban/urban. Pasture and rangeland is concentrated in the Lavaca River watershed, while cultivated crops are generally located along the eastern tributaries of the Navidad river. The Lavaca and Navidad River watersheds are a combined 5,966 km², or approximately 73% of the entire Lavaca Bay watershed area. Discharge from the Navidad River is regulated by Lake Texana which has been in operation since 1980. Lake Texana provides 0.210 km³ of water storage and discharges into the tidal section of the Navidad River which ultimately joins the tidal section of the Lavaca River 15 km upstream of the confluence with the Lavaca Bay.

Daily discharges for the Lavaca River (USGS-08164000, Fig. 1) were obtained from the United States Geologic Survey (USGS) National Water Information System (https://waterdata.usgs.gov/nwis) using the dataRetrieval R package (De Cicco et al., 2022). Gaged daily discharges from the outlet of Lake Texana on the Navidad River (USGS-0816425) were provided by the Texas Water Development Board (TWDB) (April 21, 2022 email from R. Neupane, TWDB).

Water quality data for the two freshwater and three estuary locations were obtained from the Texas Commission on Environmental Quality (TCEQ) Surface Water Quality Monitoring Information System. Data submitted through the system are required to be collected under Quality Assurance Project Plans and lab method procedures outlined by the TCEQ’s procedures manual. These operating procedures ensure consistent collection and laboratory methods are applied between samples collected by different entities and under different projects. All sites had varying time periods and availability of water quality data. For freshwater locations, TP from January 2000 through December 2020 and nitrate-nitrogen (NO₃) data from January 2005 through December 2020 were downloaded (Table 1). Less than 5-years of total nitrogen and TKN concentration data were available at the freshwater sites and deemed insufficient to develop load estimation models (Horowitz, 2003; Snelder, McDowell & Fraser, 2017). The three estuary sites included an upper-bay site near the outlet of the Lavaca River system (TCEQ-13563), a mid-bay site (TCEQ-13383), and the lower-bay site near the mouth of the Bay (TCEQ-13384). For estuary locations, we obtained data for TP, Nitrite+Nitrate(NO_X), TKN, chlorophyll-a, and DO concentrations from January 2005 through December 2020 (Table 2). Freshwater nutrient samples were collected 0.3 m below the surface or halfway between the surface and bottom if depth was insufficient. Estuary nutrient samples were all collected 0.3 m below the surface. Estuary DO measurements were made at 3 locations along the water column (0.3 m, mid-depth, and bottom) and averaged. All water quality data was collected approximately quarterly.

Estimating watershed based nutrient loads

Estimates of NO₃ and TP loads at the Lavaca River (USGS-08164000) and the outlet of Lake Texana on the Navidad River (USGS-08164525) were developed using GAMs relating nutrient concentration to river discharge, season, and time. Separate models were fit at each station for each parameter and used to predict nutrient concentrations for each day in the study period. GAMs were fit using the mgcv package in R which makes available multiple types of smooth functions with automatic smoothness selection (Wood, 2011). The general form of the model related NO₃ or TP concentration to a long term trend, season, streamflow, and two different antecedent discharge terms, shown in Eq. (1):

${\displaystyle g\left(\mu \right)=\alpha +f_1\left(ddate\right)+f_2\left(yday\right)+f_3\left(log\left(1+Q\right)\right)+f_4\left(ma\right)+f_5\left(fa\right),}$ (1)${\displaystyle y\sim \mathcal{N}\left(\mu ,{\sigma }^2\right)}$

where μ is the conditional expected NO₃ or TP concentration, g() is the log-link, α is the intercept, f_n() are smoothing functions. y is the response variable (NO₃ or TP concentration) modeled as normally distributed with mean μ and standard deviation σ. ddate is the date converted to decimal notation, yday is numeric day of year (1–366), and Q is the mean daily streamflow.

Moving average (ma) is an exponentially smoothed moving average that attempts to incorporate the influence of prior streamflow events on concentration at the current time period (Wang, Kuhnert & Henderson, 2011; Kuhnert et al., 2012; Wang & Tian, 2013; Zhang & Ball, 2017), using Eq. (2): (2)${\displaystyle ma\left(d\right)=\frac{Q_j+d{Q}_{j-1}+\dots +d^{j-1}{Q}_1}{1+d+...+d^{j-1}}}$

where ma at discount factor d is calculated using flows (Q₁, Q₂, …, Q_j) for days 1 through j. Here, as d approaches zero, the time series becomes the daily observed values and as d gets closer to one, the time series becomes the mean flow. Although Kuhnert et al. (2012) and Wang & Tian (2013) suggest using multiple covariates with d set at different values from 0.75 to 0.99, Robson & Dourdet (2015) and Zhang & Ball (2017) found substantial improvements in statistical nutrient models with d equal to 0.95 which was adopted in this study.

Flow anomaly (fa) is a dimensionless term that represents how wet or dry (or how anomalous) the current time period is from previous time periods (Vecchia et al., 2009; Zhang & Ball, 2017). Long-term flow anomaly (ltfa) is the streamflow over the previous year relative to the entire period (Eq. 3, Zhang & Ball, 2017) and and the short-term flow anomaly (stfa) calculated as the current day flow compared to the preceding 1-month streamflow (Eq. 4, Zhang & Ball, 2017):

(3)${\displaystyle ltfa\left(t\right)={\bar{x}}_{1year}\left(t\right)-{\bar{x}}_{entireperiod},}$ (4)${\displaystyle stfa\left(t\right)=x_{currentday}\left(t\right)-{\bar{x}}_{1month}\left(t\right)}$

where x are the averages of log-transformed streamflow over the antecedent period (1-year, 1-month, etc.) for time t. We used ltfa in NO₃ models and stfa in TP models based on previous work demonstrating major improvements in NO_X regression models that incorporated ltfa and moderate improvements in TP regression models that incorporated stfa (Zhang & Ball, 2017). Moving averages and flow anomalies were calculated with the adc R package (Schramm, 2023).

The calculation of model terms for the Lake Texana site were modified because daily loads are not a function of natural stream flow processes alone, but of dam releases and nutrient concentrations at the discharge point of the lake. Q, ma, and fa terms were calculated based on total gaged inflow from the 4 major tributaries to the lake. Thin-plate regression splines were used for ddate, log(1+Q), fa, and ma. A cyclic cubic regression spline was used for yday to ensure the ends of the spline match (day 1 and day 366 are expected to match). First order penalties were applied to the smooths of flow-based variables which penalize departures from a flat function to help constrain extrapolations for high flow measurements.

Left-censored data were not uncommon in this dataset (4% of TP samples and 20% of NO₃ samples). Several methods are available to account for censored data. We transformed left-censored nutrient concentrations to one-half the detection limit. Although this simple approach can introduce bias (Hornung & Reed, 1990), we considered it acceptable because high concentrations and loadings are associated with high-flow events and low-flow/low-concentration events will account for a small proportion of total loadings (McDowell et al., 2021).

Daily loads were estimated as the predicted concentration multiplied by the daily streamflow. For the Navidad River (USGS-08164525) site, daily loads at the dam were calculated from the discrete daily concentration at the discharge point of the lake and corresponding reported daily discharge from the dam. Flow-normalized loads were estimated similar to WRTDS by setting flow-based covariates on each day of the year equal to each of the historical values for that day of the year over the study period (Hirsch, Moyer & Archfield, 2010). The flow-normalized estimate was calculated as the mean of all the predictions for each day considering all possible flow values. Standard deviations and 90% credible intervals were obtained by drawing samples from the multivariate normal posterior distribution of the fitted GAM (Wood, 2006; Marra & Wood, 2012; McDowell et al., 2021). GAM performance was evaluated with repeated 5-fold cross validation (Burman, 1989) using average Nash-Sutcliffe Efficiency (NSE), Pearson sample correlation (r) and percent bias (PBIAS) metrics.

Linking estuary water quality to hydrology and nutrient loads

To test if changes in freshwater inflow and nutrient loading had explanatory effect on changes in estuary water quality a series of GAM models were fit at each site relating parameter concentration to temporal trends (Eq. 5), temporal trends and inflow (Eq. 6), and temporal trends, inflow, and nutrient loads (Eq. 7): (5)${\displaystyle g\left(\mu \right)=\alpha +f_1\left(ddate\right)+f_2\left(yday\right),}$ (6)${\displaystyle g\left(\mu \right)=\alpha +f_1\left(ddate\right)+f_2\left(yday\right)+f_3\left(Q\right),}$ (7)${\displaystyle g\left(\mu \right)=\alpha +f_1\left(ddate\right)+f_2\left(yday\right)+f_3\left(Q\right)+f_4\left(Load\right)}$

where μ is the conditional expected response (parameter concentration), g() is the log link, and response variable was modeled as Gamma distributed with mean μ and scale λ. f₁(ddate) is decimal date smoothed with a thin-plate regression spline, f₂(yday) is the numeric day of year smoothed with a cyclic cubic regression spline. Prior work has shown that many water quality parameters may have lagged effects lasting days or even months following storms and large discharge events (Mooney & McClelland, 2012; Wetz & Yoskowitz, 2013; Bukaveckas et al., 2020; Walker et al., 2021). Walker et al. (2021) showed the impact of large flow events like Hurricane Harvey had somewhat short water quality impacts in Lavaca Bay and adjacent estuary systems, ranging somewhere around 20 days to a few months. We incorporate this lag effect by using cumulative totals for f₃(Q) and f₄(Load). f₃(Q) is the cumulative of the previous 20 days of inflow (combined measurements from Lavaca River and Navidad River) and f₄(Load) is the cumulative of the previous 20 days total NO₃ or TP watershed load, both smoothed with thin plate regression splines. The set of models specified for each water quality response are in Table 3.

Because streamflow and nutrient loads are tightly correlated, freshwater inflow can mask signals due to changes in nutrient loads alone. Following the methodology implemented by Murphy et al. (2022), both freshwater inflow and nutrient loads were prepossessed to account for season and streamflow respectively. Cumulative inflow values were replaced by seasonally adjusted cumulative inflow obtained from the residuals of a GAM model fit between season (day of year) and log transformed cumulative inflow. Nutrient loads utilized the flow-normalized loads estimated in the previous section.

This study used an information theoretic approach to evaluate evidence of model covariate effects on Lavaca Bay water quality. Model probabilities were calculated and compared using the AIC_c scores between each group of temporal, inflow, and inflow+load models (Burnham, Anderson & Huyvaert, 2011). Improvements in model probabilities provide evidence that the terms explain additional variation in the response variable.

Watershed nutrient loads

Predictive performance of nutrient loads ranged from “satisfactory” to “very good” based on standardized evaluation metrics of NSE, r, and PBIAS (Moriasi et al., 2015) calculated using 5-fold cross validation. Median goodness-of-fit metrics for NO₃ models in the Lavaca River were 0.34 NSE, 0.70 r, and 2.00 PBIAS. Navidad River NO₃ models appeared to perform slightly better with 0.48 NSE and 0.87 r but with higher bias at 10.90 PBIAS. Generally, TP models performed better than NO₃ models. Median goodness-of-fit metrics for TP in the Lavaca River were 0.81 NSE, 0.93 r, and −7.20 PBIAS. Navidad River TP models had similar performance with 0.91 NSE, 0.99 r, and −3.30 PBIAS. Density plots of metrics show similar distribution of values between sites for the same parameter, with the exception r values for NO₃ loads where Lavaca River had a much larger variance in values compared to the Navidad River (Fig. 2). TP GAMS had higher average NSE and r values and lower variance in metric values compared to NO₃.

Annual NO₃ and TP loads show considerable variation, generally following patterns in discharge (Figs. 3, 4). Flow-normalized TP loads at both sites and flow-normalized NO₃ loads in the Lavaca River indicated watershed based loads did not change much over time when accounting for variation driven by streamflow (Fig. 3). Flow-normalized loads in the Lavaca River showed small variation over time with some decreases in NO₃ loads since 2013.

Aggregated across both sites, the mean annual NO₃ load from 2005 through 2020 was 205,405 kg (126,867 kg–341,569 kg, 90% CI). Annual NO₃ loads ranged from 12,574 kg in 2011 to 794,510 kg in 2007. Total annual TP loads ranged from 7,839 kg in 2011 to 595,075 kg in 2007. Mean annual TP loading from 2005 through 2020 was 182,673 kg (152,227 kg–219,310 kg, 90% CI). On average, the Navidad River accounted for 68% of NO₃ loads and 59% of TP loads from 2005 through 2020. However, during periods of extreme drought the Lavaca River became the primary source of nutrient loading in the watershed with the Navidad River only accounting for 15% and 25% of NO₃ and TP loads in 2011 (Fig. 4).

Linkages between water quality and watershed flows and loads

There was no evidence of long-term changes in TP or DO concentrations at any Lavaca Bay site (Fig. 5). The upper-bay site(TCEQ-13563) had evidence of a long-term linear increase in NO_X while chlorophyll-a decreased from 2005 through 2014 (Fig. 5). NO_X concentration at the mid-bay site (TCEQ-13383) displayed a periodic pattern that is indicative of a strong influence from inflow or precipitation. The temporal GAMs did not provide evidence of long-term trends in any of the water quality constituents at the lower-bay site (TCEQ-13384).

Freshwater inflow provided additional explanation for changes in NO_x and chlorophyll-a at all of the Lavaca Bay sites according to AIC_c and model probability values (Table 4). Freshwater inflow also explained additional variation in TP at the upper- and mid-bay sites but not the lower-bay site. Freshwater inflow did not explain additional variation in TKN concentrations at any of the sites. The upper- and lower-bay sites saw improvements in explanation of DO concentration with the inclusion of inflow.

Inclusion of TP loads provided additional explanation of TP concentrations at the upper- and mid-bay sites. Inclusion of NO₃ loads did not provide model explanatory improvements at any of the sites. The addition of nutrient loads (both TP and NO₃) terms did not provide additional explanation for changes in DO concentrations but did provide model improvement for the chlorophyll-a model at the mid-bay site.

Increases in aggregated freshwater inflow resulted in increases in TP at the upper- and mid-bay sites and increases in NO_X and chlorophyll-a concentration at all three sites (Fig. 6). Linear increases in DO concentration were observed with increasing flow at the upper- and lower-bay sites. TKN showed no response to changes in inflow at any of the sites.

Increased TP loads resulted in nearly linear increases of TP concentration at the upper- and mid-bay sites (Fig. 7). The relative effect size appeared much smaller than the effect of freshwater inflow alone. Increased NO₃ loads did effect NO_x concentrations at any site. Linear increases in chlorophyll-a were observed in response to increased TP loads at the mid-bay site, but not the other sites or in response to NO₃ loads.

Nutrient loads

TP and NO₃ loadings from the Lavaca Bay watershed showed high inter-annual variability driven primarily by fluctuations in discharge. Notably, there were no indications of trends in flow-normalized NO₃ and TP loads in the Navidad River. In comparison, there was weak evidence of more recent decreases in flow-normalized NO₃ (but not TP) loads in the Lavaca River watershed. While the dominant agricultural land uses differ between the Lavaca (primarily grazed pasture and rangeland) and Navidad(mix of pasture and row crops) catchments, we did not have a reason to expect different flow normalized trends between the two systems from land use alone. Freshwater discharges in the Navidad River are regulated by the Palmetto Bend Dam forming Lake Texana at the lower extent of the river. Lentic nitrogen uptake and cycling may have regulating effects that mask changes in upstream nitrogen loadings. Additional nutrient data collection in the tributaries of Lake Texana is needed to fully assess the role of Lake Texana in regulating nutrient delivery to the Lavaca Bay system. However, these results suggest that there have been no changes in the NO₃ or TP loading from the Navidad River system at the Lake Texana discharge point when accounting for variations in year to year discharge.

The evidence of decreased Lavaca River NO₃ loading, although weak, is a potential positive sign for water quality managers working to implement practices that improve water quality in the freshwater sections of the Lavaca River watershed. Planning and implementation efforts to increase agricultural producer participation in water quality protection practices have been ongoing in the watershed since 2016 (Schramm et al., 2018; Berthold, Olsovsky & Schramm, 2021), however little work has been conducted to directly link these efforts with water quality outcomes. The decrease in flow-normalized NO₃ loads could be a reflection of those collective efforts but the lack of evidence for similar changes in flow-normalized TP loads provide contrary support. The inconsistent flow-normalized trends may also reflect some of the weakness of the water quality dataset that is primarily composed of ambient water quality measurements. The issues associated with the lack of flow-biased measurements is further discussed later in this section.

Some prior studies have generated estimates of mean annual TP yields in the Lavaca River watershed (Table 5, Dunn, 1996; Rebich et al., 2011; Omani, Srinivasan & Lee, 2014; Wise, Anning & Miller, 2019). Although these studies differ in time periods and methodologies, they provide a sanity check for the reasonableness of the annual estimates generated in the current study. In a regional assessment of nutrient loading in river’s along the Gulf of Mexico, Dunn (1996) used the LOADEST model to develop an estimated mean annual yield of 28.9 kg/km². LOADEST is a multiple linear regression model that fits log transformed pollutant concentrations to long term, seasonal, and flow based predictors and includes methods for bias correction when exponentiating the response variable. Rebich et al. (2011) and Wise, Anning & Miller (2019) used SPARROW to provide a more recent assessment of regional catchment based loadings to the Gulf of Mexico(Table 5). SPARROW is a hybrid statistical-process model with the underlying nutrient load estimation methods based on the previously described LOADEST (Schwarz et al., 2006). The functional form of the LOADEST regression model is similar to the terms applied in the GAMs used in the current study. The only study to apply a mechanistic watershed model (SWAT) to estimate nutrient loadings in the Lavaca River watershed Omani, Srinivasan & Lee (2014) developed estimated yields (42 kg/km²) similar to the two SPARROW models. Although direct comparisons are complicated by varying time periods, the estimates in this study do fall within the range of of previously developed estimates. To evaluate changes in long-term trends in discharge might be associated with the nutrient yield estimates covering different time periods, we fit a GAM relating log-transformed daily discharges on the Lavaca River to season and time (Fig. 8). The long-term trends in discharge indicate watershed discharges were at or above average from 1972 though the early and mid-1980s. In comparison watershed discharges since the mid-2000s are at or below average. It is probable that the lower than average discharges observed from 2010 through 2021 (Fig. 8) bias our estimates downward compared to studies that included higher than average streamflow periods (1995–2005). Overall, the ranges of estimated yields among different studies along with the apparent large variability in streamflow driven loadings (Fig. 3, Fig. 4) suggest that the current estimates of TP loading are reasonable.

Estuary water quality

The non-linear estuary water quality trends identified in the current study differed slightly from previously identified trends (Bugica, Sterba-Boatwright & Wetz, 2020). This is not unexpected due to the different time periods, different methodology, and generally small slopes identified for most of the significant water quality parameters in prior work. Both DO and cholorophyll-a concentrations at all three Lavaca Bay sites were stable from 2005 through 2020. This is a positive outcome in comparison to other Texas estuaries that are facing larger demands for freshwater diversions, higher population growth, and more intense agricultural production which have resulted in more direct signs of eutrophication (Wetz et al., 2016; Bugica, Sterba-Boatwright & Wetz, 2020). Despite the stability of DO and cholorophyll-a, there are concerning site specific increases in NO_X and TKN concentration over the same time period. These trends are especially concerning due to the nitrogen limitation identified in many Texas estuaries (Gardner et al., 2006; Hou et al., 2012; Dorado et al., 2015; Paudel, Montagna & Adams, 2019; Wetz et al., 2017) and the relatively low ambient concentrations observed in Lavaca Bay.

The strong positive effect of freshwater inflow on NO_X and TP concentration are suggestive of nonpoint watershed sources, consistent with watershed uses and with other studies relating freshwater inflow with nutrient concentrations in Lavaca Bay and other estuaries (Russell, Montagna & Kalke, 2006; Caffrey et al., 2007; Peierls, Hall & Paerl, 2012; Palmer & Montagna, 2015; Cira, Palmer & Wetz, 2021). Inflow had a non-linear relationship with NO_x, with NO_x increasing as freshwater inflow transitioned from low to moderate levels. At higher freshwater inflows, the NO_x decreased at the mid- and lower-bay sites, possibly indicating a flushing effect at higher freshwater inflow. No relationship between TKN and freshwater inflow was observed at any of the sites. The results coincide with previous work that suggest processing of organic loads in the tidal portions of the Lavaca River reduce transport of nutrients to the lower reaches of Lavaca Bay (Russell, Montagna & Kalke, 2006).

Freshwater inflow also displayed positive effects on chlorophyll-a at each site, with the largest effect at the lower-bay site. The lower-bay site also showed an attenuation in in chlorophyll-a at the highest freshwater inflow volumes. Freshwater flushing or increases in turbidity are associated with limitations or decreases in chlorophyll-a in other estuaries (Peierls, Hall & Paerl, 2012; Cloern, Foster & Kleckner, 2014). No relationships between inorganic nitrogen and chlorophyll-a were observed. However, a small positive effect between flow-normalized TP loads and chlorophyll-a concentrations were detected at the mid-bay site. Although, this is not suggestive of a phosphorus limitation in the Bay, it is supportive of prior work emphasizing the importance of controlling both nitrogen and phosphorus runoff to protect water quality (Conley et al., 2009). Due to the lack of TKN loading information, no assessment between organic nitrogen loads and chlorophyll-a were possible.

Although other studies have identified complex relationships between estuary nutrient concentrations, nutrient loading and chlorophyll-a concentrations in Texas estuaries (Örnólfsdóttir, Lumsden & Pinckney, 2004; Dorado et al., 2015; Cira, Palmer & Wetz, 2021; Tominack & Wetz, 2022), this study specifically used flow-adjusted freshwater derived nutrient loads to parse out contributions from changes in nutrient loadings while accounting for variations in load due to flow. Loading GAMs indicated no evidence of changes in flow-normalized TP loads in either river, and no changes in flow-normalized NO₃ loads in the Navidad River. The small changes in flow-normalized NO₃ loads in the Lavaca River are probably masked under most conditions by discharge from the Navidad River. Given the relatively small variation in flow-normalized loads, it can be expected that they would contribute little to the variance in downstream water quality.

There was no evidence that adjusted freshwater inflow and nutrient loads had effects on DO concentration in Lavaca Bay. The seasonality term in the temporal GAM models explained a substantial amount of DO variation at all of the sites. Responses of estuary metabolic processes and resulting DO concentrations can be quite complicated and often locally specific (Caffrey, 2004). While the lack of total nitrogen or TKN loading data hinders interpretation, the large seasonal effect on DO concentration indicates physical factors (such as temperature, wind, and turbidity) play an important role and should be included in future models. Prior work suggests that Lavaca Bay may not be limited by nutrients alone, with high turbidity or nutrient processing in upper portions of the Bay or intertidal river limiting production (Russell, Montagna & Kalke, 2006). Finally, it is reasonable to assume that fluctuations in DO may not occur immediately in response to nutrient pulses or freshwater inflow. Work has shown that many water quality parameters may have lagged effects lasting days or even months following storms and large discharge events (Mooney & McClelland, 2012; Wetz & Yoskowitz, 2013; Bukaveckas et al., 2020; Walker et al., 2021). However, this study only evaluated responses to 20-day cumulative loading and inflows.

Overall, this study suggests that DO and chlorophyll-a concentrations have been relatively stable in Lavaca Bay. Site-specific increases in TKN and NO_X concentrations are a cause of concern for increasing risks of eutrophication within Lavaca Bay which might be currently attenuated by changes in freshwater flow, turbidity, and other physical processes. While loading models indicate that there are large annual fluctuations in NO₃ loads, these changes have been largely driven environmental conditions (changes in runoff and river discharge). These models also provide evidence that estuary NO_X and TP concentrations are strongly driven by freshwater inflow and to a lesser extent fluctuation in flow-normalized riverine loadings. Site-specific changes in the relationships between freshwater inflow and responses in both chlorophyll-a and NO_x concentrations are indicative of nutrient processing and or tidal flushing effects moving from the river discharge point to the mouth of Lavaca Bay. This study does not completely explain site specific increases in NO_X and TKN concentrations in Lavaca Bay. The freshwater study sites did not quantify nutrient loadings from tidal contribution areas or ungauaged watersheds. Nutrient contributions from wastewater facility discharges, septic systems, and stormwater could be considerable contributors to nutrient loadings in Lavaca Bay since they are not processed by a tidal river reach prior to entering the Bay. The Garcitas-Arenosa Creek, Placedo Creek, and Cox Bay subwatersheds are currently undersampled but compose approximately 27% of the watershed area. The contribution of nutrient loadings from these undersampled areas is unknown.

Future management and research needs

The GAM approach proved useful for both estimating loads and assessing downstream responses in water quality. Although we did not compare other models, it is likely similar estimates of loadings would be obtained by methods such as LOADEST, WRTDS, or SPARROW given the functionally similar dependent variable structures. The underlying weakness in the estimates of loading in the current study is the reliance on ambient water quality data used for statewide water quality assessments. Cross-validation of the nutrient loading models highlights that predictions are prone to high bias, owing to the lack of targeted storm or flow biased measurements. The high biases are indicative that subsets of values were unable to capture the functional relationships with the flow based dependent variables. It was beyond the scope of the current study to evaluate the subsets of cross-validation data and scores. However, the cross-validation procedure is indicative that more robust sampling is needed. Supplementary flow-biased monitoring targeting storm- or high-flow conditions is critical to improve model performance and strength of evidence produced by these models (Horowitz, 2003; Snelder, McDowell & Fraser, 2017). Although there is existing work on the samples sizes and sample design required for reliable performance of both LOADEST (Park & Engel, 2014) and WRTDS (Kumar et al., 2019) models, similar work does not appear to have been extended to water quality applications of GAMs.

Due to the concerning increases in eutrophication associated parameters in Lavaca Bay and other Texas estuaries (Bugica, Sterba-Boatwright & Wetz, 2020), and the desire to quantify linkages between environmental outcomes and on the ground management actions (Schramm, Gitter & Gregory, 2022) there is a strong need for reliable estimates of pollutant loadings and responses along the Texas coast. Within Texas, statewide water quality monitoring programs have focused on collection of ambient condition data. A framework for establishing pollutant load monitoring programs across catchments that explicitly incorporate flow biased data is needed for assessing nutrient loading and estuary health along the Texas coast.

Additional efforts focused on identifying relevant effect sizes, sampling designs, and funding mechanisms that can support long term efforts are also needed to adequately design such a framework. Large long-term monitoring programs in and around the Chesapeake Bay, San Francisco Bay, and along the Mississippi River have proven extremely effective at informing management actions and tracking progress towards long-term pollutant reduction goals. Similar coordinated efforts across Texas coastal watersheds would prove useful for resource management efforts intended to protect the biological and water quality integrity of Texas’s estuaries.