Identification and Analysis of Storm Events

Storm events were extracted from about 10 years of 15-minute stage data at 14 monitoring stations to evaluate hydrologic, nutrient, and sediment conditions during storms. This effort expands on work by Porter and others (2020b), in which hydrologic and water-quality responses were characterized during all flow conditions. Streamflow data are commonly used to extract stormflow hydrographs for the purpose of event-based analyses (Blume and others, 2007; Hopkins and others, 2020), and such methods are well established (Eckhardt, 2005; Nathan and McMahon, 1990; Sloto and Crouse, 1996); however, stage data, which measure water level above an arbitrary datum, were used in this report because continuous streamflow data were not available at all stations.

Before extracting storm events, data were adjusted to correct two issues unique to stage data that would have reduced the accuracy of results: shifting channel geomorphology and diel patterns evident during low flows in some watersheds. The first issue results from changes in channel geomorphology that are common in small urban watersheds and can cause gradual drifts or, after storms, sudden shifts in stage that may not reflect changes in streamflow. Records were corrected for such channel instability by subtracting the i^th observation from a 5-day sliding minimum value. The second issue, diel waves in the dry weather stage data, was likely caused by variability in barometric pressures and (or) effects of evapotranspiration (Bond and others, 2002; Schwab and others, 2016). Because these fluctuations were not indicative of storm conditions, the diurnal signal was removed by applying a smoothing filter that replaced the i^th value with an average of the previous 24 hours, representing 96-unit values. These adjusted data were used only for storm extractions, not for subsequent calculations of storm-event metrics.

Storm events were first extracted from the adjusted data as periods when stage exceeded a local minimum, a technique based on approaches used by Nathan and McMahon (1990). The following modifications were made to this published method: minima of 3-day nonoverlapping windows were computed for the entire study period and values less than 1.05 times the outer values (representing values between each minima) were located to serve as turning points. The timeseries was then computed by connecting all minima values. Next, storm events with a maximum stage range less than 0.5 foot were removed because these periods typically represented non-storm related fluctuations in water level. Finally, storms with multiple peaks were removed to avoid bias in several event-based metrics.

Multiple metrics were calculated to represent the number, duration, hydrologic characteristics, and water-quality conditions of extracted storm events. Holistically, these metrics provide insight into how each watershed transforms precipitation into event-based runoff. Rising- and falling-limb duration was calculated as the time from the beginning of a storm to peak stage, and from peak stage to the end of a storm, respectively. Measures of duration were the only characterization of storms made from stage data at all 14 stations. Three additional measures of storm-event hydrology were calculated using streamflow data at four intensive-monitoring stations. The magnitude of each event (“peak”) describes the peak streamflow volume normalized by watershed area. The volume of stormflow (“volume”) describes the area-normalized amount of total streamflow delivered during the event. The rate of streamflow change during the rising limb (“rise rate”) describes peak streamflow divided by the rising duration. Event-mean concentrations (EMCs) of TN, TP, and SS were calculated for each storm event at the four intensive-monitoring stations as the storm-event constituent load divided by the total storm-event streamflow volume. Unlike constituent loads, which are a function of concentration and runoff volume, EMCs represent the average flow-weighted water-quality concentrations during a storm event. These metrics commonly are used to assess the quality of stormwater runoff, exceedance of water-quality standards, and the efficacy of management-practice treatment and prevention processes (Maniquiz and others, 2010; U.S. Environmental Protection Agency, 2002).

Trends over time in all storm-event metrics were evaluated by a Mann-Kendall test, and the magnitude of trend was summarized from the slope of a Theil-Sen line, a non-parametric approach that is commonly used in environmental statistics and was used in previous Fairfax County water-quality analyses (Helsel and others, 2020; Porter and others, 2020b). Trends in the number and duration of extracted storm events represent conditions from 14 stations from August 2009 through March 2018, whereas trends in storm-event metrics derived from streamflow and water-quality data represent conditions from four intensive-monitoring stations from water years 2009 through 2017. Trends were not adjusted for precipitation because spatially explicit, sub-daily precipitation data were not available to this study. All analyses were conducted in R v.4.0.3 (R Core Team, 2020).

Computation and Analysis of Water-Quality Loads

Nutrient and SS loads describe the total mass of material delivered downstream from a monitoring station over a period of time. Annual loads are often of interest because they are commonly used for regulatory purposes (Fairfax County, 2017) and to quantify expected nutrient and sediment management-practice effects. TN, TP, and SS loads were reported for the five intensive-monitoring stations from 2009 through 2018 using published surrogate regression models (Porter and others, 2020b) following previously documented methods (Jastram and others, 2009; Jastram, 2014; Rasmussen and others, 2008). These models quantify a fixed relation between observed nutrient and sediment concentrations and measurements available at sub-daily time steps, such as streamflow, WT, TB, and time of year. After estimating a continuous time-series of concentrations, loads were calculated by multiplying modeled concentrations by measured streamflow. Additional analyses based on these loads are described below.

Base-flow nitrogen, phosphorus, and sediment loads describe the amount of material transported past a station during non-storm-event hydrologic conditions. Many methods can be used to identify base flow (Raffensperger and others, 2017); the general intent of these methods is to separate the amount of streamflow associated with groundwater discharge from surface-water runoff (referred to as quickflow). This report used the “BaseflowSeparation” function within the EcoHydRology R package to define base flow and quickflow; a method used in related urban-stormwater research (Hopkins and others, 2020). This method uses an automated approach to identify base flow from a station hydrograph, with the separation of base flow and quickflow dependent on a digital filtering parameter (Duncan, 2019). A three-pass filtering parameter of 0.975 was selected for all stations after visually inspecting hydrograph separation results from a range of recommended values. Base-flow nutrient and sediment loads were then calculated by multiplying modeled concentrations by the proportion of streamflow delivered during base flow.

Trends in nutrient and sediment loads over time were calculated to describe how water-quality conditions changed from 2009 through 2018 in four intensive-monitoring stations. Nutrient and sediment loads in Fairfax County streams are heavily affected by hydrologic conditions (Jastram, 2014; Porter and others, 2020b); therefore, trends in load were calculated by a non-parametric method that removes variability associated with streamflow (Helsel and others, 2020). The resulting trends are likely more closely associated with effects from landscape activities, including the role of management practices. For each station, monthly nutrient and sediment loads were computed using surrogate regression models and regressed against total monthly streamflow using a LOESS smooth relation. Trends in the residuals of these relations were calculated with a Mann-Kendall test, which describes patterns in monthly load after adjusting for streamflow. Trend magnitude was determined from the slope of a Thiel-Sen line, and monthly residual values were retransformed into units of mass. These results are referred to as “flow-adjusted loads” throughout the report. A non-parametric test recommended for water-quality datasets was used to determine the timing and significance of a change point in the distribution of monthly flow-adjusted loads (Helsel and others, 2020; Pettitt, 1979).

Nitrogen and phosphorus retention rates describe the percentage of nutrient inputs in a watershed that are not delivered to streams. Annual retention rates of nitrogen and phosphorus were calculated for the intensive-monitoring stations as the difference between estimated inputs, reflecting the sum of all nonpoint sources and calculated loads. Retention rates are low when most nutrient inputs in a watershed are delivered to streams, which typically occurs during high-streamflow years, as these conditions generate large loads. Therefore, changes in nutrient retention over time are explored as a function of computed loads and flow-adjusted loads.

Linear Mixed-Effect Models

Spatial and temporal relations between water-quality responses and watershed predictors were evaluated using linear mixed-effect (LME) models. In particular, LME models permitted management-practice effects to be evaluated after accounting for climatic and landscape factors: a critical step to identify the effects of management practices on water-quality and ecological responses. Linear mixed-effect models, also known as multilevel models or hierarchical models, are an extension of ordinary least-squares regression and are useful for datasets with a nested structure, where data within groups are more similar than data between groups. Such datasets violate ordinary least-squares assumptions of independence and, without accounting for their nested structure, other models can produce Type I errors (incorrect rejection of a null hypothesis) and biased parameter estimates (Peugh, 2010). LME models allow researchers to explore the effect of variables on differences within groups (referred to as level-1 variance) and on differences between groups (referred to as level-2 variance). Longitudinal studies, where repeat observations are collected from multiple groups over time, are often analyzed using LME models to understand changes within and between groups (Grimm and others, 2016). Similar methods have been used to characterize conditions affecting biological communities (Wenger and others, 2016; Xu and others, 2010) and water-quality conditions (Flint and McDowell, 2015; Funes and others, 2018; Goyette and others, 2019; Li and others, 2013).

LME models were used to evaluate how annual predictor variables related to annual measures of TN, TP, SS, SC, and benthic-macroinvertebrate IBI scores collected from 14 watersheds over a 10-year period. Watersheds were used to define the grouping structure of each LME model. For TN, TP, SS, and SC, relations between predictors and responses were evaluated from 2009 through 2018 using median-annual values summarized from monthly samples. Median-annual values allow for joint consideration of data collected across 14 watersheds, a comparison against annual predictor data, and a reduction of the influence of samples collected during stormflow conditions. For IBI scores, models were built using 132 samples collected annually across 14 watersheds from 2008 through 2017. As samples were collected during March or April of each year, models were built to relate annual IBI scores with predictors representing conditions in the year before sample collection. For landscape predictors in 2007, the 2008 value was used under the assumption that no change occurred over the 2007–08 period. All models were built using original units of measurement, that is, mg/L for TN, TP, and SS; microsiemens per centimeter (µS/cm) at 25 degrees Celsius (°C) for SC; and points for IBI scores.

LME models, described in equations 1–6, were used in this study to (1) partition the overall variability of each response variable to within-watershed and between-watershed variability, and (2) account for both sources of variability to the maximum extent possible using explanatory data. For example, TN concentrations from 2009 through 2018 varied annually within a watershed (the variability of observations about a watershed-specific regression line), and regression intercepts and slopes differed between watersheds (fig. 3). A model for this partitioning, described in equations 1–3, is referred to as a linear-growth model (Grimm and others, 2016). The “level-1” equation (eq. 1) in the model specifies a watershed-specific intercept ($b_{1i}$, eq. 1) as the sum of two “level-2” components: an overall fixed-effect constant (${\beta }_{01},$eq.2)and a watershed-specific random effect ($d_{1i},$eq. 2) watershed-specific slopes ($b_{1i}$, eq. 1) were specified similarly (${\beta }_{02} \text{and} d_{2i},$eq. 3). The fixed effects ${\beta }_{01}$ and ${\beta }_{02}$ can be interpreted as the overall average response value at time zero and the overall average linear rate-of-change over time, respectively, whereas the random effects reflect watershed-specific deviations from the corresponding fixed effect. The linear-growth model in equations 1–3 includes six estimated parameters: estimates of the fixed-effect intercept (${\beta }_{01})$ and slope (${\beta }_{02}$), variance of the random intercept ($d_{1i}$) and slope ($d_{2i}$), variance of the residuals ($u_{ti}$), and the covariance between the random intercept and slope. Collectively, the random intercept, random slope, and residual variance in the linear-growth model (eqs. 1–3) represent the unaccounted-for variability of all observations about the fixed-effect population intercept and slope.

The second objective of using LME models, as outlined in the previous paragraph, is achieved by augmenting the linear-growth model with explanatory covariates that reduce one or more of these three variance components in a manner that improves, or at least does not degrade, the fit to the observed data. Equations 4–6 describe a linear-growth model with capacity to include predictors as fixed effects that can vary over time, space, or both (table 2, table 3, table 5, and table 6), while retaining the same random-effect structure as the linear-growth model with no covariates, as described in equations 1–3. Time-invariant predictors, for example watershed area, have no influence over the variability of individual observations about a watershed-specific regression line; they only can affect between-watershed differences in slope and intercept. Thus, such variables and their corresponding coefficients (notated as $X_1$ though $X_c$ and ${\beta }_{1x}$ through ${\beta }_{cx}$, respectively) can only appear in equations 5 (intercept; x=1) and 6 (slope; x=2), which are level-2 equations; the benefit of their inclusion is evidenced by the reduction in random slope and intercept variances. Time-varying predictors vary over time and, potentially, space. Those that vary over time but not space can affect within-watershed variability over time but not between-watershed differences in slope and intercept. These variables, which exclusively represent county-wide climatic predictors in this study, appear only in equation 4 (level-1), notated as $Z_1-Z_{\text{f}}$ with corresponding coefficients $b_3-b_{f+2}$. The benefit of the inclusion of time-varying predictors would be evidenced in a reduction of residual variance in the models.

