Catchment and climatic influences on spatio-temporal variations in suspended sediment transport dynamics in rivers

Climate change is predicted to affect the flux of suspended sediment in river systems. Changes in rainfall, temperature and CO₂ concentrations are expected to alter the forces (i.e. rainfall intensity, river discharge and wind) that drive erosion and transport but also the factors that moderate these forces, such as vegetation growth (Naylor et al. 2017; Lowe et al. 2018; IPCC 2019). These moderating factors operate naturally and via anthropogenic activity, including changes to land use and land cover (LULC) in response to or as part of mitigating actions against climate change, e.g. shifts in crop types or conversion of agricultural land to forest (Murphy 2020). Considerable research has been conducted on the impact of LULC and climate change on annual sediment load (i.e. sediment yield) at catchment scales (Dai et al. 2009; Hung et al. 2014; Camilo & Restrepo 2017; Kemper et al. 2019) and at continental scales (Syvitski 2003; Walling & Fang 2003; Wang et al. 2011; Syvitski et al. 2014; Moragoda & Cohen 2020). However, there has been less focus on how these factors and their interactions affect suspended sediment transport dynamics, i.e. short-term changes in suspended sediment concentration (SSC) at a location over days to months (Mount & Abrahart 2011; Roman et al. 2012; Murphy 2020).

Uncertain future trajectories in SSC dynamics exist because numerous factors affect SS transport. Climate, LULC and human activities all affect the magnitude, frequency and timing of SS events (Vercruysse et al. 2017). Climate change is shifting temperatures and rainfall, but is also altering the frequency and magnitude of extreme weather events (van de Pol et al. 2010). These climate perturbations can influence SS generation and transport directly, through rainfall-induced detachment of soil particles and resuspension of stored sediment on land and in river channels. They can also alter SS flux at catchment scale through indirect effects on vegetation growth, both positive and negative. The combined effect can lead to multiple possible future trajectories in SSC dynamics in rivers with climate change. Some studies suggest that increased vegetation growth triggered by higher temperatures and carbon dioxide concentrations in the air will reduce soil erosion, despite a higher frequency and magnitude of precipitation events and runoff (Tietjen et al. 2017). Other studies conclude that rainfall has less effect on sediment generation than the catchment hydrological regime, controlled by lithology and geological characteristics, which determines the partitioning of rainfall on the land surface (Fortesa et al. 2021). Natural responses to climate change could be compounded or mitigated by human activity, such as LULC changes, engineering infrastructure and erosion control methods implemented across the landscape. However, suspended sediment transport is not linearly responsive to natural processes or anthropogenic activities (Bussi et al. 2016). While rainfall and discharge are often used to predict suspended sediment transport, relationships between these drivers and SSC are variable in space and time due to impacts by different processes (Vercruysse et al. 2017). Thus, to better predict future changes to SSC transport dynamics (i.e. sediment regime), further research is needed on the climate and catchment factors that interact to affect the magnitude, frequency and timing of high magnitude SSC events in rivers.

Suspended sediment transport dynamics are often investigated explicitly with discharge, either in the form of ratings curves or sediment loads (Vercruysse et al. 2017), but an approach that considers the magnitude, frequency and timing of high SSC events is important for sediment management (Reckendorfer 2019; Shin et al. 2023). A magnitude, frequency, timing approach can also be analysed using regression approaches to determine the contribution of climate and catchment factors, as done in previous research to predict long term mean annual SS loads (Roman et al. 2012; Grauso et al. 2021) and total SS loads (Ayadi et al. 2010). A deeper understanding of climate and LULC controls on SSC and their interactions with catchment characteristics is crucial to predicting changes in the magnitude, frequency and timing of high SSC in rivers in the future.

The aim of this study is to determine the influence of climate and catchment characteristics on suspended sediment transport dynamics in rivers. The objectives are to (i) develop statistical and spatial models of indicators of SSC dynamics (magnitude, frequency and timing) at annual and seasonal scales in order to (ii) assess the relative contribution of climate and catchment characteristics. The aims and objectives are achieved through analysis of long time series of SSC data from river gauging stations distributed across the contiguous United States of America (USA). This statistical approach can generate empirical models capable of differentiating the relative contributions of climate (Roman et al. 2012; Eurich et al. 2021) that have the potential to be more transferable than process-based modelling (Girolamo et al. 2015).

SSC data were obtained from the United States Geological Survey (USGS), which has a large database of SS monitoring data, collected over a long time period, and encompassing a range of climatic regions and catchment characteristics (U.S. Geological Survey 2016). The study started with daily SSC (mg/l) data from 1,667 sites and daily discharge (m³/s) from 1,605 sites from USGS distributed across the continent and overseas territories of the US. Text files containing each SSC (Parameter code: 80154) and discharge (Parameter code: 00060) were downloaded from USGS Surface-Water Daily Data for the Nation. More information on the collection and processing of the SSC data can be found in Shin et al. (2023). The spatial extent of the study was reduced to the contiguous U.S., due to limitation of the site characteristics datasets in Geospatial Attributes of Gages for Evaluating Streamflow (GAGES) (e.g. environmental features and anthropogenic influences) in the remaining areas of the USA (Falcone et al. 2010). The number of sites was then reduced further due to SS data availability and the application of site selection criteria.

