Indicators of suspended sediment transport dynamics in rivers

Suspended sediment (SS) is an important part of river systems but is challenging to manage due to the high variability in SS transport dynamics (Cohen et al. 2022), caused by a wide range of climate and catchment factors (Vercruysse et al. 2017). SS concentrations and their temporal dynamics affect aquatic ecological communities and the geomorphic evolution of river systems (Bilotta & Brazier 2008; Neachell 2014; Wharton et al. 2017; Poeppl et al. 2019) and influence decision-making for water resources management, such as the control of intakes for water abstraction and hydropower operations to prevent SS inflows into the systems (Tsyplenkov et al. 2021). With decadal trends in SS flux and increased uncertainties with climate change (Moragoda & Cohen 2020), a deeper understanding of SS transport dynamics is needed.

SS transport dynamics studies focus typically on the calculation and analysis of sediment load, which is the mass of sediment transported per unit time (month or year) based on SS concentration and discharge data (Roman et al. 2012; Yilmaz et al. 2018). Sediment load is useful for calculation of total flux of SS, and in conjunction with the hysteresis in SS rating curves, it can be related to catchment hydrological and geomorphic processes (Horowitz 2003; Misset et al. 2019). Studies on sediment load have identified the relative impacts of factors operating at different spatial and temporal scales, such as different vegetation cover and erosion factors, e.g. gullies (Vercruysse et al. 2017), and relationships between the two main drivers of SS transport, discharge and rainfall (Khan et al. 2016; Kemper et al. 2019). However, it remains difficult to identify the spatio-temporal patterns in, and drivers of, shorter time scales of transport dynamics (days to months) due to the paucity of data.

Most SS transport studies use data from a single temporal scale, often with substantial spatial and temporal gaps, which limits our understanding of drivers (Francke et al. 2014; Vercruysse et al. 2017). A recent study has helped to address this challenge through a multi-scale analysis of the variability in the rating curve relationship between SS transport and river discharge, using turbidity as a surrogate (Wang & Steinschneider 2022). The approach allowed the identification of catchment characteristics that affect the stability of turbidity-discharge relationships, such as stream network complexity, perennial snow coverage, and saturation-excess overland flow. The incorporation of transport dynamic indicators into these multi-scale analyses, i.e. magnitude, frequency and timing (MFT) of events, as is done for other areas of hydrology and water quality, has the potential to further specify the contribution of different factors to SS transport dynamics (Mize et al. 2018; Misset et al. 2019).

An MFT approach has the potential to generate a richer representation of temporal variability of SS transport, which can be combined with analytical approaches, like multivariate data mining techniques, to relate to underlying environmental or human factors (Vercruysse et al. 2017). Magnitude indicators, especially those that capture extreme values, like 95th percentile, are useful to link with extreme values of precipitation, especially flows and sediment yield (Stryker et al. 2018; Whitfield & Shook 2020). Indicators of frequency relate to both the number of events per unit of time and the duration of those events, which can help to identify the underlying factors affecting the hydrological responsiveness of the catchment and sediment generation and transport mechanisms. For example, long-duration SSC events can indicate large catchments and large weather-related drivers, such as hurricanes or monsoons, while short events can relate to convective storms that might be frequent (continental climate) or infrequent (arid climate) (Alexandrov et al. 2007; Yao et al. 2022). The timing of high SSC can be influenced by factors such as climate or land use that influence the generation, mobilization and delivery of SS (Rose & Karwan 2021). These underlying factors can be better understood when a hierarchy of temporal scales and frequency is considered. Previous research has used monthly (Girolamo et al. 2015) and event-based scales (Nosrati et al. 2021) to investigate sediment transport dynamics within watersheds (Kemper et al. 2019).

The study aim is to use an MFT approach to determine the spatial patterns and temporal trends in short-term (daily to monthly) SS transport dynamics that exist at a continental scale. The objectives of the study are to (i) develop and evaluate new indicators of SS dynamics based on MFT, (ii) determine how the new indicators are distributed across the continent and identify spatial patterns and heterogeneity, and (iii) evaluate the indicators to provide recommendations for future studies. Spatial patterns in SS indicators could provide insights into catchment and climate controls on SS generation and transport that could inform environmental water management practices.