Several predictors vary over space and time. For example, the percentage of impervious cover has generally increased in the Captain Hickory Run at Route 681 near Great Falls, Va. (CAPT HICK) and Castle Creek at Newman Road at Clifton, Va. (CASTLE) watersheds from 2009 through 2018, but impervious cover in CAPT HICK remained consistently higher than in CASTLE over this 10-year period (fig. 4). To isolate their effects on random variance components, such predictors were decomposed following recommended centering techniques (Enders and Tofighi, 2007; Grimm and others, 2016; Hoffman and Stawski, 2009). In brief, raw predictor values were split into two separate components. First, time-invariant average levels of the predictor were computed for each watershed, then centered relative to the overall population mean, such that the mean of the watershed averages was zero. Second, the watershed-specific values of the predictor over time were centered relative to each watershed’s mean value, such that the average in each watershed was zero (fig. 4); hereafter, the former component is referred to as “grand-mean” centered, and the latter as “watershed-mean” centered. Time-invariant predictors were centered only to their grand-mean; time-varying predictors were centered to their grand-mean and watershed-mean. Grand-mean-centered components of time-invariant predictor and time-varying predictor variables, included through the level-2 slope and intercept equations (eqs. 5 and 6), have the capacity to reduce random intercept and slope variances. Watershed-mean time-varying predictor components (eq. 4) have the capacity to reduce residual variance.

The overarching goal of modeling water-quality responses with these additional parameters is to identify relations among watershed predictors and response variables. Specifically, models were constructed to address the following three relations: (1) predictors and average responses between watersheds, (2) predictors and the linear rate-of-change in responses between watersheds over the study period, and (3) annual changes in predictors and responses within watersheds over the study period. Except where the random-effect variance-covariance matrix converged to a perfectly linear relation (referred to as a singularity error in LME-model literature) and required simplification, all models were fit following recommendations advising the use of a maximal random-effect structure (Barr and others, 2013) and used a random intercept and random-linear relation with time for each watershed.

A combination of predictor variables was selected for consideration in LME models, largely guided by the evaluation of individual predictor-response relations and by applying professional judgement guided by previous Fairfax County studies (Hyer and others, 2016; Jastram, 2014; Porter and others, 2020b). After appropriate centering and decomposition, all possible combinations of selected predictors were modeled. Models that resulted in the largest variance reduction with the fewest additional parameters were prioritized because the linear-growth model described in equations 1–3 already contains six parameters for a response dataset of about 140 observations. Likelihood-ratio tests and model Akaike information criterion (AIC; Akaike, 1974) were used in model comparisons to determine the value of an added predictor. The usefulness of a predictor was also evaluated by examining the 95-percent confidence intervals of a coefficient to ensure that these ranges do not include zero. For all selected models, standardized coefficients are reported following methods described in Hoffman (2015) to determine the relative-effect size of predictors on modeled responses. Conditional and marginal R² values are reported for selected models and respectively describe the proportion of total variance explained by fixed and random effects and by only fixed effects (Nakagawa and others, 2017). Tables used to evaluate relations between predictor and response variables and guide the selection of LME models are contained in appendix 1. Tables and figures that support the interpretation of selected and alternative LME models are contained in appendix 2.

Uncertainties and limitations of the predictor and response datasets affect the identification and interpretation of predictor-response relations. Except for climatic and water-quality data, most predictor data represent little change or a linear rate-of-change over the 10-year period; therefore, they do not help to explain the interannual differences in responses within watersheds. Therefore, temperature, precipitation, and water-quality terms, such as WT and DO, may be preferentially selected in LME models for their unique ability to reduce within-watershed variances. Causal links between response variables and annual measures of climate and water quality may exist, but the inclusion of these terms might be affected by the limited temporal accuracy of other predictor datasets.

Limitations of the LME model analysis also affect the identification and interpretation of predictor-response relations. This study did not consider every possible combination of predictors, predictor interactive effects, alternative random-effect structures, or models that may offer additional insights. LME models evaluate the strength of linear predictor-response relations and likely underrepresent non-linear patterns that may be present in the data. Strong correlations among predictors challenge the ability of LME models to identify the most influential predictors of a given response. It is convenient to discuss the interpretation of a single model; however, alternative models with similar performance often exist and are acknowledged in this report. Finally, the legacy effects of a predictor may heavily influence contemporary water-quality responses, but, generally, were not considered.

All LME model development was performed in R v.4.0.3 (R Core Team, 2020) using the “lme4” package (Bates and others, 2015). Standardized model coefficients were produced with the “standardize_parameters” function in the “effectsize” package (Ben-Shachar and others, 2020).

Patterns of Observed Responses

TN concentrations differed among study watersheds, varied seasonally, and were strongly related to streamflow (Porter and others, 2020b). Median-annual TN concentrations below about 2 mg/L typically occurred in most watersheds, but concentrations above 3 mg/L were observed in CAPT HICK and SFLIL (fig. 34). More than 90 percent of the TN in the monthly samples was typically dissolved nitrogen, mostly nitrate. Concentrations of TN and nitrate were commonly lower in warmer months than cooler months, likely reflecting nitrate loss through plant uptake and denitrification. TN concentrations typically increased with streamflow, whereas nitrate concentrations decreased (fig. 35). Higher stormflow TN concentrations likely result from surface-water runoff that carries leaves, grass clippings, woody debris, dry atmospheric nitrogen deposition, and other materials containing particulate forms of nitrogen to streams.

Nitrogen-concentration trends were evaluated using monthly water-quality samples collected from April 2008 through March 2018 to describe how conditions changed over time in 14 watersheds with 10 years of data, analyses, and interpretations originally conducted by Porter and others (2020b; table 20). After adjusting for streamflow, TN concentrations declined in 5, increased in 1, and had no trend in the remaining 8 watersheds. No network-wide trend in TN was detected with a regional Seasonal-Kendall analysis. TN decreases ranged from about −2 percent to −4 percent per year; the largest 10-year decline of about 0.5 mg/L was in DEAD. The TN increase in SFLIL represented a change of about 2 percent per year, or an increase of about 0.6 mg/L. The increase in this watershed is notable as SFLIL had higher average nitrogen concentrations relative to the rest of the monitoring network. At all stations, TN trends were greatly affected by nitrate concentrations. Nitrate decreased network-wide by -1.3 percent per year. Flow-adjusted nitrate concentrations increased only in SFLIL and CAPT HICK (by 3 percent and 1 percent per year, respectively): the two watersheds with the highest average nitrate concentrations. Nitrate decreases occurred in four watersheds; the largest percentage of decrease was −6 percent in BRR, and the largest total decrease was −0.6 mg/L in DEAD.

Between-watershed differences in median-annual TN concentrations were typically larger than year-to-year changes within watersheds (fig. 34). These two sources of variability, representing spatial differences among study watersheds and temporal differences from 2009 through 2018, were described by a linear-growth model (table 21). This model indicated that the average TN concentration across all watersheds in 2009 was 1.92 mg/L and that the linear rate-of-change could not be differentiated from zero. Random-effect variances describe a dataset where group-level variance (0.929), representing between-watershed differences in intercept, accounted for most total model variance (0.961). The variances representing TN change over time, differences in between-watershed slope (0.001) and within-watershed annual changes (0.031), were small. The median-annual TN changes that occurred from 2009 through 2018 included TN decreases in DEAD and increases in SFLIL (fig. 36). The AIC of the linear-growth model (12.06) was compared against models that included additional predictors; more parsimonious models reduced the random-effect variances.

Relations to Predictor Variables

Between-watershed TN differences were likely affected by a combination of human and natural factors, patterns that were explored by comparing average predictor values to TN concentrations from 2014 through 2018 at 20 study watersheds (app. 1, table 1.1). Septic-system density (SEPTIC DEN, r=0.74; fig. 37A), the estimated nitrogen input from septic systems (N SEPTIC, r=0.74), the sum of all nitrogen inputs (N TOT IN, r=0.66) and the length of sanitary-sewer lines (SAN SEWER LEN, r=−0.62) had the strongest correlations with TN, suggesting that between-watershed TN differences may be related to wastewater infrastructure. Septic-system effluent is likely an important factor in some watersheds with the highest TN concentrations, an influence previously discussed in this report, in previous Fairfax County studies (Hyer and others, 2016; Porter and others, 2020b), and in related water-quality investigations (Burns and others, 2005; Kaushal and others, 2006; Kaushal and others, 2011; Valiela and others, 1997; Wollheim and others, 2005). Negative correlations between TN and predictors such as impervious cover (IMP, r=−0.48; fig. 37B), developed land use (DEV ALL, r=−0.47), roadway length (TOT ROAD, r=−0.54), housing units (TOT HU, r=−0.57), and the amount of stormwater infrastructure (PIPE LEN, r=−0.52) demonstrate that more urbanized watersheds had lower average TN concentrations. Such conditions may occur because these watersheds contained the largest amount of infrastructure to control stormwater runoff, had the smallest fraction of wastewater served by septic systems, or otherwise limited the sources or delivery of nitrogen to streams.

Denitrification, the permanent removal of nitrogen through conversion to atmospheric gas, can represent a dominant nitrogen loss term in mid-Atlantic watersheds (Van Breemen and others, 2002) and may affect between-watershed differences in TN concentrations. Considering this process, the inverse relation between TN and urbanization seems unusual, because more impervious watersheds commonly deliver increased amounts of water through surface runoff and limit opportunities for denitrification and other microbial pathways that affect nitrogen availability and mobility (Burgin and Hamilton, 2007). Similar rates of denitrification can occur, however, in urban areas as in lesser developed watersheds; urban denitrification rates are strongly affected by near-stream environments such as floodplains and detention basins (Groffman and Crawford, 2003). Therefore, even in highly impervious watersheds, denitrification opportunities exist if riparian zones allow for infiltration and hyporheic exchange, an important process in Fairfax County where many stream valleys are designated as park land and are afforded protection from urban development. Denitrification potential is often greatest in areas with saturated soils, abundant organic matter, and warmer stream and air temperatures (Bowden, 1987; Burgin and Hamilton, 2007; Pilegaard, 2013). Correlations among predictors representing these features and median TN concentrations offered mixed evidence for the role of denitrification (app. 1, table 1.1). These patterns cannot be used to discount the importance of denitrification because denitrification rates are highly variable in urban areas (Groffman and Crawford, 2003) and predictor data used to summarize soil properties do not represent spatial variability within a watershed or throughout soil horizons.

Some temporally averaged measures of wastewater and urbanization were related to the 2009 through 2018 percentage change in median-annual TN concentrations (app. 1, table 1.2). Average impervious cover had the strongest correlation to TN percentage change (IMP, r=−0.84; fig. 38A), indicating that declines in TN concentration were larger in more impervious watersheds. Conversely, the correlation of TN to septic-system density (SEPTIC DEN, r=0.70; fig. 38B) showed that concentration increases were larger in watersheds primarily served by septic systems. Responses at OCSB and SFLIL demonstrated these patterns. OCSB, one of the most impervious watersheds, was served by a sanitary-sewer network and had a TN decline of about 3 percent per year. In contrast to these stations, SFLIL had a smaller amount of impervious cover and a high density of septic systems and had a TN increase of about 2 percent per year.

Correlations between annual measures of urbanization and median-annual TN concentrations from 2009 through 2018 varied by watershed; most predictors had no consistent relation to TN changes over time (app. 1, table 1.3). For example, impervious cover increased in OCSB and SFLIL from 2009 through 2018, but TN concentrations declined in OCSB (r=−0.71) and increased in SFLIL (r=0.68). Annual changes in the number of septic systems were not correlated with TN changes at most watersheds (SEPTIC DEN, median r=0.00); however, most watersheds were predominantly served by sanitary sewers, and septic-system increases may affect TN concentrations only in watersheds with elevated septic-system density. Median TN concentrations and the number of septic systems increased in SFLIL and CAPT HICK from 2009 through 2018, watersheds with hundreds of septic systems per square mile (table 8), resulting in positive correlations (r=0.63 and r=0.33, respectively). However, increased septic effluent contributed by new systems may not be the primary driver of increased TN concentrations in these watersheds because fewer than 30 new systems were installed in each watershed over the 10-year period. Not to discount the potentially important relations between septic systems and TN, it is possible that stream conditions were responding to past changes in septic infrastructure, nitrogen removal rates have declined in aging systems, or that effluent was delivered to streams more efficiently.