The daily SSC dataset was used, rather than turbidity, despite concerns raised by previous research regarding its use for calculating annual sediment yield, in which a lack of strong correlation between SSC and turbidity or suspended solids was noted (Sommerfield 2016). However, analyses for this study found generally high R² and, importantly, coherence in the timings of rises and falls between SSC and FNU (Formazin Nephelometric Unit) turbidity datasets (Supplementary material, Figure A1). The observed correlations for these sites provide confidence in the use of the SSC daily data for this study to investigate catchment and climate influences on suspended sediment dynamics.

The study used summary data on climatic and catchment characteristics from USGS for all 120 sites (Table 1). This dataset aligns with the SSC data pre-processing period; climate data was compiled from the 30 year time period (1971–2000) (Daly et al. 2008) and LULC calculated for the year 2001 (GAGES) (Falcone et al. 2010). The site attribute data are broadly grouped into climate, land use, land cover, soils, and river hydrology data. Ecoregion distinguishes areas with internally consistent LULC and climate in consideration of natural vegetation community (Dodds & Whiles 2004; Omernik & Griffith 2014). The research used Ecoregion level-I (7.0, 11.0, 6.0, 10.0, 9.0, 8.0), which for this study included Marine West Coast Forest (MWCF), Mediterranean California (MC), Northwestern Forested Mountains (NFM), North American Deserts (NAD), Great Plains (GP) and Eastern Temperate Forest (ETF).

The research used a multivariate statistical analysis approach to determine the contribution of catchment and climate explanatory factors to spatial and temporal variations of suspended sediment dynamics. The following steps were used: (i) calculation of SSC dynamics indicators, (ii) data exploration and a principal component analysis (PCA) to inform the reduction in the number of SSC dynamics indicators to model, (iii) regression analysis using ordinary least squares regression and the removal of spatial autocorrelation with geographically weighted regression (GWR), and (iv) cluster analysis and statistical analysis of GWR coefficients by ecoregions.

Indicators to describe SSC dynamics in magnitude (e.g. percentiles and rates of changes in mg/l), frequency (durations and counts) were calculated at annual and seasonal timescales (Shin et al. 2023). The indicators were calculated from the daily SSC time series, separately for each year of record, and then averaged annually and seasonally. A total of 80 SS dynamics indicators were initially calculated (Table 2).

Intra-annual variation in SSC dynamics was investigated seasonally so that the potential impacts of seasonally varying land management practices for the different LULC units (%) could be explored: winter (December, January and February; DJF), spring (March, April and May; MAM), summer (June, July and August; JJA), and autumn (September, October and November; SON) (Table 2). All pre-processed time series were divided into seasons, and magnitude and frequency indicators were recalculated. Visual checks of time series were done at each step to verify whether there were impacts from missing months as well as yearly gaps in time series.

A PCA was conducted, using the ‘factoextra’ package in R (R Core Team 2023), to explore the multivariate dataset and reduce the number of response variables in the subsequent statistical models (i.e. the SSC dynamics indicators) (Haddadchi & Hicks 2021). The analyses used all continuous and discrete (e.g. C.of.e. and P.p.e.) response variables to identify SSC indicators that correspond to different dimensions in the PCA plots, aided by K-means clustering (Fang et al. 2012). The redundant variables with low contribution were removed.

OLS linear regression analysis was used to determine the relationship between climatic and catchment characteristics (Predictor variables) and the SSC dynamics indicators (Response variables) for the selected 120 sites, annual and seasonal. Data transformation was conducted to produce more symmetric residuals with linearity and homoscedasticity. A log 10 transformation was conducted on the response variables to achieve normality (Cho & Lee 2018). The OLS linear regression analysis was conducted using a best subset selection approach with the reduced number of response variables identified in the PCA and all predictor variables. In addition to the identified significant predictor variables, the OLS analysis was deemed appropriate to bypass detrending of the input variables spatial prior to exploring the spatial patterns in the GWR. Due to the number of predictor variables (<30) (Nakariyakul & Casasent 2007), the ‘leaps’ package with Branch and Bound method (Fouad et al. 2020; Lumley 2022) was used to develop parsimonious models in R (R Core Team 2023). The subset selection process started with five independent variables for each model, which was one more than the number of categories of independent variables (Table 1). Variables were added or removed to maximise Adjusted R² and minimise Mallow's Cp (Cp, a variant of AIC – Akaike Information Criterion) and the Bayesian information criteria (BIC).

Cross validation of the models (K-fold) was done to determine predictor variables with the least prediction error using the ‘caret’ package in R (R Core Team 2023). A five-fold model was used (Max et al. 2023) as informed by recent research in environmental sciences and earth surface processes (Coker et al. 2021; Su et al. 2021), with the data divided into test (20%) and training (80%) subsets. The number of independent variables in each model was determined based on the least cross validation errors. Then, the variance inflation factor (VIF) was checked for multicollinerarity, removing independent variables with VIF greater than 10 (Arabameri et al. 2020; Wang & Wang 2022), followed by checking interaction terms between independent variables. Finally, predictor variables with statistical significance greater than p-value of 0.05 were removed.