Filtering and processing were kept to a minimum to retain the greatest number of sites to investigate spatial patterns in SSC. Sites were selected that had greater than 365 days of record and a minimum of 200 days of SSC data. In total 1,425 sites were selected with an average record length of 9 years (min = 1 year, max = 80 years) and SSC data covering an average of 83.6% of days per site and period of record (min = 6.2% and max = 100% of available daily data). Please see the supplementary material for more information on the sites, catchment areas, and length and completion of SSC records (Figure A-1).

Previous research identified limitations in USGS SSC data that can affect annual SS load calculations (Sommerfield 2016), which has led some researchers to use turbidity as a proxy for SSC (Wang & Steinschneider 2022). A strong correlation was observed between SSC and turbidity, with coherence in the timing and duration of SSC events, for a subset of these sites, which were selected for further analysis in a study on catchment and climate drivers (Shin et al. 2023). Those correlations provided confidence in the use of the SSC data to develop and test the indicators, with the large number of SSC sites providing greater coverage to explore potential spatial patterns in SSC transport dynamics.

The methodology comprised several steps: (i) developing MFT indicators that characterise short-term variations in daily SSC; (ii) evaluating covariation in the MFT indicators using multivariate approaches (Principal Component Analysis, PCA) to reduce the number of indicators; (iii) visualising the spatial patterns in MFT indicators; and (iv) analysing their spatial variability using spatial statistics, clustering and temporal trends using trend analysis.

Magnitude indicators were developed to capture intra- and inter-annual changes in SSC (Figure 2). Average annual values were calculated for each indicator in each water year (Oct–Sep). The magnitude of different percentiles (95th, 75th, 50^th and 25th) was calculated for each year of record and averaged across the number of years per record (M95, M75, M50, M25). Rising and falling rates (i.e. positive and negative changes in SSC per day, respectively) were calculated by subtracting the previous day values from current day values, as the time series consisted of daily data, and then averaged across the period of record (rising rate – R.r., falling rate – F.r.). Peaks in the SSC record were identified from the rising/falling rate time series by detecting positive to negative inflection points, with a maximum peak duration set at 7 days (Yu & Disse 2017). The number of peaks per year was calculated and then averaged across the number of years per record. The date of each peak was exported for use as a timing indicator, more information is given in the following. Finally, the standard deviation of monthly minimum and maximum values was calculated (Sd.min and Sd.max, respectively).

Four α values (0.975, 0.9, 0.6, 0.2) were used to capture SSC events lasting a range of durations, from months to days (see the results section for an example). A value of 30 days was used for the parameter ‘n.reflected’, based on initial testing and previous studies (Ladson et al. 2015). While the filtering process worked well overall, it encountered problems when applied to the longest frequencies due to the presence of data gaps in some sites. Consequently, frequency indicator values were not able to be calculated for the 0.975 and 0.9 α values at 29 and 6 of the 1,425 sites, respectively.

The event duration was calculated based on the contact points identified in the ‘baseflow’ filtering, i.e. the dates when the SSC data and filtering curves had the same value, subtracting the end date from the start date. The analysis was conducted for all four α values (D0.975, D0.9, D0.6, D0.2). As α values become smaller, the filtering curves plot closer to the SSC curve, and short-duration events are identified. The number of events per year was counted and then averaged across the number of years per record (count of events – C.of.e.). The analysis was conducted for all four α values (C.of.e0.975, C.of.e0.9, C.of.e0.6, C.of.e0.2). Finally, SSC peaks that occurred during an SSC event were counted, as described above in the magnitude indicators (P.p.e.). The analysis was conducted for all four α values (P.p.e0.975, P.p.e0.9, P.p.e0.6, P.p.e0.2). Duration and counts were recorded at the event end dates and averaged across the number of years.

Two indicators of the timing of high SSC were developed and calculated, relating to magnitude and frequency (Figure 2). First, the month with the highest SSC (i.e. magnitude) was identified by totalling the daily SSC values for each month across the time series of records at 1,425 sites. Second, the month with the highest number of SSC peaks (i.e. frequency) was determined by totalling the number of contact points (i.e. inflection points at local minima where SSC data and the filtering curves have the same value for the date), as described in the frequency section above, for different α values (0.975, 0.9, 0.6, 0.2) for each month across the timeseries of record. Sites having the same monthly values in different seasons were removed in the end to avoid selecting median month value (α = 0.975, 1,008 sites; α = 0.9, 1,085 sites; α = 0.6, 1,213 sites, α 0.2 = 1,248 sites). The results were aggregated and reported based on the season: winter – December, January & February (DJF); spring – March, April, May (MAM); summer – June, July and August (JJA); and autumn – September, October and November (SON).