Annual changes in some spatially averaged climatic predictors were correlated with annual changes in median TN concentration across the 14 watersheds (app. 1, table 1.3). Years with more total precipitation and, a related measure, more days with total precipitation greater than or equal to about 0.4 in. (hereafter referred to as “heavy rainfall days”), were associated with higher annual TN concentrations in most watersheds. Thirteen of the fourteen watersheds had a positive correlation between TN and heavy rainfall days (R10D, median r=0.54), and only CAPT HICK had a negative relation (fig. 39). Years with more heavy rainfall produced more storm responses and stormwater runoff, hydrologic conditions that can produce elevated TN concentrations (fig. 35) and can result in the lowest retention of nitrogen inputs on the landscape (fig. 25). These patterns, which are similar to patterns reported from other small urban watersheds (Kaushal and others, 2008a), may reflect increased delivery of stored nitrogen to streams from soils and groundwater or the influence of atmospheric deposition, a nitrogen input that was higher in wetter years. Annual changes in atmospheric deposition had a positive relation with annual changes in TN concentration in most watersheds. Correlations were stronger with wet deposition (N WET DEPO, median r=0.50; fig. 39) than total deposition, which includes the dry deposition fraction (N TOT DEPO, median r=0.35). Atmospheric deposition can be an important nitrogen input in urban watersheds, especially during high-flow events that deliver surface runoff to streams (Kaushal and others, 2011): a process previously documented at Fairfax County streams (Hyer and others, 2016).

Years with higher average air temperature had lower median TN concentrations in most watersheds (AIR TEMP, median r=−0.31; app. 1, table 1.3). This association was also in some measures of temperature extremes, where colder temperatures during the coldest night of the year (TNn, median r=−0.45) and more frost days per year (FD, median r=0.47) were correlated with lower median TN concentrations. Thirteen of fourteen watersheds had a positive correlation between the number of frost days and TN, but only CAPT HICK had a negative correlation (fig. 39). The unique negative correlation at CAPT HICK between annual TN and annual climatic predictors suggests that drivers of TN concentration in this watershed may differ from others in the monitoring network. Relations between warmer air temperature and lower TN concentrations are consistent with effects of denitrification (Pilegaard, 2013) and are similar to relations between air temperature and flow-normalized total nitrogen load reported for streams throughout the Chesapeake Bay watershed (Chanat and Yang, 2018).

Explanation of Variability

Based on the information described above, six predictors were included in LME models to describe the spatiotemporal variability of median-annual TN concentrations: (1) median-annual streamflow, (2) number of heavy rainfall days, (3) number of frost days, (4) impervious cover, (5) wet atmospheric nitrogen deposition, and (6) septic-system density. Median-annual streamflow was used to evaluate whether median-annual TN was affected by varying flow conditions because varying flow conditions can mask the role of other factors. Combined with the six parameters of the linear-growth model, all possible model combinations of the six additional predictors were considered. Models with the lowest AIC included a similar combination of terms and performed similarly (app. 2, table 2.1).

Many similarly performing models described the spatiotemporal variability of median-annual TN concentrations at the 14 study watersheds from 2009 through 2018, but interpretations from a single model can help inform watershed management. A 9-parameter model that included fixed-effect terms for the number of heavy rainfall days and septic-system density was the simplest model form among the top-performing models; more complex models did not offer a substantial improvement over this combination of predictors (app. 2, fig. 2.1). With the exception of a fixed effect to represent the linear rate of TN change over time, the slope term included in all models, the 95-percent confidence intervals of all fixed-effect parameters in the selected model did not include zero, suggesting an interpretable positive or negative effect on TN (table 22).

The fixed-effect terms in the selected model explained 58 percent of the spatiotemporal variability in median-annual TN concentration (table 22). Compared to the linear-growth model (table 21), the selected model explained an additional 30 percent of the annual within-watershed variance, an additional 49 percent of the between-watershed intercept variance, and an additional 79 percent of the between-watershed slope variance. Observed and predicted responses plotted against time were used to visually assess the selected model’s performance in each watershed (app. 2, fig. 2.2). The selected model residuals for most stations appear evenly distributed around zero in most years: an improvement compared with predictions from the linear-growth model (fig. 40). The distribution of residuals around zero in the linear growth model is related to extremes in the annual record of heavy rainfall days (fig. 41A).

The selected model estimated two effects of septic-system density on median-annual TN concentrations. In watersheds with above average septic-system densities (1) concentrations were higher than average and (2) concentration increases from 2009 through 2018 were larger. These effects are consistent with the observed relations between septic-system density and median-annual TN concentrations (fig. 37A and fig. 38B). The influence of septic systems on differences in TN concentration between watersheds has been reported in Fairfax County streams (Hyer and others, 2016; Porter and others, 2020b), but their potential effect on concentration changes over time has not been previously explored. The number of septic systems generally increased over time in each study watershed, but annual changes in septic density were not a predictor of median-annual TN concentrations. A few hypotheses can be deduced from these patterns: (1) TN concentrations are responding to the expansion of septic infrastructure that occurred in previous years, (2) TN concentrations are responding to changes in the delivery of groundwater concentrations of septic effluent, and (or) (3) TN concentrations are responding to changes in other landscape or climatic drivers, which are co-located or processed differently in areas with high septic-system density.

The selected model estimated that, on average, median-annual TN concentrations were higher during years with more heavy rainfall days (table 22), consistent with observed patterns in most watersheds (fig. 39). The annual number of heavy rainfall days was the only modeled effect that caused estimated values to deviate from a perfectly straight line (app. 2, fig. 2.2). The number of heavy rainfall days may affect annual TN concentrations in Fairfax County streams for a variety of reasons. First, it should be noted that relations between TN and the number of heavy rainfall days were similar to TN and average-annual precipitation (PRCP; app. 1, table 1.3), so, more broadly, precipitation may be an important driver of TN concentrations. Wetter years may deliver more nitrogen to streams because of (1) the direct runoff of inputs from the landscape, (2) the flushing of inputs that were previously stored in soils or groundwater, and (or) (3) changes in biogeochemical processes that limit removal rates of nitrogen. TN concentrations were typically higher during elevated streamflow (fig. 35), but it is unlikely that the modeled precipitation effect reflects the effect of direct runoff, because the model was built from samples mostly collected during non-storm conditions. An additional driver of TN during wet years is the increased delivery of atmospheric nitrogen. Annual relations between TN and wet atmospheric nitrogen deposition were similar to the annual relations between TN and the number of heavy rainfall days (fig. 39), and many alternative models included an effect for wet-nitrogen deposition (app. 2, table 2.1).

Management-Practice Effects

Stormwater-retrofit and channel-restoration practices to reduce nitrogen were installed in seven study watersheds from 2009 through 2018 and, normalized by watershed area, are expected to lower TN yields by tens to thousands of pounds. These credited reductions do not coincide with decreases in median-annual TN concentrations in most study watersheds (fig. 42). Credited TN reductions were considered as a predictor of median-annual TN concentrations in the previously described LME model to assess their potential effect on observed conditions, relative to other environmental and anthropogenic drivers. Specifically, the model considered the annual effect of cumulative TN reductions on differences in median-annual TN concentrations from 2009 through 2018 in 14 watersheds. A model including the effect of credited management-practice TN reductions did not reduce the AIC of the previously described model (–47.25 compared to –48.89), and the coefficient was not significantly different from zero (table 23).

Patterns of Observed Responses

TP concentrations, a measure of all dissolved and particulate phosphorus forms in the water column, differed among study watersheds, varied seasonally, and were strongly related to streamflow (Porter and others, 2020b). Median-annual TP concentrations were typically below about 0.03 mg/L in most watersheds, and concentrations as much as 0.06 mg/L were commonly observed in FLAT and Frog Branch above Flatlick Branch at Chantilly, Va. (FROG) (fig. 43). Dissolved phosphorus, predominately in the bioavailable form of orthophosphate, represented most of the non-storm TP concentration in most of the watersheds. Higher orthophosphate concentrations were observed in some streams during warmer months, likely because of temperature-dependent processes that affect the solubility of phosphorus minerals and the desorption of orthophosphate from sediments (Duan and others, 2012). Despite a similar seasonal pattern and range of WT observed at all stations, seasonal differences in orthophosphate were largest in FROG, FLAT, and HPEN (fig. 44). TP concentrations at all stations were highest during high flows because the increased delivery of sediment carries more phosphorus (fig. 45).

Phosphorus concentration trends were evaluated using monthly water-quality samples collected from April 2008 through March 2018 to describe how conditions changed over time in 14 watersheds with 10 years of data, analyses and interpretations originally conducted by Porter and others (2020b; table 24). After adjusting for streamflow, TP concentrations increased in 4, declined in 1, and had no trend in 9 watersheds. A network-wide increasing TP trend of 0.008 mg/L, an increase of 3.9 percent over the 10-year period, was determined by a regional Seasonal-Kendall analysis. The largest TP concentration increases of about 0.02 mg/L occurred in INDIAN and PSB. The only significant decrease occurred in FROG, a decline of about 0.04 mg/L. That decline is notable because this watershed typically had the highest median annual TP concentrations in the monitoring network (fig. 43). At stations other than FLAT, changes in TP were consistent with changes in orthophosphate (table 24). In FLAT, increases in TP at the end of the study period were likely caused by increases in particulate material and coincided with the completion of a stream restoration near the monitoring station (Porter and others, 2020b; fig. 46).

Between-watershed differences in median-annual TP concentrations were typically larger than year-to-year changes within watersheds (fig. 43). These two sources of variability, representing spatial differences among study watersheds and temporal differences from 2009 through 2018, were described by a linear-growth model (table 25). The TP linear-growth model, and all subsequent TP LME models, represent TP concentrations multiplied by a factor of 100, which allows fixed- and random-effect estimates to be interpreted using fewer decimal places. The TP linear-growth model indicated that the average TP concentration across all watersheds in 2009 was 0.025 mg/L and that the linear rate-of-change could not be differentiated from zero. Random-effect variances describe a dataset where group-level variance (3.700), representing between-watershed differences in intercept, accounted for most of the total model variance (4.059). Variances representing TP change over time, differences in between-watershed slope (0.004) and within-watershed annual changes (0.355), were small. The median-annual TP and orthophosphate changes that occurred from 2009 through 2018 included large declines in FROG and smaller increases in most other watersheds (fig. 47). Changes in orthophosphate were consistent with changes in TP in all watersheds except FLAT, suggesting that changes in particulate phosphorus may have affected TP trends in this watershed. The AIC of the linear-growth model (335.42) was compared against models that included additional predictors; more parsimonious models reduced the between- and within-watershed random-effect variances.

Relations to Predictor Variables

Between-watershed TP differences were likely affected by a combination of human and natural factors; patterns that were explored by comparing average predictor values and TP concentrations from 2014 through 2018 at 20 study watersheds (app. 1, table 1.4). The percentage of a watershed within the Triassic Lowlands (TRIASSIC, r=0.75) had the strongest correlation with TP: a pattern previously observed and reported (Jastram, 2014; Porter and others, 2020b). Aside from differences in sediment transport, likely to be a minor influence on these median-annual TP concentrations, TP concentrations may be higher in Triassic Lowland streams for one or a combination of reasons: (1) they contain more nonpoint-phosphorus inputs, (2) they contain more naturally occurring phosphorus minerals that are contributed to streams as dissolved phosphorus through dissolution pathways, and (or) (3) their geologic, ecological, or water-quality properties result in greater rates of dissolved phosphorus mobilization to streams. Some evidence for this last point can be drawn from the previous discussion of nutrient input retention rates in the five intensive-monitoring stations, where FLAT had a slightly lower percentage of nitrogen and phosphorus inputs retained during base-flow conditions than SFLIL, DEAD, or DIFF (table 17).

Phosphorus concentrations may be higher in Triassic Lowland watersheds because of contributions from phosphorus-bearing minerals to streams. Phosphorus is naturally occurring in many minerals, most commonly in apatite and feldspars, and although most of this material is insoluble and not delivered to streams, the weathering of phosphorus-bearing minerals can contribute to elevated phosphorus concentrations in streams (Denver and others, 2010). Fairfax County streams in the Triassic Lowlands have higher naturally occurring dissolved solid concentrations, likely contributed by soluble carbonate minerals (Larson, 1978); however, the phosphorus content of these minerals is uncertain. Not all weathered phosphorus reaches streams because inorganic dissolved phosphorus can quickly attach to surfaces of other solids (“adsorption”), be converted to organic material (“immobilization”), or form complexes with metal compounds (“precipitation”; Denver and others, 2010). These complex biogeochemical interactions are reversible and variable across space and time based on local soil and environmental conditions.