Spatial autocorrelation was calculated using Moran's I for OLS and GWR using the coordinates available for each site and distance of derived bandwidths.

The spatial distribution of GWR coefficients was investigated directly and indirectly using K-means clustering (Kopczewska & Ćwiakowski 2021). Patterns in GWR coefficients were related to continental scale variations in precipitation. Those were summarised by Ecoregion Level-I. To aid visual comparison of results, only the two predictor variables with the first and second highest relevant weights are presented in figures in the results (Johnson 2010). The significant differences in mean values of the GWR coefficients by Ecoregion were evaluated using one-way ANOVA and Tukey's test (Oh et al. 2019).

Finally, given the large number of SSC indicators investigated, K-means clustering was conducted to identify regions of the U.S. with distinct combinations of GWR coefficients (ArcGIS Pro 2.8). All SS dynamics indicators and statistically significant predictor variables were included in the cluster analysis with cluster size determined by the spatial coverage of Ecoregion level-I. Side-by-side boxplots were used to check variability between predictor variables. Maps of standard error (s.e.) of each GWR model were produced to check locally different reliability of the model and possible local multicollinearity due to implicit constraints tied to the local coefficients, i.e. larger deviations in SSC magnitude (Wheeler & Tiefelsdorf 2005).

For the seven SSC dynamics indicators identified in the PCA analysis, a total of 35 OLS regression models (annual and seasonal) were initially developed for regression analyses. The regression models included different predictor variables and interactions for the SSC indicators, but the most common predictor variables were LULC (Urbanperc, Agperc, Forestperc), BFI, clay content and PPT30 (Table 3). The GWR model with fixed bandwidths accounted for spatial autocorrelation in the regression models, thus yielded models with better R² and AIC (Supplementary material, Table A1). Following tests for significant spatial autocorrelation in the residuals of the GWR model, 23 GWR models had successfully removed autocorrelation and were studied further (Table 3).

Table 3

Final statistical models without significant autocorrelation in the residuals (n = 120, p-value > 0.05)

Response Variable .	Season .	Intercept .	OLS .															GWR (Fixed bandwidths) .
			Model .
			PPT30 .	Rfact .	Kfact .	BFI .	Urbanperc .	Agperc .	Forestperc .	Clay .	Sand .	BFI × Clay .	Rfact × Clay .	R.storage .	Catch.A .	Catch.A × BFI .	R² .	R² .	Observed .
																			Moran's I .
M95	Annual	5.80	−1.00			−6.00	−1.10	−1.00		−4.20		2.00					6.30	8.30	−7.20
	DJF	5.44	−5.00			−6.30		−1.20		−4.30		2.00					5.90	6.80	−2.90
	MAM	3.73	−1.00			−1.70	−9.00	−8.00		4.00							5.50	8.20	−7.40
	JJA	2.29	−2.80	1.10					1.00	7.60			−2.00				5.90	8.40	−5.50
	SON	3.71	−1.10			−2.20	−1.20	−1.10		4.60							5.10	8.70	−5.70
M50	DJF	4.25	−8.00			−4.90	−1.10			−4.40		2.00					5.10	7.60	−4.60
	MAM	3.06	−1.30			−1.30	−1.20	−6.00		3.80							6.30	7.90	−3.80
R.r	Annual	3.25	−9.00			−2.30		−1.30		3.50							4.60	7.70	−6.60
	DJF	4.60				−5.60		−1.60	−7.00	−4.20		1.00					5.40	5.90	−1.30
	MAM	3.00	−8.00			−2.10		−1.00		3.90							4.60	7.50	−6.80
	SON	4.96	−7.00			−6.10		−1.20		−4.60		2.00					5.20	7.80	−6.10
D0.9	DJF	1.12				2.00			−1.00	6.00				8.70			2.70	6.90	−7.90
	MAM	1.22				2.00		1.00							−3.90	9.50	2.70	4.20	−3.10
	JJA	1.20				2.00		2.00									1.60	2.70	−3.30
	SON	1.36				3.00			−2.00		−4.00						2.70	5.10	−5.20
D0.6	JJA	1.20	1.00		−8.09				−3.00								2.00	5.80	−6.20
	SON	1.01							−2.00								1.40	7.00	−1.58
C.of.e0.9	Annual	2.81			4.58				3.00	−4.00				−7.40			3.80	7.10	−7.80
	MAM	9.60			6.88	−2.00	2.00		3.00		4.00						3.40	5.60	−6.80
	JJA	1.35			6.15		2.00		2.00								2.10	4.80	−4.20
C.of.e0.6	Annual	1.35	−2.00		1.35			4.00	8.00								1.90	4.50	−2.70
	DJF	7.59							2.00								6.00	4.30	−6.90
	SON	8.84							2.00		−3.00						1.50	4.80	−7.30