A PCA was conducted with the SS transport dynamics indicators using ‘factoextra’ package in R as an initial data simplification step (Hegde 2016). Indicators were selected for further analyses based on being associated with different axes to remove indicators with high degrees of covariation. Following removal of sites with gaps in indicators, a total of 1,364 sites were included in the analysis.

A Mann–Kendall trend analysis from R package ‘Kendall’ was used to detect trends in all average annual magnitude and frequency indicators, for each site (Kar & Sarkar 2021). Sites with significant trends (p value < 0.05) were identified and selected for visualization and discussion.

Large-scale spatial patterns in the SS indicators across the study area were visualised with block statistics in Arc GIS 3.1.0 (Yang et al. 2021). Raster files were created for each indicator based on the location of each monitoring station (cell size: 0.087 radian). The coarse resolution raster output had an average of 1.17 sites per cell (min: 1 and max: 7). The outputs of the block statistics were mean and standard deviation of indicators summarised over a larger area (20 × 20 cells, rectangle), resulting in a coarser-resolution raster. Block statistics were not conducted on the sites with statistically significant trends, due to the low point density.

K-means cluster analysis was used to integrate the findings across the indicators to find general patterns that might relate to natural and anthropogenic factors (Plexida et al. 2014; Celestino et al. 2018; Fang et al. 2021). K-means clustering was conducted in ArcGIS Pro 3.1.0 using the multivariate clustering function. Magnitude and frequency indicators were log-transformed (Tomlinson et al. 2018), and timing indicators were left untransformed. The number of clusters was informed by Elbow and Silhouette methods (Masud et al. 2018), with a final decision based on interpretability (Ballabio et al. 2017).

This study analysed timeseries of SSC data from river monitoring stations using newly developed indicators of SS transport dynamics based on MFT. Clear spatial patterns were observed for some indicators and indicator classes, which suggest that these SSC indicators are driven by large-scale factors, such as climate or LULC. Other indicators had more complex, and heterogenous spatial distributions, suggest more localised influences in SS transport dynamics. In this section, we summarise these findings on the spatial patterns in SS transport dynamics indicators and make recommendations for their future use in studies to determine the influence of catchment and climate factors.

Of the SSC dynamics indicators, M95 had the most distinctive spatial pattern. M95 and the standard deviation of M95 were low across most of the US (Figure 4). However, elevated values were observed in dry and mountainous regions, especially the Rocky Mountain, North American desert and western Great Plains regions (Li et al. 2020a). The high SSC values detected near Mount St. Helens are a clear illustration of the influence of volcanic and other tectonic events on SS transport dynamics, which also impact local ecosystems (Uhrich et al. 2021). It highlights the importance of sediment retention systems near possible eruption areas (Hoess & Geist 2021). There was a strong correlation between the magnitude indicators (Figure 3), thus future studies may only require the inclusion of a single magnitude indicator (e.g. M50 or M95), or perhaps a reduced number with slightly different alignment along the PCA axes (e.g. M95 and M25), to investigate climate and catchment controls.

Frequency indicators had continental scale spatial patterns that were distinctly different to magnitude (Figure 6). SSC events were longer in the upper Midwest and lower in the eastern US at all α values (Figure 6(a) and 6(c)). These results could be caused by spatial patterns in climate and precipitation (Akinsanola et al. 2020). The duration of shorter events (i.e. 1- to 2-week duration) had slightly more defined clusters than long-duration events, especially in the upper Midwest, upper Great Plains and west coast (Figure 6(c)). P.p.e. was generally low across the US, and high values were widely distributed. However, there is a suggestion of a cluster of moderate to high values along the Mississippi River, especially in the upper Midwest (Figure 6(e)), which may be due to short-duration intense precipitation in the summer (Alexandrov et al. 2007). Again, the strong covariation between the duration indicators suggests that a single indicator (e.g. D0.6) may be appropriate for further research (Figure 3).