Estimates of naturally occurring phosphorus contributions to streams can be made from data that describe concentrations of bed-sediment samples collected from small, undeveloped watersheds (Nardi, 2014). Phosphorus concentrations estimated from these data are highest in the Triassic Lowland study watersheds (FLAT, FROG, HPEN, and Sugarland Run Tributary below Crayton Road near Herndon, Va. [SGRLND]) and some Piedmont watersheds and lowest in the Coastal Plain (fig. 48A). An alternative representation of mineral phosphorus contributions is available from data that describe phosphorus concentrations in A and C soil horizons (Terziotti, 2019). Soil phosphorus concentrations estimated by these data do not differ by geologic terrane (fig. 48B). Estimates of bed-sediment and soil phosphorus concentrations were not strongly correlated with between-watershed differences in median-annual TP concentrations (app. 1, table 1.4; P GEOL, r=0.39; P SOIL A, r=0.24; and P SOIL C, r=−0.08).

Historical water-quality samples offered additional insights about the contribution of phosphorus from natural or anthropogenic sources. Water-quality samples were collected from 49 Fairfax County streams during a period of low flow in August 1977 (Larson, 1978). Fairfax County’s population and number of housing units has more than doubled since the late 1970s (fig. 5; Han and others, 2018), so these water-quality samples may reflect lower urban impacts. Orthophosphate concentrations in the historical data were lower than concentrations measured in August of the current study (fig. 49). Four of the five orthophosphate concentrations collected in 1977 from Triassic Lowland watersheds were below 0.02 mg/L, and one sample had a concentration of 0.08 mg/L, a value similar to conditions in FLAT, FROG, and HPEN during the study period. No stations in the current study were sampled in 1977, but a concentration below the 0.01 mg/L detection limit was measured in 1977 from a station on Flatlick Branch (USGS station ID 01656900), about 400 ft downstream from the current FLAT monitoring station. Compared with higher concentrations in FLAT during the study period, this historical concentration suggests that orthophosphate concentrations may have increased in FLAT because of human-derived inputs in addition to possible contributions from natural-mineral dissolution.

The shallow soils and fractured bedrock common in Triassic Lowland watersheds may offer fewer opportunities for the retention of phosphorus in the soil and allow nonpoint-phosphorus inputs to reach streams quickly and efficiently. Compared to Piedmont watersheds with thicker soil layers, Triassic Lowland watersheds have shallow fractured bedrock that is poorly buffered from surface inputs (Froelich and Zenone, 1985b). In the Triassic Lowland study watersheds (FLAT, FROG, HPEN, and SGRLND), the average soil depth was about 45 in., compared with about 60 in. in other watersheds (table 26). Soil depth had a strong negative correlation with between-watershed, median-annual TP concentrations (app. 1, table 1.4; SOIL DEPTH, r=−0.70; fig. 50A).

Some additional properties that describe physical soil structure and soil chemical composition were related to between-watershed differences in median-annual TP concentrations (app. 1, table 1.4) and may represent conditions that promote the mobilization of phosphorus to groundwater and streams. Phosphorus concentrations in streams can be high in areas with sandy soils because these areas have low phosphorus adsorption rates (Withers and Jarvie, 2008). The negative correlation of median-annual TP concentrations and soil sand content (SAND, r=−0.71), however, suggests that sand content, which is lowest in Triassic Lowland watersheds (table 26), may not be a driver of between-watershed TP differences. Watershed soils containing aluminum and iron can increase phosphorus retention, because dissolved phosphorus binds and forms insoluble complexes with these ions (Borggaard and others, 2004; Denver and others, 2010). Stream concentrations of aluminum and iron were not measured during the study period, but historical samples collected in Fairfax County between 1954 and 1979 (Froelich and Zenone, 1985b) indicate that some Triassic Lowland watersheds have lower dissolved iron groundwater concentrations than Piedmont or Coastal Plain watersheds (fig. 51). Watersheds in the Triassic Lowlands contained more soil organic matter than other study watersheds (table 26), a property that can increase dissolved phosphorus availability (Deb and Datta, 1967; Nziguheba and others, 1998; Yang and others, 2019; Ye and others, 2006). Despite this influence, differences in soil organic matter percentage were not correlated with between-watershed differences in median-annual TP concentrations (OM, r=0.26).

The capacity of soils to retain water may affect phosphorus retention in multiple ways. Fluctuations between wet and dry conditions can increase rates of phosphorus leaching from soils (Kinsman-Costello and others, 2016; Young and Ross, 2001); lower rates of leaching are typically associated with soils that have slower water transport times (Djodjic and others, 2004). Soils that remain saturated and form reducing conditions, however, can increase the leaching of phosphorus bound to iron (Shenker and others, 2005; Willett, 1989). Properties that describe the storage and retention of soil water were correlated with between-watershed differences in average-annual TP concentrations (app. 1, table 1.4). TP concentrations were typically higher in watersheds with higher soil field capacities (FC, r=0.70), lower rates of vertical hydraulic conductivity (KV, r=−0.74), and lower amounts of plant-available water (AWC, r=−0.67). Triassic Lowland watersheds generally had the highest FC and the lowest AWC and KV values of all study watersheds.

Biogeochemical transformations of phosphorus between soluble and insoluble forms affect phosphorus delivery to streams and are influenced by WT, DO, stream acidity, and microbial communities. Median-annual TP concentrations were higher in streams with higher WT at the time of sample collection (WT, r=0.52). This relation is consistent with studies demonstrating that temperature affects the release of phosphorus from sediment and that biotic processes enhance phosphorus availability (Duan and others, 2012; Duan and Kaushal, 2013; Withers and Jarvie, 2008). Median TP concentrations were higher in streams with lower median DO concentrations at the time of sample collection (DO, r=–0.61; fig. 50B). This relation may be related to the previously described relation between TP and WT because warmer water contains less DO. However, differences in oxidizing or reducing conditions can also affect phosphorus availability, with lower DO concentrations favoring phosphorus release from sediments (Duan and others, 2012; Withers and Jarvie, 2008). The range of DO concentrations commonly observed in Fairfax County streams reflects oxic conditions (fig. 44C), but steep declines in DO can occur within some shallow streambeds (House, 2003), resulting in lower DO concentrations in groundwater than surface water (Denver and others, 2010). Even temporary anoxic conditions that develop near the sediment surface can increase the release of phosphorus from sediments (James and Barko, 2004). Phosphorus availability can also change through pH-dependent processes, but the correlation between median-annual TP concentrations and pH values was not significant (pH, r=0.18).

The phosphorus input of fertilizer was positively correlated with between-watershed differences in median-annual TP concentrations (P FERT, r=0.65). Properly maintained turfgrass can efficiently retain fertilizer inputs, but fertilizer leaching and delivery to streams can occur in response to over applications or applications that precede storm events (Bachman and others, 2016). Phosphorus loading to streams can also increase after years of fertilizer applications that result in phosphorus-saturated soils, which describe soils without available phosphorus-bonding sites (Soldat and Petrovic, 2008). Phosphorus inputs from turfgrass fertilizer may have declined in recent years because of legislative restrictions on the sale of most residential fertilizers containing phosphorus (Harper, 2011), but it can take years to reverse the effects of phosphorus-saturated soils (Kleinman and others, 2011; Sharpley and others, 2013).

Turfgrass area also was positively correlated to median-annual TP concentrations (TURF, r=0.65; fig. 50C) and had the same correlation as fertilizer inputs, because amounts of applied fertilizer were derived directly from the turfgrass cover. Phosphorus inputs may be underrepresented in FLAT if total fertilizer applications are higher on golf courses than residential lawns, because a 200-acre golf course is within the watershed about 2 mi upstream from the monitoring station. Nutrient export from golf courses, however, can be highly variable and, like other managed turfgrass systems, is dependent on proper application and irrigation practices (Bock and Easton, 2020). In addition to representing phosphorus-fertilizer applications, the correlation between median-annual TP concentrations and turfgrass cover may result from other nonpoint inputs of phosphorus associated with residential activities. Sanitary-sewer densities were higher in watersheds with more turfgrass and were correlated with between-watershed differences in median-annual TP concentration (SAN SEWER LEN, r=0.53). Sanitary-sewer leaks in urban areas can deliver nutrients to streams (Kaushal and others, 2011), but their contributions are difficult to quantify (Amick and Burgess, 2000; Fenz and others, 2005) and may be less common in Fairfax County, which has a very low rate of reported sanitary-sewer overflows and backups (Fairfax County, 2021).

Some physical watershed properties were correlated with the 2009 through 2018 percentage change of median-annual TP concentrations (app. 1, table 1.5), which may describe the effects of phosphorus being stored and released from soils. Larger TP increases occurred in watersheds with lower soil cation-exchange capacities (CEC, r=−0.67; fig. 52A), which measures ability of a soil to hold positively charged ions. Watersheds in the Triassic Lowlands had the highest cation-exchange capacities (table 26) because these soils contained higher amounts of organic material, silts, and clays that contribute to negatively charged soil structures. Therefore, larger TP increases occurred in watersheds with more sandy soils (SAND, r=0.59) and less organic matter (OM, r=−0.60). The positive relation between TP percentage change and average soil depth (SOIL DEPTH, r=0.55; fig. 52B) would have been much stronger but for the TP increase that occurred in FLAT.

Annual climatic conditions were correlated with changes in median TP concentrations from 2009 through 2018 (app. 1, table 1.6; and fig. 53). Median TP concentrations were higher in warmer years (AIR TEMP, median r=0.38) and in years with less precipitation (PRCP, median r=−0.33) in most watersheds. The relation with air temperature is consistent with the previously described effects of temperature-dependent phosphorus releases, driven by abiotic and biotic processes. The inverse relation between TP and annual precipitation may be driven by phosphorus geochemical transformations. Although years with more precipitation produce many storm events that can deliver high phosphorus concentrations to streams through sediment runoff and erosion, the phosphorus concentrations considered here are reflective of base-flow conditions and are not likely affected by such processes. Under these longer-term low-flow conditions, differences in TP concentration may be affected by differences in precipitation, as streambed sediments wet and dry and cycle between anoxic and oxic conditions. Although the net effect of these processes on phosphorus mobilization are difficult to predict (Kinsman-Costello and others, 2016), the release of phosphorus from organic matter through microbial mineralization can increase when anoxic sediments are dried (Baldwin and Mitchell, 2000; Dieter and others, 2015). Conversely, dry years and oxic conditions may be associated with the increased storage of phosphorus in bed sediments and lower stream concentrations of phosphorus (House, 2003).

Annual changes in some measures of urban land use were correlated with changes in median-annual TP concentrations from 2009 through 2018 (app. 1, table 1.6). Increases in the amount of impervious cover were associated with higher TP concentrations (IMP, median r=0.37; fig. 53) in most watersheds. Increased impervious cover may account for increased applications of phosphorus inputs, more efficient delivery of phosphorus inputs to streams, or effects of salinization on phosphorus mobilization. Stream salinity, measured as SC, was higher in the more impervious study watersheds, and large pulses of SC were observed during some winter storm events, which was consistent with the effects of road-salt application (Porter and others, 2020b). SC increased in the study watersheds over the past 10 years (Porter and others, 2020b); these increases can mobilize phosphorus (along with other cations and nutrients) to urban streams through various ion-exchange and biogeochemical pathways (Duan and Kaushal, 2015; Haq and others, 2018). Sodium, a common component of road salt, can preferentially displace some trivalent and divalent cations, including aluminum, iron, and calcium (Amrhein and others, 1992; Shanley, 1994), which commonly bind with phosphorus to form insoluble compounds. The stability of these phosphate complexes decreases with increasing salinity and can result in higher dissolved phosphorus concentrations in streams (Duan and Kaushal, 2015; Haq and others, 2018). Watersheds with low soil cation-exchange capacities, where TP concentration increases were largest (fig. 52A), may be more susceptible to these ion-exchange processes. The correlations between median-annual SC and median-annual TP concentrations were inconsistent across the study watersheds (app. 1, table 1.6), suggesting that potential salinization effects on phosphorus mobilization are locally variable, poorly represented by median-annual conditions, offset by other predictors, affected by salinity in previous years, or not strongly influential.

Explanation of Variability

Based on the information described above, seven predictors were included in LME models to describe the spatiotemporal variability of median-annual TP concentrations: (1) total-annual precipitation, (2) average-annual air temperature, (3) DO, (4) soil depth, (5) bed-sediment phosphorus concentrations, (6) impervious cover, and (7) turfgrass cover. Combined with the six parameters of the linear-growth model, all possible model combinations of the seven additional predictors were considered. Models with the lowest AIC included a similar combination of terms and performed similarly (app. 2, table 2.2). This model-building process was repeated after excluding FLAT because the phosphorus increases in this watershed were unique; the structure of top-performing models was similar to those including FLAT.