Response Variable .	Season .	Intercept .	OLS .															GWR (Fixed bandwidths) .
			Model .
			PPT30 .	Rfact .	Kfact .	BFI .	Urbanperc .	Agperc .	Forestperc .	Clay .	Sand .	BFI × Clay .	Rfact × Clay .	R.storage .	Catch.A .	Catch.A × BFI .	R² .	R² .	Observed .
																			Moran's I .
M95	Annual	5.80	−1.00			−6.00	−1.10	−1.00		−4.20		2.00					6.30	8.30	−7.20
	DJF	5.44	−5.00			−6.30		−1.20		−4.30		2.00					5.90	6.80	−2.90
	MAM	3.73	−1.00			−1.70	−9.00	−8.00		4.00							5.50	8.20	−7.40
	JJA	2.29	−2.80	1.10					1.00	7.60			−2.00				5.90	8.40	−5.50
	SON	3.71	−1.10			−2.20	−1.20	−1.10		4.60							5.10	8.70	−5.70
M50	DJF	4.25	−8.00			−4.90	−1.10			−4.40		2.00					5.10	7.60	−4.60
	MAM	3.06	−1.30			−1.30	−1.20	−6.00		3.80							6.30	7.90	−3.80
R.r	Annual	3.25	−9.00			−2.30		−1.30		3.50							4.60	7.70	−6.60
	DJF	4.60				−5.60		−1.60	−7.00	−4.20		1.00					5.40	5.90	−1.30
	MAM	3.00	−8.00			−2.10		−1.00		3.90							4.60	7.50	−6.80
	SON	4.96	−7.00			−6.10		−1.20		−4.60		2.00					5.20	7.80	−6.10
D0.9	DJF	1.12				2.00			−1.00	6.00				8.70			2.70	6.90	−7.90
	MAM	1.22				2.00		1.00							−3.90	9.50	2.70	4.20	−3.10
	JJA	1.20				2.00		2.00									1.60	2.70	−3.30
	SON	1.36				3.00			−2.00		−4.00						2.70	5.10	−5.20
D0.6	JJA	1.20	1.00		−8.09				−3.00								2.00	5.80	−6.20
	SON	1.01							−2.00								1.40	7.00	−1.58
C.of.e0.9	Annual	2.81			4.58				3.00	−4.00				−7.40			3.80	7.10	−7.80
	MAM	9.60			6.88	−2.00	2.00		3.00		4.00						3.40	5.60	−6.80
	JJA	1.35			6.15		2.00		2.00								2.10	4.80	−4.20
C.of.e0.6	Annual	1.35	−2.00		1.35			4.00	8.00								1.90	4.50	−2.70
	DJF	7.59							2.00								6.00	4.30	−6.90
	SON	8.84							2.00		−3.00						1.50	4.80	−7.30

View Large

M95 had the greatest number of models, with statistically signficant models created for annual and all four seasons (Table 3). For M50, models were generated for winter (DJF) and spring (MAM). For rising rate, models were created for annual and all seasons except summer (JJA). For these magnitude indicators, there were similar predictor variables and directions of effects (positive, +, or negative, –), generally PPT30 (–), BFI (–), Urbanperc (–), Agperc (–), Clay(+/–), and BFI × Clay (+), though Urbanperc was not included in the model for rising rate. There is some suggestion of seasonal variations in the impacts of predictor variables; Forestperc appears in two seasonal models, the direction of the impact of clay content varies considerably between seasons, and BFI × Clay appears in models for winter (DJF).

For the duration of longer events (α = 0.9), statistically significant models were created for all four seasons, but few predictor variables were common amongst them, except BFI (+). For the duration of shorter events (α = 0.6), models were significant for summer (JJA) and autumn (SON) with Forestperc exhibiting a negative effect in both seasons. An annual model and two seasonal models were created for count of events for both longer and short events showing generally positive associations with LULC and, for α = 0.9, a positive association with K-fact.

GWR coefficients varied spatially and seasonally, suggesting differences and shifts in dominant factors affecting sediment erosion and transport. For example, PPT30 had a negative effect on M95 in all ecoregions, except MWCF, at annual resolution (Figure 3). However, the GWR coefficient for PPT30 was positive only in spring (MAM) and autumn (SON) (Figure 3). Similarly, GWR coefficients for PPT30 were higher in some of the ecoregions, e.g. MC and NFM, in the autumn (SON), with some shifting to positive, and were much lower in other ecoregions in summer (JJA), e.g. NAD. In general, clay content had the 2nd highest weighting in the linear regressions for M95. At annual resolution, clay content was negatively associated with M95 in all ecoregions, but there were substantial differences seasonally. GWR coefficients for clay content were strongly positive in spring (MAM) and autumn (SON) (Figure 3), but negative in winter (DJF) (Supplementary material, Figure A10). In spring (MAM) and autumn (SON), clay content has the greatest positive effect on M95 in the NFM, NAD and GP ecoregions. Finally, the predictor variables with the highest weighting differed for winter (DJF). BFI had a negative relationship with M95, with the lowest GWR coefficient value in ETF and highest in NAD. Agperc had very different relationships with M95 by ecoregions, with a low negative GWR coefficient in the western MWCF ecoregion rising to a high positive in the eastern ETF ecoregion.