The spatial pattern of the timing of SSC events suggests a strong climatic signal for magnitude but not frequency. The highest SSC values were found to occur almost exclusively in winter (DJF) in the west coast region of the US and in autumn in Puerto Rico (Figure 7(a); purple). These results correlate to periods of high rainfall and availability of erodible materials (Guirguis et al. 2019), for example, the late summer and autumn Atlantic hurricane season for Puerto Rico. The majority of sites across the US had their highest SSC in spring and summer, but there was no clear spatial pattern, suggesting influences of local factors that affect sediment availability (LULC) overlain onto climatic patterns. For example, the sites with summer (JJA, blue) could be caused by higher and more intense rainfall in this period, while the spring sites might be influenced by agricultural management practices (e.g. ploughing and sowing). While some clusters of high frequency SSC events are apparent (e.g. upper Midwest in summer) that may relate to climate (e.g. convective storms), more research is needed to determine the utility of this indicator (Figure 7(b)).

The trend analysis identified largely decreasing trends in magnitude (Figure 8(a) and 8(c)) in accordance with the overall decreasing trend of SSC (Li et al. 2020b), which is believed to be due to SS management, e.g. soil conservation (Paolo & Minella 2020). However, increasing trends were identified for duration and frequency indicators over large areas of the country (Figure 8(d)–8(f)). Thus, while SSC magnitude is decreasing, periods of elevated SSC (relative to the annual average) are lasting longer in rivers. These results align with a hypothesis that improved soil management is decreasing rainfall-induced erosion and transport to river channels, resulting in a less anthropogenically influenced, flashy SS regime. The localised areas of increasing magnitude might be due to changes of land uses (Murphy 2020). However, due to the limitations and uncertainties in the daily SSC dataset (Sommerfield 2016), interpretation of new indicators about trend analysis requires consideration of uncertainty.

The spatial patterns for MFT indicators were integrated in the cluster analysis, which, despite considerable heterogeneity, showed some important large-scale spatial patterns (Figure 9). Sites in the Great Plains, central Rocky Mountain, and North American Desert regions were associated with the highest M50, R.r. and P.p.e 0.6, which may be associated with LULC and availability of erodible materials, perhaps due to agriculture. Sites with a low M50 and low R.r. with the greatest number of SSC events in summer (i.e. a high month of max i.p. in the hydrological year) were found across the Midwest, Appalachian and northern Rocky Mountain regions. Sites with a moderate M50 and R.r., low D0.6, and a timing of highest SSC in autumn (green) were identified across the US, but a pronounced cluster was found in the Appalachian region, northern California and in pockets on the west coast. The spatial patterns are likely due to a mixture of human impacts and LULC, e.g. higher impacts of magnitude in the driest and mountainous regions, longer duration of events from urban areas, and P.p.e. from areas where receives frequent precipitation events (Li et al. 2020b; Taszarek et al. 2020). For future studies we recommend the use of a combination of magnitude (e.g. M50 or M95, R.r.), frequency (e.g. D0.6 and/or P.p.e 0.6) and timing (e.g. Max sum months and Max i.p.), especially when climate is a factor.

New indicators of SS transport dynamics were used to characterise the MFT of SS transport dynamics. Large-scale patterns in magnitude and frequency indicators were identified, which could be caused by climate and land use. Sites with the highest SSC events (M95) were located in the desert and mountainous regions, while sites with the longest duration of SSC events were found in the upper Midwest. Timing indicators showed varying patterns with clear dominant seasons for higher SSC in the western US and Puerto Rico but a high spatial heterogeneity for the season with the most events (Max i.p.). The trend analysis confirmed the decreasing trend in SSC magnitude reported by other studies but provides new evidence of increasing trends for the duration of SSC events and the number of P.p.e. These indicators provide new insights into the spatial variability of temporal dynamics in SSC in rivers and can be used by studies to determine the catchment and climatic controls. We recommend the use of an MFT indicator for future studies, particularly α = 0.6 to capture short-term dynamics, which may be useful coupled with M95 and R.r. for climate change impact studies. A deeper understanding of the drivers of SSC dynamics will support the prediction of climate change impacts on SSC and develop mitigation measures to reduce the ecological and societal impacts of SS transport.

The authors would like to thank the anonymous reviewers for their careful and thoughtful review of earlier versions of the manuscript. The first author acknowledges bursary support from Cranfield University for this research.

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

The authors declare there is no conflict.