Many similarly performing models described the spatiotemporal variability of median-annual TP concentrations in the 14 study watersheds from 2009 through 2018, but interpretations from a single model can help inform watershed management. An 11-parameter model was selected for interpretation after considering the results of models with and without FLAT and the most likely factors that may be influencing TP concentrations. This model included fixed-effect terms for total-annual precipitation, soil depth, stream DO concentration, and turfgrass cover (table 27). Alternative models with the same or fewer parameters had similar AIC values, but the inclusion of more parameters did not offer meaningful improvements in model performance (app. 2, fig. 2.3). The 95-percent confidence interval of all model coefficients, except the effect of soil depth on the linear rate of TP change over time (SOIL DEPTH, b₂), did not include zero. The effect of soil depth on the linear rate of TP change over time was consistently identified in models that did not include FLAT and was therefore included in the selected model of 14 watersheds. Additional predictors not included during initial model development were iteratively considered but did not improve the selected model’s performance. These predictors included effects of median-annual streamflow, cation-exchange capacity, organic matter, available water capacity, WT, SC, and the number of dry stormwater ponds.

The fixed-effect terms in the selected model explained 85 percent of the spatiotemporal variability in median-annual TP concentrations (table 27). Compared to the linear-growth model (table 25), the selected model explained an additional 96 percent of the between-watershed intercept variance, an additional 8 percent of the annual within-watershed variance, and an additional 29 percent of the between-watershed slope variance. Observed and predicted responses plotted against time were used to visually assess the selected model’s performance in each watershed (app. 2, fig. 2.4). Notably, the model did not accurately estimate the increase in TP concentrations in FLAT that occurred towards the end of the study period. The selected model residuals for most stations were evenly distributed around zero in most years: an improvement compared with predictions from the linear-growth model (fig. 54).

The fixed effects in the selected model describe the average effect of average-annual precipitation, soil depth, turfgrass cover, and stream DO concentration on median-annual TP concentrations in 14 watersheds (table 27): predictors that vary over space and time (fig. 55). The largest model effect sizes corresponded to effects of soil depth; watersheds with deeper soils had lower TP concentrations (SOIL DEPTH b₁, −0.77) and a more positive TP increase from 2009 through 2018 (SOIL DEPTH b₂, 0.31). Watersheds with higher average DO concentrations and lower average turfgrass cover were associated with lower TP concentrations. The effect of annual precipitation on TP changes within watersheds over time (PRCP b₃, −0.23) was similar in magnitude to that of the model slope (Year b₃, 0.23), highlighting that, on average, TP concentrations were higher in years with less precipitation and in later years of the study period. The positive slope coefficient suggests that TP increases may be affected by additional factors not represented in the model. Annual precipitation was the only modeled effect that caused estimated values to deviate from a perfectly straight line (app. 2, fig. 2.4). The interpretation of these terms was similar to a model of 13 study watersheds that excluded FLAT (app. 2, table 2.3).

The modeled effects of soil depth, turfgrass cover, DO concentration, and annual precipitation may be better understood by considering the retention of phosphorus on the landscape and phosphorus inputs. The modeled effect of turfgrass cover may suggest that phosphorus concentrations are higher in watersheds that receive more fertilizer applications. The modeled effects of soil depth describe lower average phosphorus concentrations that increased more positively over time in watersheds with deeper soils. Because deeper soils may store more phosphorus, changes in TP over time in Fairfax County streams may be affected by the retention or release of soil-bound phosphorus. The storage and release of phosphorus in a watershed is also controlled by environmental conditions such as temperature, DO, acidity, wetness, and microbial communities (House, 2003; Withers and Jarvie, 2008). The selected model identified two such effects: that TP concentrations were lower in watersheds with higher stream DO concentrations and that years with less precipitation were associated with higher TP concentrations. Lower DO concentrations can produce anoxic conditions that favor phosphorus release from sediments (Duan and others, 2012; Withers and Jarvie, 2008). The association between annual changes in precipitation and TP concentrations may be explained by the mineralization of organic matter, which is expected to occur more readily in drier conditions, releasing dissolved phosphorus from soils (Baldwin and Mitchell, 2000; Dieter and others, 2015). However, the net effect of precipitation on phosphorus mobilization is difficult to predict (Kinsman-Costello and others, 2016) and these effects are likely overwhelmed by the delivery of particulate phosphorus to streams during high-flow events. The positive slope term (Year b₃) suggests that the modeled landscape and climatic factors do not completely account for an average increase in TP concentrations from 2009 through 2018.

Management-Practice Effects

Stormwater-retrofit and channel-restoration practices to reduce phosphorus were installed in seven study watersheds from 2009 through 2018 and, normalized by watershed area, were expected to lower TP yields by tens to hundreds of pounds. Despite the credited reductions, median annual TP concentrations typically did not decrease in these watersheds (fig. 56). Credited TP reductions were considered as a predictor of median-annual TP concentrations in the previously described LME model to assess their potential effect on observed conditions, relative to other environmental and anthropogenic drivers. Specifically, the model considered the annual effect of cumulative TP reductions on differences in median-annual TP concentrations from 2009 through 2018 in 14 watersheds. The AIC of a model including TP reductions credited to management practices was nearly 10 units lower than the previously described model (table 28). The modeled effect of management practices was positive, suggesting that credited reductions increase stream phosphorus concentrations; however, this relation was heavily affected by patterns in FLAT. Among the seven watersheds with credited phosphorus reductions, the modeled relation between phosphorus concentrations and credited reductions can be best seen at FLAT, where increases in credited TP reductions coincided with increases in median TP concentrations from 2016 through 2018 (fig. 56). After removing FLAT and refitting the model to 13 study watersheds, the credited phosphorus reduction term did not improve the modeled representation of TP concentrations (AIC of 247.16 versus 245.21 in a simpler model; app. 2, table 2.3), and the term was not significantly different from zero (app. 2, table 2.4).

The modeled effect of management practices on median-annual TP concentrations in FLAT may have been related to the construction of a stream-restoration practice immediately upstream from the monitoring station (fig. 46). Stream restorations that stabilize streambanks may offer some of the best opportunities to reduce phosphorus (Lammers and Bledsoe, 2017), but their construction may cause a short-term increase in particulate water-quality conditions. Construction sites are one of the largest contributors of sediment to surface waters in urban areas (U.S. Environmental Protection Agency, 2005); however, proper sediment and erosion controls can limit the impact of construction on downstream water quality (Houser and Pruess, 2009). Minimal water-quality effects have been associated with construction projects in larger river systems (Colangelo and Jones, 2005; Jastram and others, 2015), and water-quality impacts from the construction of other Fairfax County stream restorations (table 13) were not apparent in the monitoring data, including a stream restoration completed farther upstream from the FLAT monitoring station in 2016. These patterns may suggest that particulate material mobilized during construction is stored in-channel or in floodplains before being delivered downstream. Additionally, changes in flow-adjusted TP loads at FLAT did not show an increase during the 2016–18 period (fig. 31), suggesting that potential construction impacts were not influential during high streamflow conditions, which generate the most TP load.

Patterns of Observed Responses

SS concentrations differed among study watersheds, varied seasonally, and were strongly related to streamflow (Porter and others, 2020b). Median-annual SS concentrations below about 10 mg/L were commonly observed (fig. 57), though concentrations quickly increased by multiple orders of magnitude during stormflow events (fig. 58). SS concentrations and their relation to streamflow are influenced by a combination of factors that affect how quickly runoff enters streams, how much sediment is delivered from streambanks or upland sources, and how in-channel sediment is stored and remobilized. In general, most high streamflow SS concentrations in Fairfax County consist of fine materials, such as silts and clays (Porter and others, 2020b) that originate from stream banks (Cashman and others, 2018). Summer months were typically associated with more frequent and higher-intensity storm events (fig. 15), which result in greater SS concentrations and loads (Porter and others, 2020b).

SS concentration trends were evaluated using monthly water-quality samples collected from April 2008 through March 2018 to describe how conditions changed over time in 14 watersheds with 10 years of data, analyses, and interpretations originally conducted by Porter and others (2020b; table 29). After adjusting for streamflow, SS concentrations increased in FLAT and INDIAN, decreased in BRR, and had no trend in the remaining 11 watersheds. No network-wide trend in SS concentrations was detected. The increases in FLAT occurred from 2016 through 2018 and were consistent with a period of stream-restoration construction immediately upstream from the monitoring station (fig. 46).

Within any watershed and year, a large range of SS concentrations was observed because the monthly samples were collected during a range of streamflow conditions; however, median-annual values removed most of this variability (fig. 57). Between- and within-watershed differences in median-annual SS concentrations, representing spatial differences among study watersheds and temporal differences from 2009 through 2018, were described by a linear-growth model (table 30). This model indicated that the average SS concentration across all watersheds in 2009 was 3.67 mg/L and that the linear rate-of-change could not be differentiated from zero. Random-effect variances describe a dataset where within-watershed annual changes (3.898) represented most of the total model variance (4.437). Group-level variance (0.534), representing between-watershed differences in intercept and differences in between-watershed slope (0.005), was small. Interannual differences in median-annual SS concentration were larger than the linear rate-of-change over time in most watersheds (fig. 59). The AIC of the linear-growth model (615.02) was compared against models that included additional predictors; more parsimonious models reduced the between- and within-watershed random-effect variances.

Relations to Predictor Variables

Between-watershed SS differences were likely affected by a combination of human and natural factors, patterns that were explored by comparing average predictor values and SS concentrations at 20 study watersheds from 2014 through 2018 (app. 1, table 1.7). Most predictors had weak relations or no relations with between-watershed SS differences, likely because the spatial variability of median-annual SS concentrations was small (fig. 57). Average TB conditions during the time of water-quality sample collection had the strongest correlation with SS (TURB, r=0.80; fig. 60A): an expected and well-documented pattern (Jastram and others, 2009; Porter and others, 2020b).

Between-watershed differences in median-annual SS concentrations were positively correlated with stream density, reflecting the length of perennial and non-perennial streams normalized by watershed area (TOTSTREAM, r=0.48; fig. 60B). Watersheds with a higher stream density may have more sinuous stream channels, more opportunities for stream-bank erosion and bed-sediment remobilization, and a faster pathway for particulate material in stormwater runoff to be delivered downstream. The length of streams in small Piedmont watersheds has previously been identified as an important predictor of SS loads throughout the Chesapeake Bay watershed (Brakebill and others, 2010). Greater rates of stream-bank erosion and floodplain deposition have been associated with more sinuous stream channels in Fairfax County streams (Gellis and others, 2017; Schenk and others, 2013).

Sediment sources vary over time and space based on local conditions. Stream banks can represent an important sediment source in mid-Atlantic urban watersheds and are typically a dominant source of sediment in first- and second-order streams (Devereux and others, 2010; Donovan and others, 2015; Noe and others, 2020). In Difficult Run, Fairfax County’s largest watershed (fig. 1), bank-derived sediment is the dominant source during storm events and accounts for most of the SS load (Cashman and others, 2018; Gellis and others, 2017). Much of this bank-derived material represents legacy sediment, that is, sediment that accumulated in the floodplains and channel beds during the colonial period in the 17th to 19th centuries (Noe and others, 2020). Legacy sediment accounts for most of the current bank height in Difficult Run and likely represents a substantial fraction of eroded material transported by Fairfax County streams (Hupp and others, 2013). Bank-derived sediment can be remobilized from in-channel storage, but most material is contributed from the erosion of stream banks (Cashman and others, 2018).

Changing stream-channel shapes and the lateral migration of stream banks can be an important erosional process in urban watersheds (Donovan and others, 2015) and occurs throughout Difficult Run, particularly in response to large storm events (Gellis and others, 2017; Hupp and others, 2013). Urbanization and expanded impervious cover can increase the delivery of stormwater runoff to streams, hydrologic changes that can cause stream-channel widening, migration, and incision (Paul and Meyer, 2001). The historical streamflow record and patterns of sediment erosion and deposition in Difficult Run indicate that such stream-channel changes are likely occurring, and that the watershed is not in geomorphic equilibrium (Hupp and others, 2013). Rates of stream-bank erosion or aggradation, the deposition of sediment along stream banks, vary throughout Difficult Run (Gellis and others, 2017). Difficult Run headwaters are typically erosional and downstream reaches are typically net depositional (Hupp and others, 2013). Patterns in Difficult Run may be indicative of the geomorphic changes occurring in other Fairfax County streams.