Urbanperc, BFI, Agperc, clay cover and PPT30 were the catchment variables affecting M50, mainly in winter (DJF) and spring (MAM) (Table 3 and Figure 4 and Supplementary material, Figure A11). PPT30 had a consistently negative association with M50 for winter (DJF) and spring (MAM), though there were significant differences across the ecoregions. The general west–east negative unimodal pattern was observed for both seasons, with a greater range observed in winter (DJF). Clay cover was positively related to M50 but significantly lower in MWCF and ETF than in other ecoregions. BFI was negatively associated with M50 overall (Table 3), especially in the winter (DJF) when GWR coefficients are significantly lower for GP and ETF (Figure 4). However, BFI was positively associated with M50 in MWCF and MC in spring (MAM) (Supplementary material, Figure A11).

The models predicting rising rates were composed of PPT30 and LULC components, e.g., Agperc, Forestperc, clay content and BFI × Clay (Table 3). However, the predictor variables with the highest weighting were PPT30 and BFI (Figure 5). GWR coefficients for PPT30 displayed the negative unimodal distribution across the ecoregions, identified earlier, with the lowest values (Negative) in NAD and highest (Positive) in MWCF and ETF. BFI was generally negatively associated with rising rates with greater differences by ecoregion in winter (DJF) and spring (MAM). Agperc identified solely with negative coefficients without impact from PPT30.

For duration of events, a much wider range of predictor variables were identified as having high weighting across the seasonal models than for the other magnitude indicators, including Clay, Forestperc, Agperc, BFI, Sand, and K-fact (Figure 6). For the duration of events (α = 0.9), Forestperc and soil types (Clay or Sand) had the highest weightings in the linear regression variables in autumn (SON) and winter (DJF), while Agperc and BFI had the highest weightings in spring (MAM) and summer (JJA). Forestperc was generally negatively associated with duration of longer events in autumn (SON), though there were mixed effects (positive and negative) in winter (DJF). Clay cover was generally positively associated with duration, while sand content was negative. Agperc showed a pattern by ecoregion that was opposite to M95 (DJF), with the highest GWR coefficients in MMCF and lowest in GP and ETF. BFI had a positive unimodal pattern across the ecoregions, with the highest positive contribution to the duration of long events in NAD and GP. Duration of events (α = 0.6) were more related with Forestperc and K-factor. Significant differences by ecoregion were only noted for summer (JJA), when Forestperc and K-factor were negatively associated with duration of shorter events, with significantly lower GWR coefficients in MWCF and NFM than ETF, and GP for K-factor.

For the count of events (α = 0.9 and 0.6), Forestperc was identified as a high weighted predictor variable in all models, annual and seasonal, along with K-factor, BFI, Sand and R.storage (Table 3 and Figure 7). It was generally positively associated with count of events with strong spatial variations in some seasons, i.e. lower in the west and higher in the east in the spring (MAM) for longer duration events (α = 0.9). The association became negative for count of shorter events (α = 0.6) in the western MWCF, MC and NFM ecoregions in winter (DJF). K-fact had a strongly declining impact from west to east, with GWR coefficients highest in MWCF and lowest in NAD, GP and ETF, turning negative for the short duration events (α = 0.6). Sand content in the soil was generally negatively associated with the duration of shorter events (α = 0.6, SON), except in NAD and ETF ecoregions. Finally, both BFI and R.storage were negatively associated with the duration of longer SSC events (α = 0.9).

Sites located in the western U.S. (e.g. California) grouped together consistently, where M95 and M50 were more strongly associated with PPT30, duration of event (α = 0.6) was less strongly associated with Forestperc and K-fact, and count of event was more associated with Agperc and K-fact, relative to the coefficient mean. Sites in the Rocky Mountain region grouped together, except for sites in Montana (Figure 8(c)–8(f)). BFI and Clay had a stronger relationship with M95 in this region in winter (DJF), relative to the coefficient mean. Forestperc with Rfact had a stronger relationship with M95 in summer (JJA). Count of event (α = 0.6) was also affected by Forestperc and Agperc in this region. Sites in the Rocky Mountain and the NAD region grouped together for many of the indicators and seasons, especially M95 and the frequency indicators (Figure 8(a)–8(g)). Southern sites group separately in spring, when M50 is associated with Agperc (Figure 8(d)), and again in the autumn when Urbanperc has a greater influence on M95 (Supplementary material, Figure A17). Sites in the Great Plains, Midwest and Mississippi River regions group together for most of the indicators, especially duration of event (α = 0.6) and count of event (α = 0.9), which were closely associated with K-fact. The regions grouped separated for some of the magnitude indicators. For example, clay content was more closely associated with M95 (JJA) in the Midwest than lower Mississippi. Model coefficients are generally closer to the average than for sites in the western US. Sites (M95) on the Mississippi cluster together all year, except in summer (JJA) (Figure 8(b)). Similarly, the Mid Atlantic and Midwest sites group together, except duration and count of events. For example, in the GWR models for count of events (α = 0.6), the Mid Atlantic has lower coefficient values for Agperc and Forestperc, in comparison to sites located further west (Figure 8(g)). The sites in New England and Mid Atlantic always grouped together for SSC indicators without relatively significant influences from catchment characteristics, except M95 (Agperc and BFI × Clay), event duration (α = 0.6, Forestperc) and event count (α = 0.6, PPT30).