Changes in stream geometry were not directly measured in this study, but stream-channel dynamics can be partly represented from the maintenance of real-time streamflow records computed from the relation of water level to measured streamflow (Messinger and Burgholzer, 2018). As the stream channel changes shape, so too does the relation between water level and streamflow, which requires the development of a new streamflow rating curve (Nolan and others, 2005; Turnipseed and Sauer, 2010). From 2008 to 2019, 3 to 6 new rating curves were established at the intensive-monitoring stations. The most changes occurred in DEAD and the fewest in SFLIL, corresponding with differences in watershed impervious cover (table 7) and storm-event hydrologic responses (table 15). Over this same period, only one new rating curve was developed for a nearby headwater stream that drains mostly forested land in Prince William Forest Park and is considered a reference watershed for Fairfax County ecological assessments: South Fork Quantico Creek near Independent Hill, Va. (fig. 1; USGS station ID: 01658500).

Between-watershed differences in median-annual SS concentration were correlated with watershed turfgrass cover (TURF, r=−0.49; fig. 60C). Turfgrass cover may affect how much sediment reaches streams from upland sources by slowing stormwater runoff velocities and increasing infiltration rates, resulting in less sediment loss to streams (Gross and others, 1990; Krenitsky and others, 1998). These benefits can increase as turfgrass density increases (Easton and Petrovic, 2004).

Correlations between watershed predictors and the 2009 through 2018 percentage change in median-annual SS concentrations (app. 1, table 1.8) were influenced by the large SS increase in FLAT. After removing FLAT and reassessing these relations, stronger correlations were identified with stream-bank height (r=0.49), watershed slope (r=0.54), and wet stormwater detention ponds (r=−0.57; fig. 61). Correlations with stream-bank height and watershed slope suggest that median-annual SS concentration increases were larger over the study period in watersheds with taller stream banks and steeper slopes. Both metrics may affect sediment supply and storage: taller stream banks provide more opportunities for erosion and steeper slopes affect the delivery of upland material to streams. Conversely, watersheds with a greater area of wet stormwater detention ponds had a more negative SS rate-of-change over the 10-year period. Wet stormwater detention ponds are permanent pools of water designed to reduce nutrient and sediment loads in downstream receiving waters through chemical processes, the physical settling of particulate material, and attenuating stormflows. Stormwater detention ponds can reduce sediment concentrations and loads in many urban and suburban settings (Hossain and others, 2005; Krajewski and others, 2017; Marsalek and others, 2002), including Fairfax County (Schwartz and others, 2017), though effects can be highly variable.

Annual differences in median SS concentration were correlated with measures of stream TB (TB, median r=0.75; app. 1, table 1.9; fig. 62), an expected relation because SS increases light scatter and reduces water clarity. Stream turbidity measured at the time of sample collection was correlated with DO; lower DO concentrations were typical during periods of elevated TB. Therefore, DO was not likely a factor influencing sediment delivery to streams, despite a correlation between annual changes in DO and SS concentrations (DO, median r=−0.44).

Median-annual SS concentrations were higher during years with colder minimum air temperatures in most watersheds (TNn, median r=−0.38; app. 1, table 1.9; fig. 62), which typically occurred in January or February (fig. 63). The pattern of colder years having higher SS concentrations may be related to freeze-thaw processes, which have been hypothesized as a driver of elevated spring-time concentrations in Fairfax County streams (Porter and others, 2020b). Previous studies have documented high stream-bank erosion rates during the winter and that frost is one of the strongest predictors of stream-bank erosion (Lawler, 1986; Wolman, 1959; Wynn and others, 2008). Soil water expands when frozen, trapping sediment within ice crystals and weakening the cohesive structure of a stream bank. When frozen soil water melts, sediment stored within the ice is carried downstream and the loosened stream bank is more susceptible to erosion during a subsequent stormflow event (Lawler, 1986). Colder minimum temperatures have been associated with greater rates of erosion because these conditions likely result in more widespread frost conditions (Stott, 1997). More frequent freeze-thaw conditions may occur in urban watersheds because stormwater runoff is typically warmer than stream temperatures (Paul and Meyer, 2001; Walsh and others, 2005) and, in recent years, as climatic extremes became more prevalent (Inamdar and others, 2018). Stream banks contributed to a higher percentage of SS concentrations during the spring in Difficult Run, possibly as a result of freeze-thaw processes (Cashman and others, 2018); however, the effect of frost on stream-bank erosion can be overshadowed by large storm events (Gellis and others, 2017).

Explanation of Variability

Based on the information described above, seven predictors were included in LME models to describe the spatiotemporal variability of median-annual SS concentrations: (1) median-annual streamflow, (2) minimum-annual air temperature, (3) stream density, (4) turfgrass cover, (5) stream-bank height, (6) watershed slope, and (7) the area of wet stormwater detention ponds. Median-annual streamflow was considered to evaluate whether median-annual SS was affected by varying flow conditions because varying flow conditions can mask the role of other factors. Combined with the six parameters of the linear-growth model, all possible model combinations of the seven additional predictors were considered. Models with the lowest AIC had common terms and performed similarly (app. 2, table 2.5). This model-building process was repeated after excluding FLAT because the SS increases in that watershed were unique; the structure of the top-performing models was similar to those including FLAT.

Many similarly performing models described the spatiotemporal variability of median-annual SS concentrations in the 14 study watersheds from 2009 through 2018, but interpretations from a single model can help inform watershed management. An 8-parameter model was selected for interpretation after assessing model-performance metrics and considering the interpretations of previous studies on sediment sources and transport in Fairfax County (Cashman and others, 2018; Gellis and others, 2017; Hupp and others, 2013; Schenk and others, 2013). This model included fixed-effect terms for minimum-annual air temperature and stream density (table 31). In models that included and excluded FLAT, this 8-parameter model had one of the lowest AIC values, and the inclusion of additional terms did not greatly improve the representation of SS concentration variability (app. 2, fig. 2.5). This model was selected over an alternative 8-parameter model with similar AIC that replaced stream density with turfgrass cover. Although turfgrass can slow stormwater runoff and reduce sediment delivery to streams (Gross and others, 1990; Krenitsky and others, 1998), stream-bank erosion is an important sediment source in Fairfax County (Cashman and others, 2018; Gellis and others, 2017) and may be related to stream density. The effect of streamflow on interannual SS concentration variability was identified in many of the best-performing models (app. 2, table 2.5), but was not significant when combined with minimum-annual air temperature and stream density.

The fixed-effect terms in the selected model explained 19 percent of the spatiotemporal variability in median-annual SS concentrations (table 31). Because between-watershed SS concentrations were similar (fig. 57), most of this predictive ability was assessed against within-watershed SS concentration variability. Compared to the linear-growth model (table 30), the effect of minimum-annual air temperature in the selected model explained an additional 15 percent of the within-watershed variability in median-annual SS concentrations. Observed and predicted responses plotted against time were used to visually assess the selected model’s performance in each watershed (app. 2, fig. 2.6). Notably, the model did not accurately estimate elevated SS concentrations that occurred in some watersheds and years (for example, OCSB in 2014 and PHCT in 2015). These years had higher median concentrations because most monthly samples were collected during stormflow conditions. Also, the model did not track the increase in SS that occurred in FLAT towards the end of the study period. The selected model residuals for most stations were evenly distributed around zero in most years: an improvement compared with predictions from the linear-growth model, although there are some large under-predictions, likely affected by samples collected during stormflow conditions (fig. 64). Patterns of over- or under-prediction in the linear-growth model appear similar to patterns of years when minimum-annual air temperature was high or low (fig. 65A).

The fixed effects in the selected model describe the average effect of minimum-annual air temperature and stream density on median-annual SS concentrations in 14 watersheds (table 31): predictors that varied over space and time (fig. 65). The largest effect size in the model corresponded to stream density: stations with a greater length of perennial and nonperennial streams per watershed area had higher SS concentrations (TOT STREAM, b₁, 0.74, table 31). The negative coefficient of minimum-annual air temperature explains that SS concentrations were higher in years with colder minimum temperatures. On average, a 2.7-degree decrease in minimum air temperature was associated with a 0.7 mg/L median-annual SS concentration increase. Minimum-annual air temperature was the only modeled effect that caused estimated values to deviate from a perfectly straight line (app. 2, fig. 2.6). The interpretation of these terms was similar to a model of 13 study watersheds that excluded FLAT (app. 2, table 2.6).

The modeled effects of stream density and minimum-annual air temperature may be understood by considering how SS concentrations are affected by physical watershed properties and climatic conditions. Watersheds with a greater stream density may have higher SS concentrations because they have more sinuous stream channels with erosive stream banks, more opportunities for bed-sediment remobilization, and (or) a more efficient pathway for particulate material in stormwater runoff to be delivered downstream. Previous research has documented that bank erosion intensity increases with stream sinuosity in Difficult Run (Hupp and others, 2013), that stream banks are important as a source of sediment in Fairfax County streams, and that lateral streambank migration has erosional effects (Cashman and others, 2018; Gellis and others, 2017; Hupp and others, 2013). The modeled effect of higher SS concentrations in colder years may be related to freeze-thaw processes that increase stream-bank erosion rates. Previously, such effects have been hypothesized to explain differences in seasonal SS conditions in Fairfax County streams (Cashman and others, 2018; Porter and others, 2020b).

Management-Practice Effects

Stormwater-retrofit and channel-restoration practices to reduce sediment were installed in seven study watersheds from 2009 through 2018 and, normalized by watershed area, were expected to lower TSS yields by thousands of pounds. These credited reductions did not coincide with decreases in median-annual SS concentrations in most study watersheds (fig. 66). Credited TSS reductions were considered as a predictor of median-annual SS concentrations in the previously described LME model to assess their potential effect on observed conditions, relative to other environmental and anthropogenic drivers. Specifically, the model considered the annual effect of cumulative TSS reductions on differences in median-annual SS concentrations from 2009 through 2018 in 14 watersheds. The AIC of a model including TSS reductions credited to management practices was about 2.5 units lower than the previously described model (table 32). The modeled effect of management practices was positive, suggesting that credited reductions increase stream SS concentrations; however, this relation was heavily affected by patterns in FLAT. Among the seven watersheds with credited sediment reductions, the modeled relation between sediment concentrations and credited reductions can be best seen at FLAT, where increases in credited sediment reductions coincided with increases in median SS concentrations from 2016 through 2018 (fig. 66). After removing FLAT and refitting the model to 13 study watersheds, the credited sediment reduction term did not improve the modeled representation of SS concentrations (AIC of 546.10 versus 544.10 in a simpler model; app. 2, table 2.6) and was not significantly different from zero (app. 2, table 2.7). These models suggest that the positive relation between SS concentrations and credited sediment reductions is caused by responses at FLAT, a finding similar to the previously described relation between TP and management-practice effects, where sediment may have been mobilized during stream restoration construction.

Patterns of Observed Responses

SC differed between study watersheds, varied seasonally, and was strongly related to streamflow (Porter and others, 2020b). Median-annual SC values below about 400 µS/cm commonly were observed in most watersheds (fig. 67). The ionic composition of Fairfax County streams is primarily influenced by the solubility of naturally occurring geologic materials and urban activities that add salts to the environment and (or) affect their delivery to streams (Larson, 1978; Porter and others, 2020b). CASTLE, LIL DIFF, and SFLIL had the lowest amounts of developed land among the study watersheds (table 7). Median-annual SC values in these three watersheds were typically lower than watersheds with more extensive impervious cover or more soluble geologic materials. Low-ionic strength rainwater often dilutes SC during storm events. However, runoff can carry salts from the landscape that cause large SC spikes, most commonly observed during winter months and associated with roadway deicing and (or) anti-icing chemicals (fig. 68).

SC trends were evaluated using monthly water-quality samples collected from April 2008 through March 2018 to describe how conditions changed over time in 14 watersheds with 10 years of data, analyses and interpretations originally conducted by Porter and others (2020b; table 33). After adjusting for streamflow, SC increased in 10 watersheds and throughout the network by 2.3 percent per year, as determined by a regional Seasonal-Kendall analysis. The largest annual SC percentage increase of 3.6 percent occurred in INDIAN, and the smallest increase of 1.1 percent occurred in FLAT. In most watersheds, the largest SC increases typically occurred during spring and fall months: a possible indication that ions are stored in the watershed and delivered to streams after their initial application to the landscape (Porter and others, 2020b).

SC and major cations and anions were analyzed in water-quality samples collected from Fairfax County streams during base-flow conditions in August of 1977 (Larson, 1978) and can be used to assess how surface-water chemistry has changed over a multi-decadal period. These historical samples were compared to base-flow samples collected during summer months (June–August) between 2011 and 2018 at a station within Difficult Run (n=22 samples, USGS station ID 01646000, fig. 1). On average, calcium and bicarbonate were the most prevalent ions in the historical samples, whereas sodium and chloride were the most prevalent ions in the recent samples from Difficult Run (fig. 69). The recent samples from Difficult Run had significantly higher measures of SC, sodium, and chloride than the historical samples (p-value ≤0.05, as determined through a Kruskal-Wallis test; fig. 70). Average-annual SC values measured at Difficult Run from 1979 through 2018 show that salinity has increased steadily over time; average-annual SC values nearly tripled over this multi-decadal period (fig. 71).

Between-watershed differences in median-annual SC were typically larger than year-to-year changes within watersheds (fig. 67). These two sources of variability, representing spatial differences among study watersheds and temporal differences from 2009 through 2018, were described by a fixed-growth model (table 34). The fixed-growth model differs from a linear-growth model because the fixed-growth model excludes a random-effect slope term ($b_{2i}$ in eq. 1). The random-effect slope term was removed because most watersheds had a similar rate of SC increase over the 10-year period (fig. 72), resulting in a singularity condition caused by a random-effect variance component being reduced to nearly zero. The fixed-growth model indicated that the average SC value across all watersheds in 2009 was 263 µS/cm and that all watersheds had an average SC increase of about 6.9 µS/cm per year. The AIC of the fixed-growth model (1,482.4) was compared against models that included additional predictors; more parsimonious models reduced the between- and within-watershed random-effect variances.

Relations to Predictor Variables

Between-watershed SC differences were likely affected by a combination of human and natural factors, patterns that were explored by comparing average predictor and SC values from 2014 through 2018 at 20 study watersheds (app. 1, table 1.10). Median SC was typically higher in watersheds with higher median streamflow, WT, and pH. Watersheds that have higher median streamflow and WT may receive more runoff than watersheds with lower streamflow and WT. Runoff can carry elevated ion concentrations, contributing to higher median SC. The strongest correlation of median SC was with pH (pH, r=0.64; fig. 73A), where watersheds with higher pH typically had higher SC. This relation may be explained by associations between waters with high-ionic strength and increases in alkalinity. Just as the runoff of salts from impervious surfaces raises stream salinity, weathering of the built environment can contribute ions such as calcium and bicarbonate that raise stream alkalinity (Kaushal and others, 2017). As stream alkalinity increases, the concentration of hydrogen ions in solution decreases, and streams become less acidic. Increases in alkalinity and pH have been observed in many urban settings and are associated with the runoff from, and the weathering of, impervious materials (Kaushal and others, 2013; Kaushal and others, 2017). Based on these previous studies, a causal link between SC and pH in Fairfax County streams is likely, but pH is not likely the driver of higher SC conditions.

Many watershed predictors that describe geologic properties were correlated to between-watershed SC differences, including soil depth, cation-exchange capacity, and organic matter (app. 1, table 1.10). A strong correlation with soil depth (SOIL DEPTH, r=−0.56; fig. 73B) indicates higher SC values in watersheds with more shallow soils. This correlation highlights the higher SC values in Triassic Lowland watersheds (FLAT, FROG, HPEN, SGRLND, and BRR; fig. 67). The Triassic Lowlands contain some of the most soluble bedrock materials underlying Fairfax County streams, and these areas have the highest natural concentration of dissolved minerals in stream water, commonly as calcium, magnesium, and bicarbonate (Larson, 1978). Additionally, SC in some Triassic Lowland watersheds may be affected by the susceptibility of groundwater to surface-water contamination caused by shallow, permeable soils (Froelich and Zenone, 1985b).

Many predictors that describe urban-landscape activities were correlated to between-watershed SC differences. In general, watersheds with more impervious land use, more roadways, higher housing densities, and more wastewater and stormwater infrastructure had higher SC values (app. 1, table 1.10). Many of these relations have been previously documented in Fairfax County streams (Jastram and others, 2009; Porter and others, 2020b) and in other urban/suburban watersheds (House and others, 1993; Paul and Meyer, 2001; Walsh and others, 2005). The percentage of watershed classified as medium intensity developed land, which describes areas with a mixture of constructed and vegetated land cover and most commonly includes single-family housing units (Homer and others, 2015), was most strongly correlated with between-watershed SC differences (DEV MED, r=0.73; fig. 73C). SC may be highest in such areas because of the opportunity for numerous salinity sources, including wastewater, fertilizer, and roadway chemicals, to reach streams through groundwater and surface-water pathways.

Watersheds with more impervious cover tended to have a slightly higher density of stormwater detention ponds; features that were positively correlated with between-watershed SC differences (N DRY POND, r=0.62). This correlation suggests that median SC values were higher in watersheds with more stormwater detention ponds, facilities constructed to manage stormwater runoff and improve water-quality conditions. Although these correlations do not demonstrate causality, previous research has found higher year-round sodium and chloride concentrations in streams and groundwater that are draining stormwater detention ponds (Snodgrass and others, 2017), and that some stormwater detention ponds may be ineffective at reducing stream salinity levels (Minnesota Groundwater Association, 2020).

Some physical watershed properties associated with the Triassic Lowlands were correlated with the 2009 through 2018 median-annual SC percentage change, including sand content, cation-exchange capacity, and soil depth (SOIL DEPTH, r=0.62; app. 1, table 1.11; fig. 74A). In general, these correlations reflect the smaller SC increases in FLAT and FROG from 2009 through 2018, relative to other watersheds. Larger SC increases tended to occur in watersheds with less turfgrass and more high intensity developed land (DEV HIGH, r =0.60; fig. 74B), representing areas with near complete impervious cover, such as apartment complexes or commercial/industrial land uses (Homer and others, 2015). Larger SC increases were also associated with watersheds that had a higher density of roadways classified as interstates and ramps (INT ROAD, r=0.68; fig. 74C). These correlations may be attributed to the role of de-icing and anti-icing roadway salts. Interstates can be large road-salt sources (Cooper and others, 2014), possibly because they receive more salt applications than local roads.

All 14 watersheds had a negative correlation between changes in median-annual SC and changes in minimum-annual air temperature; in other words, higher SC values were observed in years with colder minimum temperatures (TNn, median r=−0.44; app. 1, table 1.12; fig. 75). January or February had the coldest air temperatures of each year (fig. 63). This correlation was weakest at CASTLE, one of the least developed watersheds. The correlation between changes in median-annual SC and changes in minimum-annual air temperature may correspond to the amount of anti-icing and (or) de-icing chemicals used in a given year. A variety of winter maintenance strategies are used to prevent or remove ice from impervious surfaces, including the use of chloride-based compounds (Fay and others, 2015). Sodium chloride is a commonly used anti-icing and de-icing agent because it is less expensive than other chemicals and can be applied in granular form (as rock salt or liquid brine); however, its effectiveness decreases at colder temperatures (Vignisdottir and others, 2019). Magnesium chloride and calcium chloride compounds are more effective under extremely cold temperatures, but higher application rates of all salt-based compounds are still required during lower temperatures (Vignisdottir and others, 2019). Additionally, precipitation events in colder years likely produce more freezing conditions that require salt treatments. As the association between freshwater salinity and road salts have been clearly demonstrated in a variety of urban areas (Cooper and others, 2014; Daley and others, 2009; Kaushal and others, 2005; Kaushal and others, 2018), the observed relations between annual temperature and SC may be driven by differences in application rates or the frequency of de-icing and anti-icing salt use.

All 14 watersheds had a positive correlation between changes in median-annual SC and changes in the percentage of impervious cover (IMP, median r=0.61; app. 1, table 1.12; fig. 75). Impervious cover expanded and median SC values increased in all watersheds from 2009 through 2018, and many predictors related to impervious cover had similar patterns (app. 1, table 1.12). Relations between urban expansion and increased stream salinity have been widely documented and are related to a combination of point and nonpoint sources (Corsi and others, 2015; Daley and others, 2009; Kaushal and others, 2018). Although changes in impervious cover were small in most watersheds over the study period, increases in freshwater salinity have been documented in urban areas during periods of stable-land cover (Bird and others, 2018). Water-quality conditions can take years to stabilize after urban growth, therefore, the multiple decades of rapid urbanization in Fairfax County (fig. 5) may still affect present-day stream conditions.

Explanation of Variability

Based on the information described above, six predictors were included in LME models to describe the spatiotemporal variability of median-annual SC: (1) median-annual streamflow, (2) minimum-annual air temperature, (3) soil depth, (4) impervious cover, (5) medium intensity developed land and (6) interstate roadway length. Median-annual streamflow was considered to evaluate whether median-annual SC was affected by varying flow conditions, because varying flow conditions can mask the role of other factors. Combined with the four estimated parameters in the fixed-growth model, all possible model combinations of the additional six predictors were considered. Models with the lowest AIC included a similar combination of terms and performed similarly (app. 2, table 2.8). Lower SC values were associated with years of higher impervious cover in many of the top-performing models, a relation not supported by observed patterns (fig. 75) and likely caused by the strong correlation between impervious cover and year, the slope term included in all models. The model selection process was repeated after removing the fixed effect of year, thereby accounting for potential multicollinearity issues, to evaluate the effects of time-varying watershed predictors on annual SC differences.

Many similarly performing models described the spatiotemporal variability of median-annual SC from 14 watersheds from 2009 through 2018, but interpretations from a single combination of predictors can help inform watershed management. A 7-parameter model was selected for interpretation after assessing model-performance metrics and considering likely factors of SC response identified in other urban studies. This model included fixed-effect terms for minimum-annual air temperature, medium intensity developed land, and soil depth (table 35). Models with more predictors generally did not explain greater amounts of SC variability and often included a within-watershed effect of impervious cover: a term highly correlated with the model slope term (app. 2, fig. 2.7). An alternative model, where the fixed effect of year was replaced with medium intensity developed land, performed similarly well (app. 2, table 2.9).

The fixed-effect terms in the selected model explained 74 percent of the spatiotemporal variability in SC (table 35). Compared to the fixed-growth model (table 34), the fixed effects in the selected model explained an additional 85 percent of the between-watershed intercept variance and an additional 10 percent of the annual within-watershed variance. The additional variance explained by this latter component, representing interannual SC differences, was small because some of the within-watershed SC variance was already explained by the slope term in the fixed-growth model (Year, b₃). Observed and predicted responses plotted against time were used to visually assess the selected model’s performance in each watershed (app. 2, fig. 2.8). The selected model residuals were an improvement over the fixed-growth model, but selected years still had some over or under predictions (fig. 76).

The fixed-effect coefficients in the selected model describe the average effect of minimum-annual air temperature, soil depth, and medium intensity developed land on median-annual SC values in 14 watersheds (table 35): predictors that vary over space and time (fig. 77). The largest effect size in the model corresponded to the effect of medium intensity developed land: watersheds with more developed land had higher SC values (DEVMI b₁, 0.79). This effect was larger than that of soil depth (ROCKH b₁, −0.39), which indicated that watersheds with shallower soils had higher average SC values. The model slope (Year b₃, 0.38) and minimum-annual air temperature (TNn b₃, −0.27) terms described effects influencing interannual SC differences. These terms indicated that, on average, SC increased over time and was higher in years with colder minimum-annual air temperatures. Minimum-annual air temperature was the only modeled effect that causes estimated values to deviate from a perfectly straight line (app. 2, fig. 2.8). An alternative model that replaced the fixed effect of year with medium intensity developed land related increases in SC to increases in developed land (DEVMI b₃, 0.36; app. 2, table 2.9).

The modeled effects of developed land use, soil depth, and air temperature suggest that variability of SC in Fairfax County streams was related to physical watershed setting and human activities. On average, SC was higher in watersheds that had shallower soils and more medium-intensity developed land. The modeled effect of soil depth likely represents (1) the lower capacity of watersheds with shallow soils to trap and retain landscape-applied salts, and (or) (2) the greater amount of geologic materials delivered to streams through natural weathering. Relations between urban infrastructure and SC have previously been documented in Fairfax County streams (Jastram, 2014; Porter and others, 2020b) and in other urban/suburban watersheds (House and others, 1993; Paul and Meyer, 2001; Walsh and others, 2005). The modeled effect of annual air temperature suggests that SC values were higher in years with colder minimum air temperatures. As the effectiveness of de-icing and anti-icing salts decreases with temperature (Vignisdottir and others, 2019) and freezing precipitation events become more common, years with very cold temperatures likely receive more salt applications. The positive slope term (Year b₃) suggests that, on average, SC increased from 2009 through 2018, though these effects may be accounted for by increases in developed land.

Patterns of Observed Responses