In this section, we discuss the results of the OLS regression models and how the geographical variations in model coefficients generated by the GWR models provide insight into the interplay of catchment and climate factors that affect sediment transport processes (i.e. fine sediment generation, transport and storage overland, and transport in river systems) (Table 4). Indicators of SS dynamics were affected by climate and catchment characteristics with additional spatial and seasonal variations that are indicative of further regional controls that were not captured in the predictor variables. These additional controls could relate to climate, physiographic setting, and anthropogenic activities, which influence the annual and seasonal models.

LULC was the most common type of predictor variable affecting the SSC dynamics indicators. In the OLS models, at least one of the three LULC classes was included in each model for every indicator. However, the direction and strength of the association varied between indicators and seasons (Table 3 and Figures 3–8).

The magnitude indicators (M95 and M50) were negatively associated with Urbanperc, but the frequency indicators (C.of.e.0.9) were dominated by positive coefficients (Table 3). Therefore, in general, higher cover of urban development was associated with lower extreme and average SSC, but a higher frequency of events (Figure 8(f)). From the GWR model, though, the negative association between Urbanperc and SSC magnitude is largely confined to the driest ecoregion (Supplementary material, Figures A10 and 11). For the rest of the USA, there is either no association (GP and ETF) or it is positive (MWCF, MC, and NFM). SS is a major contaminant in stormwater runoff from urban areas (Yin & Li 2008; Borris et al. 2016), which tends to be more problematic in winter-spring than summer (Westerlund et al. 2003; Ciupa et al. 2021; Gong et al. 2021). Increased connectivity due to storm water drains, which carry material washed from impervious surfaces in urban areas (Gellis et al. 2020), and channel incision and scouring due to urban storm water could possibly explain the positive association with M95 and M50 in western regions (Zeiger & Hubbart 2016). The increased flashiness of urban environment explains the increasing number of events. However, the reduction in the magnitude can relate with street cleaning and storm water management. Seasonal differences by ecoregion are likely affected by intra-annual variations in rainfall and discharge (Supplementary material, Figures A10, A11, A15).

Agperc was associated negatively with magnitude indicators (M95, M50 and rising rate), but positively with the duration of long SSC events (MAM and JJA; α = 0.9) and the count of short duration events (DJF and SON; α = 0.6) (Tables 3 and 4). However, there are some important variations by ecoregion, including a positive association with duration of events (Figure 6(a), MAM), and seasonally, such as a negative association with M95 in winter (DJF; Figure 3). In addition, an increasing pattern from west to east (i.e. fewer negative values) was observed for rising rate. This pattern was most pronounced in DJF (Figure 5) but was also observed in MAM and SON (Supplementary material, Figure A12). A similar west–east pattern was also identified for M95, but only for DJF, the pattern was opposite in SON (Supplementary material, Figure A10). Similarly, this opposite pattern (decreasing positive values, west–east) was also identified in the duration of event (α = 0.9; Supplementary material, Figure A13) and count of events (α = 0.6; Supplementary material, Figure A16). These results suggest that Agperc affects seasonal variations in the supply of sediment to rivers, which may be related to the timing of agricultural practices or bare soil, which would varying regionally due to the presence of different crops, irrigation, and soil erodibility that affect soil erosion and delivery, often with time lags (Cooper 2009; Darama et al. 2021; Haddadchi & Hicks 2021). The agricultural lands consist of NLCD classes 81 (Pasture/Hay) and 82 (Cultivated crops) (Falcone et al. 2010), which are LCLU types known to contribute suspended sediment to rivers (Tiecher et al. 2017). However, the generally negative GWR coefficients from agricultural lands are not supported by previous research on soil erosion and delivery to channels but could be explained by confounding factors, such as the spatial distribution of LULC (e.g. farms not in the close proximity to river networks) and erosion mitigation measures such as buffer strips and local policies for maintaining soil resources retention (Tiecher et al. 2017; Shi et al. 2021).

Forestperc was positively associated with the count of long SSC events and negatively associated with the duration of events (long and short) (Table 3). Spatially, the positive association with count varies depending on the seasons but was highest for the count of longer duration events in JJA in MWCF and lowest in the driest region (NAD; Figure 7(a)). The GWR clusters suggest that Forestperc in the Rocky Mountains area experience higher pulses of sediment (M95, Figures 8(b) & Supplementary material, Figure A10) and count of events in summer (α = 0.9 in Supplementary material, Figure A15) and autumn (α = 0.6 in Supplementary material, Figure A16). While Forestperc was mostly negatively associated with duration, which could relate to higher runoff in managed forests, some sites in MC and ETG had a positive association for longer duration (α = 0.9; i.e. seasonal) events. These variable results are related to the impacts of forests on hydrological pathways and soil erosion processes on forest rich areas (Riitters et al. 2004; Zimmermann et al. 2012), complicated by potential covariations between forest cover, slope steepness, management practices, and climate.