IBI scores varied appreciably between watersheds, in terms of average level and rate-of-change from 2008 through 2017 (fig. 78). In 2009, when samples were available from all 14 stations, IBI scores were lowest in BRR, INDIAN, and DEAD and highest in OCSB, LIL DIFF, and CASTLE. Among these last three, IBI scores in CASTLE were appreciably higher than in the other two, and no annual scores in CASTLE were rated lower than “good.” Trends in IBI scores were evaluated to describe how conditions changed over time in 14 watersheds, analyses and interpretations originally conducted by Porter and others (2020b; table 37). Four watersheds (BRR, DIFF, FLAT, and SFLIL) had a significant increase in IBI score, and the 14-station network had an average-annual increase of 1.8 points. These changes, indicating an improvement in benthic-macroinvertebrate community structure, can be observed in the upward trajectory of IBI scores in most study watersheds (fig. 78). Increased IBI scores from 2008 through 2017 primarily resulted from increases in taxon/group abundance (quantified as percentage composition), though increases in some measures of group taxa richness were also observed. The relative abundance of the most numerically dominant taxon or group in the sample, referred to as percentage dominance, is a component IBI metric for Piedmont streams that decreased in most stations, reflecting improved ecological conditions. The strongest upward IBI trend was at FLAT, where scores improved by about six points per year. In most stations, IBI scores in 2017 were at least one qualitative rating higher (for example, “fair” to “good”) than the first year that samples were collected. In general, average-annual rates of change were more variable in watersheds with lower IBI scores in 2009, whereas smaller increases occurred in watersheds with higher IBI scores in 2009 (fig. 79). These relations also suggest potential differences in response based on the amount of impervious cover in a watershed, an observation that will be discussed in more detail in the follow paragraphs.

Variability in IBI scores between watersheds from 2008 through 2017 was quantified by fitting a linear-growth model to annual observations from the 14 watersheds (table 38). The fixed intercept and slope for the pooled data indicated a network average IBI score of 34.2 in 2008 and an increase of 1.72 points per year. The random-effect variances indicated that roughly 80 percent of total IBI score variance was represented by between-watershed differences in intercept, and only 20 percent was represented by between-watershed differences in slope and within-watershed year-to-year changes. The AIC of the linear-growth model (1,056.1), along with reductions in the three unexplained variance components, were the bases for evaluating the influences of watershed predictors on IBI scores.

Relations to Predictor Variables

The influences of human and natural factors on changes in IBI scores from 2008 through 2017 and between-watershed differences in IBI scores were evaluated by iteratively augmenting the linear-growth model with each watershed predictor (app. 1, table 1.13 and table 1.14). In these and all subsequent LME models, IBI scores were related to predictors representing conditions in the year before sampling because benthic-macroinvertebrate samples were collected annually in the early spring. This process was repeated to model conditions in 2010 through 2017 because, after lagging predictor values by 1 year, median-annual water-quality parameters (Q, WT, SC, PH, DO, and TB; table 3) were not available from all stations in 2008 or 2009. Median-annual water-quality parameters offered little ability to explain spatiotemporal IBI score variability because they did not substantially reduce the AIC of a linear-growth model representing this 8-year period. Seven predictors (defined in table 2 and table 3) reduced the AIC of the linear-growth model by nine or more units (fig. 80), a threshold selected based on the distribution of reductions observed across all predictors.

A model with maximum annual 5-day rainfall reduced the AIC of the linear-growth model by about 11 points, the largest such reduction (Rx5D; fig. 80B). A related heavy rainfall indicator, R95P, reduced AIC by about 9 points. Maximum annual 5-day rainfall ranged from between about 3 to 9 in. from 2008 through 2017 and generally declined over time (fig. 81A). Except for PHCT, correlations between Rx5D and IBI scores across the 14 watersheds ranged from about –0.2 to –0.7 (fig. 81B); that is, lower IBI scores were typically observed in years following high Rx5D values, a relation that was strongest at FROG (fig. 81C). In contrast, at PHCT, prior-year Rx5D was a poor predictor of IBI scores, having a correlation coefficient of about 0.3 (fig. 81D). This distribution of correlations did not appear to be influenced by the amount of impervious cover in a watershed. Effects of antecedent flow on stream invertebrate assemblages are well-documented (Biggs and others, 2005; Greenwood and Booker, 2015; Zuellig and Carlisle, 2019). Using a nationwide set of 51 sites sampled repeatedly from 2002 through 2012, Zuellig and Carlisle (2019) found it necessary to consider sampling day-of-year, antecedent streamflow, and variability in antecedent streamflow when interpreting trends in multi-metric indices of stream ecological condition. The percent numerical representation of the single most dominant taxon is one component of the county’s IBI, and increasing dominance acts to decrease IBI (Fairfax County, 2019). It is plausible that bed disturbance associated with extreme flows might favor, for example, dominance of taxa capable of rapid colonization and reproduction, thereby reducing IBI scores. Such disturbances, however, can act in different directions depending on metric and time scale, emphasizing the importance of jointly considering a suite of covariates (Greenwood and Booker, 2015; Zuellig and Carlisle, 2019).

Watershed average vertical hydraulic conductivity (KV), which describes the vertical infiltration rate of water through fully saturated soils, decreased the AIC of the linear-growth model by about 10 units (fig. 80B). In general, higher KV values are in watersheds with more permeable soils that allow water to move up and down the soil profile with relative ease. Study watersheds with higher KV values had, on average, higher IBI scores in 2009 and lower annual rates of change from 2008 through 2017 (fig. 82A). Average KV was typically higher in less-developed watersheds (table 39), indicating (1) that regions of the county with lower soil infiltration rates had preferentially more development, (2) that extreme impervious cover influenced measurement of soil infiltration rates, or (3) that the correlation between KV and development was spurious. Although observed relations and the AIC reduction of a model including KV suggested that this soil property was an important static predictor of IBI scores, it is likely that this predictor is not causative of IBI score differences. The association between KV and IBI scores may simply result from correlations between KV and various measures of development.

Watershed area decreased the AIC of the linear-growth model by about nine units (fig. 80B). Watershed area did not affect between-watershed IBI differences in 2009, but larger watersheds tended to have larger average-annual rates-of-change from 2008 through 2017 (fig. 82B). In contrast to KV, watershed area had little correlation with the degree of urban development (table 39). The positive relation between watershed area and the rate of IBI change over time is consistent with the findings of Schmidt and others (2019), who found, but did not account for, a positive indirect relation between watershed area and IBI in midwestern streams. Walker and others (2021) found a relation between watershed area and IBI in a randomized sample of Maryland streams; however, the sign of this relation differed between two sampling periods separated by 20 years. The importance of watershed area on changing IBI scores may be related to well-established watershed hydraulic geometry relations, dictating differing non-linear rates of increase in channel characteristics such as width, depth, and velocity with increasing contributing area. Lotspeich (2009) reported such relations based on bank-full channel dimensions measured in 17 undeveloped watersheds ranging in area from 0.3 to 111 mi² in the Virginia Piedmont. Those results indicate that bank-full bed shear stress, proportional to the product of depth and water-surface slope, increases with watershed area, likely decreasing bed stability. This suggests that all else equal, the same runoff-producing rainfall depth falling on two watersheds of differing size will result in more stressful benthic habitat conditions at the outlet of the larger watershed. Thus, more steeply improving IBI slopes in larger watersheds between 2008 and 2017 could reflect a differential response to the county-wide decreasing trend in maximum 5-day rainfall (fig. 81A). Additionally, larger watersheds may contain more diverse habitats, offering refuge for a greater variety of benthic macroinvertebrates (Allan and others, 2021). With less intense rainfall, IBI scores may respond quicker in these larger watersheds as benthic macroinvertebrates recolonize the stream from upstream reaches.

Predictors representing the 2008 through 2017 average value of road density (TOT RD, LOC RD) and low intensity development (DEV LI) reduced the AIC of the linear-growth model by about 9 or 10 units (fig. 80B). Watersheds having higher average values of these predictors, which were correlated with one another (table 39), tended to have lower IBI scores in 2009 and, for DEV LI and LOC RD, more rapid improvement over time (fig. 82C and D). The relation between higher measures of development and lower IBI scores has been previously documented (Paul and Meyer, 2001); however, the effect of development on changing IBI scores over time is less well understood. Although speculative, more rapid improvement of IBI in areas with greater road density may also represent a greater response to decreasing intense rainfall between 2008 and 2017.

Explanation of Variability

All possible combinations of LME models using the 7 previously described predictor variables (fig. 80B) were considered, resulting in models with 6 to 18 parameters. The difference in AIC between the highest-ranking model and the 30th-highest-ranking model was only four AIC units (app. 2, table 2.10), a small degree of distinction relative to the reductions observed when the predictors were considered individually (fig. 80B). This similarity in performance likely results from the strong correlations among considered predictor variables (table 39).

A 12-parameter model, which had the lowest AIC of all modeled combinations, was initially selected for interpretation (app. 2, fig. 2.9). As was the case for many of the best-performing models, this model indicated a singularity condition caused by the random-effect slope variance being reduced to nearly zero. Therefore, the random slope component was removed from the selected model, reducing the number of parameters by 2 (slope variance and slope-intercept covariance), resulting in a 10-parameter model (table 40). This model represents eight fixed effects plus random intercept and residual estimates of variance and has an AIC of 1,015.3, a reduction of just over 40 points relative to the linear-growth model (table 38). The model incorporates the effects of KV and AREA on between-watershed differences in intercept and slope, TOT_RD on between-watershed differences in intercept, and Rx5D accounting for year-to-year within-watershed variability; only the latter predictor enables estimated IBI scores to deviate from a linear regression line (app. 2, fig. 2.10). These predictors reduced the between-watershed unexplained variance in IBI score by about 68 percent relative to the linear-growth model and also reduced within-watershed residual variance by 10 percent. This 10-parameter model reduced the year-to-year mean residuals, relative to the linear-growth model, but was generally less effective in reducing the largest individual residuals in any year (fig. 83).

Arithmetic signs, showing standardized model coefficients in the selected model (table 40), are consistent with the observed relations between predictors and IBI scores (fig. 81B and fig. 82). The largest effect size in the model corresponded to the effect of TOT_RD, indicating that watersheds with higher roadway densities had lower average IBI scores (TOT_RD b₁, −0.63). A similar effect was described by vertical hydraulic conductivity. Watersheds with higher KV values, which occurred in less urbanized areas, had higher average IBI scores (KV b₁, 0.58). A modeled effect of this predictor also influenced the IBI rate-of-change, indicating that IBI scores had smaller improvements in watersheds with higher KV (KV b₂, −0.52). This effect has the same magnitude but opposite sign as watershed area; on average, IBI scores had larger improvements in larger watersheds (AREA b₂, 0.52). These effects are essentially independent because of the weak correlation between KV and AREA (table 39). The standardized estimate for Rx5D, −0.27, indicates that an increase in Rx5D of one standard deviation, about 1.7 in. of precipitation, from one year to the next, was associated with a decrease of roughly 2.5 IBI points. All else equal, if the increase in Rx5D between 2 consecutive years was extreme, such as 4 standard deviations (6.7 in.), IBI would be lowered by about 10 units, plausibly sufficient to lower IBI from one 20-unit-wide qualitative stream-condition description to another (fig. 78). The positive coefficient of the model’s slope term (Year b₃, 0.30) highlights an average increase in IBI scores over time that was unaccounted for by the model’s fixed effects.

These interpretations are challenged by the limitations and uncertainties of using IBI scores to summarize ecological changes. Because it is a composite of roughly 10 simpler metrics quantifying abundance, taxa richness, and pollution tolerance of diverse taxonomic or functional groups, IBI may be influenced by different factors in different watersheds in such a way that prediction of its aggregate space-time behavior is possible in only the most general terms. Moreover, the ability of IBI scores to represent county-wide drivers of ecological change is challenged by differing component metrics for streams in the Piedmont and Coastal Plain. More consistent and interpretable models might result by considering only one of these components or even the abundance or richness of a single taxon or group. Additionally, the observed trends in IBI scores may suggest a subtle bias in responsiveness to changing conditions, because improvement in scores over time is more attainable for watersheds with more degraded initial conditions. Few precedents exist for interpreting drivers of ecological change over space and time using data such as those examined here, highlighting the importance of ongoing monitoring and research.

Management-Practice Effects

The model selected for interpretation was augmented, one predictor at a time, by three time-varying predictors representing values of cumulative management-practice credits for TN, TP, and TSS. In no case did the inclusion of a predictor improve model AIC (table 41). These results do not unequivocally rule out the possibility that one or more of these predictors improved instream ecological conditions as measured by IBI. Any improvements, however, were either not well represented by the cumulative credits time series, or improvements are manifesting (or have manifested) along a non-linear or time-lagged trajectory not readily represented in a LME framework.