Soil and hydrogeological properties (K-factor, clay content, and BFI) have significant relationships with SSC indicators. Soil clay content and BFI were identified in models for magnitude (M95, M50, and rising rate) and duration, while K-factor, and to a lesser degree sand content, featured in models for frequency (count of events).

In general, K-factor is positively associated with the number of SSC events (Figures 7, Supplementary material, Figures A15 and 16), thus areas with more erodible soils have a greater frequency of high SSC events. The spatial pattern is very pronounced for shorter duration events, showing a decreasing pattern west to east, with high coefficients in MWCF, MC, and NFM (Figure 7). For longer duration events, the lack of a clear west–east pattern in the annual models is caused by opposing patterns in different seasons; there is an increasing pattern west to east in spring (MAM) but a decreasing one in summer (JJA) (Figure 6 & Supplementary material, Figure A15). These spatial patterns may be due to the timing of high intensity regional rainfall and availability of erodible materials, especially in California and Great Plains (Figure 8(f)). K-factor is associated with soil structure (e.g. particle size distribution and organic matter) and highly variable depending on LULC. The different distribution of K-factor is likely the result of the impacts of both climate and catchment controls (Martínez-Murillo et al. 2020).

Soil clay content had a variable association with SSC dynamics indicator, varying between magnitude indicators and seasonally within an indicator. While it was negatively associated with the annual model for M95 and winter (DJF) model for M50, it was positively associated with these magnitude indicators for other seasons (Table 3). The GWR results suggest that clay content may have a greater influence on SSC in some ecoregions, especially NFM, NAD, and GP (Figures 3 and 4). However, duration of events (α = 0.9) from the clay cover (Figure 6) did not follow the pattern of Figure 1 in DJF, which may be due to seasonal differences in flashness of the sedimentograph. In cold regions, snowmelt runoff would cause soil loss in spring (Aygün & Kinnard 2021), which would be exacerbated where there is a high proportion of precipitation as snow and LCLU with bare soil or grass (Tramblay et al. 2010; Huang et al. 2020), and may contribute to sediment availability and delivery to the channel. In addition, ephemeral rivers are correlated with fine grained channels with silt or clay (Li et al. 2021).

BFI was negatively associated with magnitude indicators (M95, M50, and rising rate) and positively associated with duration of longer SSC events (α = 0.9) (Table 3). Significant variation was observed in the GWR coefficients spatially and seasonally, with stronger negative associations in winter (DJF) for M95 and M50. For duration and count of events, GWR coefficients of BFI by ecoregion broadly mirror the spatial pattern for BFI (Figure 1), being highest in NFM and NAD, with some evidence of seasonal variations (Figures 5 and 6, Supplementary material, Figure A10–A13). These spatial patterns are likely driven by a combination of climate, geology, soil matrix and groundwater flows (Price 2011). Baseflow is supported by aquifers fed from ephemeral rivers (Goodrich et al. 2018). As baseflow is affected by ground water and soil structure, the negative coefficients from magnitude indicators could be due to dilution impacts from precipitation, and the positive effects to the duration of events might be caused more by additional environmental factors that also correlate with BFI and topography (e.g. the Rocky Mountains, NFM) (Thomas et al. 2004; Price 2011).

These additional environmental factors may also underly the interaction terms included in the models. A positive association between BFI × Clay and magnitude indicators was identified in the OLS models (Table 3). The GWR model coefficients displayed a west–east variation that mirrored many of the other variables, i.e. highest along the coasts (MWCF and ETF) for M95, M50 and rising rates (Supplementary material, Figures A10–A12). The K-means clustering identified the northern Rocky Mountains, NAD and Great Plains regions grouped together (Figure 8(a) and 8(c)), and those are areas where BFI × Clay was strongly positively associated with M95 and M50 (Supplementary material, Figures A10 and A11). These interactions are likely driven by covariations in baseflow, topography and climate, which affect sediment generation, delivery overland to channel, and transport within the river network. For example, mountainous areas often have greater forest cover and a colder climate with a greater proportion of precipitation as snow, which affect rainfall partitioning, hydrological flow pathways, and baseflow generation through snowmelt (Santhi 2008; Huang et al. 2021).

GWR coefficients were significantly different by ecoregion (Level-I) and varied seasonally (Figures 3–7), resulting in large-scale spatial patterns in model coefficients (Figure 8), which could be related to climate factors. Many of the GWR model coefficients were higher in the west coast, lower in the central U.S, and higher in the eastern temperate region, similar to the west-to-east pattern in annual precipitation (PPT30) across the contiguous USA (Figure 1). For example, PPT30 is associated positively with M95 in the moist coastal region (MWCF) but switches from positive to negative seasonally in the semi-arid ecoregions and northwest forest (MC and NFM; Figure 3). PPT30 is a positive driver of event frequency in ETF (Supplementary material, Figure A16).

The annual precipitation indices used in this study were generally negatively associated with magnitude indicators of SSC transport dynamics. However, there is evidence of differences by ecoregion that could relate to the timing and intensity of the predictor variable (PPT30). For example, while GWR coefficients were generally lowest for M95 in NAD over most of the seasons, they were significantly lower in GP and ETF in DJF (Supplementary material, Figure A10). In general, the negative relationship between magnitude indicators and PPT30 is likely caused by dilution effects (Kozak et al. 2019), which is moderated regionally by different frequencies and intensities of rainfall that affect the generation and transport of sediment (Hancock 2012). However, PPT30 is a static representation of 30-year rainfall (1971–2000), which may not capture well the seasonality inherent in the SS delivery system. Rfact is positively associated with M95 in JJA (Table 3), with the highest coefficients found in the Rocky Mountains, desert and Northern Great Plains (Figure 8(b)). The Rfact values may represent well the intensity and frequency of rainfall runoff generation cycles at these locations in these months (Supplementary material, Figure A10). Although larger standard errors (s.e.) in GWR may indicate the effects of local collinearity (Fadliana et al. 2019), climate controls reacted to the catchment characteristics more in wetter seasons, e.g. positive impacts of Rfact (M95, Figure 3) and PPT30 (D0.6, Supplementary material, Figure A14) in summer.

This research identified possible climate and catchment controls on SSC dynamics at continental scale using daily SSC and annual average descriptors of potential predictor variables using a regression and spatial analysis approach. Further modelling (e.g. improvement of gap filling techniques for catchment characteristics with finer time intervals, applying AI or machine learning) (Tao et al. 2021) can help to validate these findings to support the prediction of changes in SSC dynamics in the future with climate change. While this study used a part of GAGES dataset because of the site coverage and suitability of data types, future studies using SSC data from US gauging stations could consider using both GAGES and GAGES-II. These datasets have a greater number of potential explanatory variables, including seasonal precipitation data, which would align more closely with seasonal summaries of SSC dynamics indicators.

The research applied recently developed suspended sediment dynamics indicators to investigate how catchment and climate factors influence spatial and seasonal variations in SSC. Regression models were developed for seven SS transport dynamic indicators for 120 SSC monitoring sites on rivers across the US at annual and seasonal timescales. The models included a range of statistically significant predictor variables, of which LULC, PPT30, baseflow index, and soil properties, like clay content, were the most common across the SS transport dynamics indicators. When additional spatial modelling was conducted to remove spatial autocorrelation (GWR), significant differences in coefficients were identified by ecoregion, suggesting regional influences from catchment and climatic factors. Key results are as follows:

• Agricultural percent cover is positively associated with longer duration SSC events and more frequent short events, and seasonally with SSC magnitude, likely due to effects on soil hydrology and sediment supply.

• Forest percent cover is positively associated with shorter and more frequent SSC events, and extreme high magnitude SSC in summer months, which might be related to the impacts of forest management on hydrological processes or covariations with topography and climate.

• Urban percent cover was the least common LULC variable in the models. It is positively associated with the frequency of long SSC events, especially in spring and summer, and with extreme high SSC magnitude in the western US.

• Soil erodibility (K-factor) is positively associated with the frequency of SSC events, but negatively with the duration of short events, especially in the west coast and forested mountain regions. Thus, more erodible soils yield more frequent and short-lasting high SSC events, though climate and catchment variables likely have strong influences.

• Soil clay content is positively associated with average SSC in all ecoregions over most seasons (spring, summer and autumn) and with the duration of long SSC events, especially in winter. Clay content may relate to the availability of fine sediment, but there are spatial covariations with topography and climate that suggest an influence of snowmelt in mountainous regions and ephemeral rivers in dry regions.

• Baseflow index is positively associated with the duration of long SSC events and slow rates of SSC increase, but negatively associated with average and extremely high SSC, which likely relates to reduced overland flow and more stable flow regimes that reduce sediment erosion and transport.

• Annual precipitation is negatively associated with SSC magnitude and rates of SSC increase overall, likely due to dilution effects. However there are important spatial variations with positive associations in moist forested regions (MWCF and ETF).

• The rainfall and runoff factor is positively associated with extreme high SSC in desert, Great Plains and ETF regions in summer, possibly related to intense seasonal storms.

This study is one of the first to investigate the short-term dynamics (days to months) of suspended sediment transport in rivers at continental scale and the underlying catchment and climatic variables that control spatial and seasonal variations. It provides new understanding of the varying influence of LULC on the magnitude and frequency of high suspended sediment concentrations caused by climate and catchment factors. Further research should be conducted on the impact of climate on suspended sediment transport dynamics to inform future monitoring, modelling and mitigation related to the climate change.

The authors would like to thank the NWIS-TS Application Support Team and National Water Quality Network Coordinator from the USGS for their assistance with SS and site attribute data. The authors appreciate greatly the constructive feedback from an anonymous reviewer on early versions of the manuscript, which significantly improved the final article. The first author acknowledges bursary support from Cranfield University for this research.

The original SSC data were obtained freely from USGS water data portal (https://waterdata.usgs.gov/nwis). The processed data used in the modelling are available at CORD data repository (https://dx.doi.org/10.17862/cranfield.rd.22233721).

The authors declare there is no conflict.