A Union of Dynamic Hydrological Modeling and Satellite Remotely-Sensed Data for Spatiotemporal Assessment of Sediment Yields

: (1) The existing frameworks for water quality modeling overlook the connection between multiple dynamic factors affecting spatiotemporal sediment yields (SY). This study aimed to imple-ment satellite remotely sensed data and hydrological modeling to dynamically assess the multiple factors within basin-scale hydrologic models for a realistic spatiotemporal prediction of SY in watersheds. (2) A connective algorithm was developed to incorporate dynamic models of the crop and cover management factor (C-factor) and the soil erodibility factor (K-factor) into the Soil and Water Assessment Tool (SWAT) with the aid of the Python programming language and Geographic Information Systems (GIS). The algorithm predicted the annual SY in each hydrologic response unit (HRU) of similar land cover, soil, and slope characteristics in watersheds between 2002 and 2013. (3) The modeled SY closely matched the observed SY using the connective algorithm with the inclusion of the two dynamic factors of K and C (predicted R 2 (PR 2 ): 0.60–0.70, R 2 : 0.70–0.80, Nash Sutcliffe efﬁciency (NS): 0.65–0.75). The ﬁndings of the study highlight the necessity of excellent spatial and temporal data in real-time hydrological modeling of catchments. P.P. A.A.-H.;


Introduction
Soil erosion is the process by which land surface is washed away by hydrogeological factors, meteorological agents, and human interference. It leads to long-term changes in the environment and results in reduced water quality, which ultimately affects environmental and human health. Soil erosion decreases soil fertility as well as agricultural productivity by land degradation, increased flooding and landslides, and surface and ground contaminant diffusion by the inflow of sediments to rivers [1][2][3][4][5][6][7]. In the United States, numerous studies have reported that soil erosion is a significant environmental hazard owing to its critical ecologic effects [8][9][10]. Various attempts to estimate soil loss due to erosion involve the development of models such as the Universal Soil Loss Equation (USLE) Model [11,12], Revised Universal Soil Loss Equation (RUSLE) Model [13], Modified Universal Soil Loss Equation (MUSLE) Model [14], USLE-M [15], dUSLE [16], Water and Tillage Erosion Model/Sediment Delivery Model [17,18], and the Agricultural Non-Point Source Pollution Model [19]. With increasing human interventions and potential management, as well as climatic change, there is a growing need to monitor water quality in landscapes with noticeable spatial and temporal transformations, especially on and adjacent to watersheds [20][21][22]. USLE is the well-known and most used empirical model for estimating long-term average annual soil loss globally. Supporting this, [23] showed that around 40% of soil erosion and sediment yield (SY) models adopted input parameters from the USLE model for water quality modeling. USLE estimates long-term annual soil loss and guides proper cropping, management, and conservation strategies [24,25]. However, is atmospherically corrected [61]. Such advancements elevated the uncomplicated use of remotely sensed MODIS products. These data were implemented in various kinds of soil erosion and SY models, which help to quantify the concentrations of suspended sediments in watersheds [62].
Recently, more studies have made use of the concept of remote sensing for data and modeling of C-factors and K-factors because of its easy accessibility, low cost, rapid and reliable data analysis, and lower instrumentation. Furthermore, remote sensing data can be easily integrated with Geographic Information System (GIS) and geoprocessing modules. It helps in the evaluation of land covers and the quantification of soil loss and sediments with spatiotemporal heterogeneity for estimation of water use and management [63][64][65]. Therefore, remotely sensed data, geo-computation and visualization, and geographic environmental modeling work cumulatively to determine the spatiotemporal dynamics of cover management and soil erodibility, and these factors interact in a complex fashion to produce a realistic representation of SY. Nevertheless, remote sensing products have been uncertain in the results in other SY studies owing to specific variables that were directly or indirectly connected to the spatiotemporal variation of SY such as land cover, soil moisture content, and topographical characteristics (slope gradient, slope length). Few studies have evaluated the relations of C-factor, K-factor, and SY to land use change [25], urban land covers [41], dense forestation [24], vegetated zones and crop cover management that have not been accounted for in the conventional USLE equations for SY predictions. Hence, the present study used annually varying land use land cover data for the study areas to understand the vitality of dynamic land cover modeling for realistic soil loss predictions. The enhanced vegetation index, which details the vegetative property of a hydrologic unit, was considered as the primary variable in the novel C-factor modeling. Additionally, the variables of soil moisture content [66], slope length [67], and slope gradient [68], which are relevant for soil loss estimation, were used as the independent variables affecting SY for the watersheds of the study.
Currently, there is a lack of connective model frameworks for SY predictions in different river basins with varying spatiotemporal and hydrogeological conditions that can be manageably incorporated into hydrologic models. Furthermore, studies are scarce that merge ready-to-use functionalities from remote sensing into K-factor and C-factor estimation methods. This study is the very first to incorporate a connective algorithm of C-factor and K-factor into SWAT to dynamically predict SY in the watersheds of the southeastern United States, with the aid of remotely sensed MODIS data. This union augments spatially and temporally significant estimations in watersheds when compared to the traditional USLE equations. Hence the study fills the deficiencies of the USLE formulation of K and C factors in the SWAT model and reinforces the potential of using remote sensing to facilitate hydrological modeling for better water use estimates. Ideally, the study provides a useful model enhancement utilizing readily available and unambiguous data that can be used for predictive purposes.

Materials and Methods
This study is a union of two dynamic factors of the USLE equation into SWAT to advance SY predictions using satellite remotely sensed data. The C-factor and K-factor models were developed using dynamic remotely sensed MODIS data of enhanced the vegetation index, the fraction of photosynthetically active radiation, the leaf area index, and soil moisture. A GIS and Python programming language-based connective algorithm was proposed for the incorporation of the C-factor and K-factor models into SWAT. Then sensitivity analysis was performed to evaluate their influence on SY estimation in watersheds. The spatial distributions of the K-factor, C-factor, and SY at annual time scales were generated and analyzed for every hydrologic response unit (HRU) of similar land cover, soil attributes, and elevation characteristics within the watersheds of the study. The impact of the connective framework of the K-factor and C-factor on SY was evaluated for three cases: (1) using the traditional factors of C USLE and K USLE , (2) using the dynamic C-factor and Remote Sens. 2022, 14, 400 4 of 24 the traditional K USLE , and (3) using the dynamic K-factor and the dynamic C-factor. The study thus enhances spatiotemporal predictions of SY in the SWAT model. Furthermore, agreement between the real-world observations and the spatiotemporal predictions from the satellite remotely sensed data-based connective system of K-factor and C-factor will strengthen their reliability and global usability for real-time water quality modeling. It will also help to analyze the interconnectivity of the USLE factors and watershed erosion for better water use management.
2.1. Dynamic Models of C-Factor and K-Factor 2.1.1. C-Factor Model Pedotransfer functionality of the C-factor was developed based on remotely sensed geospatial data including enhanced vegetation index (EVI), fraction of photosynthetically active radiation (SR in %), leaf area index (LAI), soil moisture content (AWC in %), slope gradient (S), and percentage area for every HRU of similar land use, soil, and slope characteristics in the watershed (A). The MODIS EVI data are of 250-m resolution and contain the best possible pixel value over 16 days. The LAI and SR data product is a 500-m resolution product on a sinusoidal grid over an 8-day observation period. The AWC data product represents the soil moisture content (%) in the surface from 0-10 cm, which determines the nutrients such as nitrogen, phosphorus, and organic carbon in soils. The values of A and S were obtained from the model outputs of watershed delineation identified for the corresponding HRUs.
The remotely sensed environmental data, including the EVI, AWC, LAI, and SR, were obtained per pixel from the processed MODIS imageries for the southeastern United States [69,70] (Table 1). The remotely sensed datasets of EVI, AWC, LAI, and SR were generated by the spatial superposition of the remotely sensed data for the spatial level of HRU in the watersheds. The C-factor model used for the vigorous estimation of C USLE in this study is given here [71]: where a and b are parameters that determine the shape of the exponential curve of C and EVI. Values of 1.5 and 1 were assigned for a and b, respectively, and the coefficients presented in Equation (2) were employed, which yielded the best fit for the C-factor model for the study areas (R 2 = 0.68, PR 2 = 0.51, p < 0.05).

K-Factor Model
The dynamic functionality of the K-factor in the study was developed using the topographic factor (LS USLE ) and the crop and cover management factor (dynamic C-factor), as well as the soil properties of moisture content (AWC in %), bulk density (BD in g/cm 3 ), and permeability (Psoil in mm/h). The topographic factor, LS USLE was calculated using the reference equation [72]. The values of the variables, such as slope length L (m) and slope steepness or slope gradient S (m/m), were obtained from the model outputs of the watershed delineation identified for the corresponding HRUs. The slope S was calculated from the DEM of the watershed. The soil properties of BD and P soil were obtained from the soil attribute characterization module (.sol) of the SWAT model. They were calculated by the spatial join of the soil map (soil type) and the HRU map (HRU ID) in the SWAT model. Thus, the developed model of the K-factor serves as a dynamic and realistic improvement of the K USLE equation in terms of capturing the HRU, as well as annual variations in soil erodibility, and is given here [73]: where C-factor is the crop and cover management factor developed using Equation (1).

SWAT Model
This study uses the basin-scale hydrologic model SWAT as a platform for the implementation of dynamic functionalities to enhance soil loss predictions in watersheds. The SWAT is a physically based, distributed parameter model that operates on an interface called ArcSWAT. It is used for long-term analysis of hydrologic components and predicts the transport of sediments and contaminants in the watershed scale with varying soils, land uses, and management conditions [74]. The concept behind the modeling of spatial units in SWAT is the assortment of HRUs, which are the portions of a sub-basin that possess unique land cover, slope gradient, and soil characteristics. SWAT uses USLE to calculate SY [75]. The SY from USLE in each HRU on a given day, denoted by SY USLE (metric ton/ha), are obtained from the SWAT model using the equation where EI USLE is the rainfall erosion index in 0.017 m-metric ton cm/cu·m h, K USLE is the K-factor in 0.013 metric ton cu·m h/cu·m metric ton cm, P USLE is the USLE support practice factor, LS USLE is the USLE topographic factor, and CFRG is the coarse fragment factor.

Development of Connective Algorithm
The study developed a connective algorithm to simulate spatiotemporal sediment yields in SWAT by incorporating two models of C-factor and K-factor. These two equations of K-factor (Equation (4)) and C-factor (Equation (2)) were used for the first time in modeling multiple watersheds and to predict sediment yields. The novel functions of K-factor and C-factor, which were developed using spatiotemporally dynamic satellite remotely sensed data, are expected to tackle the uncertainties in SY predictions arising from the conventionally used static input variables for estimating K-factor and C-factor. Thus, they will improve SY predictions. The proposed algorithm incorporating dynamic models of C-factor and K-factor into SWAT was based on modeling a yearly variation scheme for the traditional equations of C USLE and K USLE in the SWAT model. The algorithm adopts an input-oriented functionality and runs on annual temporal conditions in the level of sub-basins and HRUs. The newly incorporated components include modification of SWAT equations (SWAT-Equation), addition and editing of SWAT input files (SWAT-Edit), and extraction of SWAT output files (SWAT-Extract). SWAT-Equation was used to modify the equations in the routines and subroutines. SWAT-Edit was used to add new model input files, as well as to edit the existing model input files of SWAT such as crop (land management and vegetation characterization module) and sol. Additionally, SWAT-Edit Remote Sens. 2022, 14, 400 6 of 24 edits the model input files of SWAT simultaneously for crop and sol modules. First, the remotely sensed factors that take part in the algorithm were fed into the respective input modules of SWAT and saved using SWAT-Edit. Second, each row pointing to one HRU was modified to the new functions for C-factor and K-factor (Equations (2) and (4)) using SWAT-Equation. Then SWAT was called to simulate SY for the study watersheds for the required temporal conditions. Further, the spatiotemporal maps for the C USLE , K USLE , K-factor, C-factor, and sediment yields were processed. The described processes are shown in Figure 1.
x FOR PEER REVIEW 6 model input files of SWAT simultaneously for crop and sol modules. First, the rem sensed factors that take part in the algorithm were fed into the respective input mo of SWAT and saved using SWAT-Edit. Second, each row pointing to one HRU was ified to the new functions for C-factor and K-factor (Equations (2) and (4)) using SW Equation.
Then SWAT was called to simulate SY for the study watersheds for the req temporal conditions. Further, the spatiotemporal maps for the CUSLE, KUSLE, K-facto factor, and sediment yields were processed. The described processes are shown in F 1.

Sensitivity Analysis
The developed algorithm connecting models of the C-factor and K-factor into S was validated using sensitivity analysis of SY predictions in three different cases. The sitivity analysis involved three cases of SY predictions for annual conditions in the level using spatiotemporal C-factors and K-factors. Case 1 represented SY with the tional values of CUSLE and KUSLE. Case 2 represented SY with dynamic C-factor and tional KUSLE. Case 3 represented SY with the dynamic C-factor and dynamic K-facto spectively. All cases were correlated to observed sediments at the monitoring statio understand the effect of the C-factor and K-factor functions in realizing the spatia temporal dynamics in SY estimation.

Performance Evaluation
In this study, the performance of the developed algorithm in advancing spatio poral SY predictions was evaluated using three statistical methods. The coefficient o termination R 2 typically ranges from 0 to 1 and is expressed as

Sensitivity Analysis
The developed algorithm connecting models of the C-factor and K-factor into SWAT was validated using sensitivity analysis of SY predictions in three different cases. The sensitivity analysis involved three cases of SY predictions for annual conditions in the HRU level using spatiotemporal C-factors and K-factors. Case 1 represented SY with the traditional values of C USLE and K USLE . Case 2 represented SY with dynamic C-factor and traditional K USLE . Case 3 represented SY with the dynamic C-factor and dynamic K-factor, respectively. All cases were correlated to observed sediments at the monitoring stations to understand the effect of the C-factor and K-factor functions in realizing the spatial and temporal dynamics in SY estimation.

Performance Evaluation
In this study, the performance of the developed algorithm in advancing spatiotemporal SY predictions was evaluated using three statistical methods. The coefficient of determination R 2 typically ranges from 0 to 1 and is expressed as and SS tot = ∑ y i − y mean 2 (6) y and f represent the observed data and predicted data, respectively. SS tot and SS res represent the total sum of squares proportional to the variance of the data and the sum of squares of residuals, respectively. Next, the predicted R 2 (PR 2 ) was employed, which indicates the wellness of a model in predicting responses to new observations. One of the most used coefficients of determination in hydrology is the Nash-Sutcliffe model efficiency coefficient (NS).
where X simulated_average is the mean simulated stream flows and SY of the watersheds averaged per HRU.

Spatial Interpolation and Mapping
In the present study, the inverse distance weight (IDW) method was used for computing the spatial patterns of the remotely sensed variables of EVI, AWC, LAI, and SR for all the HRUs in the watersheds annually. The raster data was used to obtain the HRU wise spatial maps of EVI, AWC, LAI, and SR. Each HRU was organized into grids where an HRU contains values representing EVI, AWC, LAI, and SR [76,77]. Later, they were used to estimate and geospatially map the modified C-factors and K-factors.

Temporal Trend Detection
The t-test was employed as a parametric trend detection test to understand the HRU wise annual trends of C-factor, K-factor, and SY [78,79].

Case Study Areas
Two case studies were chosen as representations of the water quality constituent of SY evaluated using dynamic soil erodibility and cover management factors with changing space and time. Two basins, the Tampa Bay watershed (TBW) in Florida and the Winyah Bay watershed (WBW) in South Carolina, were used in this study ( Figure 2). The TBW is located on the Gulf coast of west-central Florida and has an area of approximately 590,522 km 2 . The second study area of WBW is approximately 575,590 km 2 and encompasses the neighboring areas of the tidal waters of the estuary. The TBW highlights the hydrological responses of a coastal zone with rapid urban development and the WBW highlights the impact of coastal development on increasing soil erosion.

Data and Modeling
The elevation data was rasterized as a 30 m resolution digital elevation model (DEM) from the National Elevation Dataset (NED) provided by the United States Geological Survey (USGS). The soil data concerning texture, depth, and drainage attributes were rasterized from vector maps supplied by the Web Soil Survey (WSS) under the United States Department of Agriculture [80]. The study area watersheds were delineated using the DEMs. The annual land cover maps of the watershed HRUs were developed using a supervised classification analysis based on remotely sensed MODIS images. The developed classifier was trained for the study area watersheds, which produced their supervised classification with fifteen land cover classes [73]. The annually classified land cover maps were employed as the land cover unit in the SWAT model for enhancing hydrological predictions in SWAT.
The models were developed in SWAT for both study areas. A watershed delineation was performed, and several sub-basins were obtained. Then the overlay of land cover, soil, and the slope was carried out, producing the HRUs within the sub-basins. The threshold values used for land cover, soil, and slope while defining the HRUs were 7%, 12%, and 12%, respectively. Later, climate data such as precipitation, temperature, solar radia-

Data and Modeling
The elevation data was rasterized as a 30 m resolution digital elevation model (DEM) from the National Elevation Dataset (NED) provided by the United States Geological Survey (USGS). The soil data concerning texture, depth, and drainage attributes were rasterized from vector maps supplied by the Web Soil Survey (WSS) under the United States Department of Agriculture [80]. The study area watersheds were delineated using the DEMs. The annual land cover maps of the watershed HRUs were developed using a supervised classification analysis based on remotely sensed MODIS images. The developed classifier was trained for the study area watersheds, which produced their supervised classification with fifteen land cover classes [73]. The annually classified land cover maps were employed as the land cover unit in the SWAT model for enhancing hydrological predictions in SWAT.
The models were developed in SWAT for both study areas. A watershed delineation was performed, and several sub-basins were obtained. Then the overlay of land cover, soil, and the slope was carried out, producing the HRUs within the sub-basins. The threshold values used for land cover, soil, and slope while defining the HRUs were 7%, 12%, and 12%, respectively. Later, climate data such as precipitation, temperature, solar radiation, relative humidity, and wind speed were incorporated, model setup was completed, and simulation runs were performed. The stream flows and sediment load from each HRU were calculated separately using input data for weather, soil, topography, vegetation, and land management practices. Later, they were merged to determine the total loadings from the sub-basin as well as from the HRUs. More details and descriptions of the water balance, soil erosion, and water quality process equations can be found in the SWAT technical documentation [72]. Twelve years (2002-2013) of daily weather data were collected from ground weather stations as well as from downscaled and projected data from climate.gov [81,82] with a spatial resolution of 250 m. The observed monthly stream flows and sediment concentrations (X observed ) obtained from USGS for the years 2002-2013 at the respective monitoring stations were used to calibrate and validate the SWAT model. All the datasets for the study were collected for the years 2002 to 2013 (Table 1). Simulations were done for land cover, soil type, and climate condition of 2002 to 2013 using the weather data with the existing management practices, which generated the simulated stream flows and SY from 2002 to 2013 (X simulated ). Tables 2 and 3 explains the descriptive statistics of remotely sensed and modeled soil properties for K-factor estimation in TBW and WBW respectively. The soil texture is stable for a relatively long period of five (2009-2013 in WBW) to eight (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) in TBW) years in certain parts of the watersheds. However, the majority of the HRUs showcased annual variations in the soil property of AWC, and some of the HRUs showed annual variations for BD and P soil .

SWAT Calibration and Validation Results
The two case studies were calibrated and validated in varying spatial and temporal conditions. The two different spatial extents (TBW in Florida and WBW in South Carolina) and two different temporal extents (2002-2009 and 2009-2013) were introduced to evaluate the efficacy of SY predictions by the dynamic C and K factors in the SWAT models when there are changes in hydrogeological and spatiotemporal scenarios. The SWAT model was initially set up for the two watersheds in the southeastern United States. TBW was developed with the land cover, soil type, and climate condition for the year 2002 with 29 sub-basins, 2286 HRUs within 29 sub-basins, 15 land cover classes, and 21 soil categories. Similarly, WBW was developed with the land cover, soil type, and climate condition for the year 2009 with 31 sub-basins, 2118 HRUs within 31 sub-basins, 15 land cover classes, and 32 soil categories. The spatial resolution of the HRUs in the study ranged from 9 km to 16 km with an annual temporal resolution. The SWAT model simulations were employed at the levels of seasons, months, and years in the study models. The soil loss predictions were estimated seasonally and monthly in the study watersheds, as opposed to conventional SY predictions on a long-term timespan [83][84][85]. Both the watershed models were calibrated and validated for the years from 2002 to 2013 with the land cover, soil type, and climate conditions for the respective years. The calibration and validation procedure were performed for stream flows as well as SY by selecting the most sensitive parameters affecting stream flows and SY through one-factor-at-a-time sensitivity analysis [86,87]. The most sensitive parameters of ESCO, CN 2 , and AWCS were used to calibrate stream flows, and ERORGN, ERORGP, and HRU_SLP were used in SY calibration (Table 4). The study performed three trials each for stream flow and SY calibration, as only six parameters were significantly sensitive to the hydrologic responses of the two basins [87]. AWCS, ESCO, and CN 2 were the most critical parameters affecting stream flow patterns in wet conditions because they were found to be highly sensitive to 70% of the stream flows. Even though the parameters such as surface runoff lag coefficient (SURLAG) and Manning's coefficient for overland flow (OV_N) area are considered to be important in surface flow dynamics, they were obtuse in this study [88,89]. The parameters such as HRU_SLP, ERORGN, ERORGP, average slope length (SLSUBBSN.hru), sediment concentration in lateral and groundwater flow (LAT_SED.hru), and a linear and exponential parameter representing the sediment re-entrained through channel sediment routing (SPCON.bsn, and SPEXP.bsn, respectively) are the conventionally crucial parameters affecting SY [90,91].
However, all the parameters except HRU_SLP, ERORGN, and ERORGP were found to be less sensitive to SY in the study. The methods of auto-calibration (ERORGN and ERORGP) and manual calibration (HRU_SLP) were collectively employed for fast-paced and manageable iterations within the SWAT model. The summary results of calibration and validation are given in Table 5

Results of Connective Algorithm Validation
The developed K-factor and C-factor and the USLE K-factor and C-factor for the HRUs of the two watersheds were utilized to estimate spatiotemporal SY in the study. The performance results of the statistical evaluation for the implemented model cases are depicted in Table 6. R 2 , PR 2 , and NS ranges between 0.5 and 1 are considered to have acceptable goodness of fit in hydrologic modeling purposes [92]. The R 2 , PR 2 , and NS values of the generated SY during Case 1 were 0.582, 0.531, and 0.603, respectively, for TBW and 0.503, 0.462, and 0.418, respectively, for WBW. In TBW, they increased to 0.719, 0.636, and 0.741, respectively, after the C-factor data were applied. For the WBW, the increases were 0.594, 0.533, and 0.527, respectively. The dynamic K-factor function (Equation (4)) constitutes the C-factor estimates developed using the remotely sensed variables of EVI, LAI, SR, and AWC (Equation (2)). It designated the repeated implementation of C-factor information in Case 3, where dynamic C-factors and K-factors are connected. PR 2 acts as a potent tool to eliminate the overfitting of the model due to the repetition of the C-factor variable [93][94][95]. The increase of PR 2 from Case 2 to Case 3 indicated that the twin use of the C-factor variable does not result in overfitting of the SY in Case 3 (TBW: 0.636 (Case 2) to 0.668 (Case 3); WBW: 0.533 (Case 2) to 0.644 (Case 3)).  Figure 4 shows the observed SY versus simulated SY for Cases 1, 2, and 3 at the Rocky Creek monitoring station of TBW and monitoring stations of the Black River, Lynches River, and Pee Dee River in WBW between 2002 and 2013. The significant deviations in actual and simulated SY in Case 1 were significantly reduced in Case 2 during the low SY conditions. This demonstrates that the dynamic C-factor model with the traditional K USLE model is a better simulator of the low SY that is generally found at the inlet of a watershed and in the dry seasons. Most estimates of lower SY converged to the actual measured values in TBW (NS: 0.603 to 0.741) with the application of the dynamic C-factor and traditional K USLE , whereas they showed a minor impact in WBW (NS: 0.418 to 0.527). Interestingly, both study areas exhibited potential and further nearness of observed and simulated SY with the connectivity of dynamic K-factor and C-factor (NS: 0.603 to 0.744 (TBW); NS: 0.418 to 0.683 (WBW)). This shows that model efficiency of SY predictions increased for every dynamic factor added; for all cases, the highest model efficiency occurred when the two dynamic factors were connected (Case 3: PR 2 : 0.60-0.70, R 2 : 0.70-0.80, NS: 0.65-0.75). These findings emphasized that Case 3 is a much better representation of the actual sediment prediction for different watersheds with varying spatial and temporal environments compared to Case 2. Hence, the connectivity of the dynamic estimates of C-factor and K-factor into SWAT realistically predicted SY when compared to the dynamic assessment using the singular C-factor model [96,97].  The correlations and predictions of Case 3 indicated that the combined application of C-factor and K-factor powerfully represented the SY estimates both for the whole watershed and for the HRUs in the watersheds. This finding is corroborated by some recent studies [98,99]. Despite the different spatial resolution of the various remotely sensed MODIS variables of the study ranging from 30 m to 500 m, the study could demonstrate the advancement of the existing algorithm in predicting the SY. The differences between the actual and estimated SY accounting for the soil erosion for the spatial and temporal extents in the watersheds can be attributed to the modeling techniques, interpolation methods, and real hydrogeological conditions involved in the study areas. These differences also include the ambiguity in the space-wise trends from the IDW technique when interpreting different factors from remote sensing, which were employed in models for catchment scales. Additionally, dynamic modification of the other USLE factors such as rainfall erosion index, support practice factor, and topographic factor might improve the precision and reliability of the spatiotemporal SY estimates in watersheds. The correlations and predictions of Case 3 indicated that the combined application of C-factor and K-factor powerfully represented the SY estimates both for the whole watershed and for the HRUs in the watersheds. This finding is corroborated by some recent studies [98,99]. Despite the different spatial resolution of the various remotely sensed MODIS variables of the study ranging from 30 m to 500 m, the study could demonstrate the advancement of the existing algorithm in predicting the SY. The differences between the actual and estimated SY accounting for the soil erosion for the spatial and temporal extents in the watersheds can be attributed to the modeling techniques, interpolation methods, and real hydrogeological conditions involved in the study areas. These differences also include the ambiguity in the space-wise trends from the IDW technique when interpreting different factors from remote sensing, which were employed in models for catchment scales. Additionally, dynamic modification of the other USLE factors such as rainfall erosion index, support practice factor, and topographic factor might improve the precision and reliability of the spatiotemporal SY estimates in watersheds.

Spatiotemporal Predictions of C-Factor and K-Factor
The annual C-factors and K-factors for each HRU in the watersheds were calculated using the dynamic functions from Equations (2) and (4) to investigate their distribution in the watersheds. The spatial C USLE and K USLE simulated by the SWAT model along with the spatially and temporally dynamic C-factors and K-factors averaged from 2002 to 2009 for TBW and from 2009 to 2013 for WBW are reported in Figure 5. For both watersheds, the C-factor range of 0.01-0.50 was almost 60% of the total watershed, which indicates the dominance of grass and shrubs, pasture, and vegetation in TBW and WBW. The proportion of urban spaces and forests (C-factor: 0-0.01) and wetlands and water (C-factor: 0.5-1) were 35% and 5%, respectively, for TBW and 10% and 30% for WBW. The spatial distribution characteristics of the C USLE deviated from the annual C-factors averaged from 2002 to 2013. C-factors were found to increase in Lower Tampa Bay, indicating the presence of deforestation as well as land management activities in the Lower Tampa Bay region [22]. The distribution of K USLE for Lower and Middle TBW was underestimated in comparison with the developed annual K-factors by 45% and 35%, respectively. In contrast, it was overestimated by around 20% and 28% in the northwest and central portions of the WBW, respectively. In both watersheds, no significant annual dynamics were observed in K-factor between 2002 and 2013, despite the noticeable deviation of annual K-factors from the K USLE .
Remote Sens. 2022, 14, x FOR PEER REVIEW 14 of 24 1) were 35% and 5%, respectively, for TBW and 10% and 30% for WBW. The spatial distribution characteristics of the CUSLE deviated from the annual C-factors averaged from 2002 to 2013. C-factors were found to increase in Lower Tampa Bay, indicating the presence of deforestation as well as land management activities in the Lower Tampa Bay region [22]. The distribution of KUSLE for Lower and Middle TBW was underestimated in comparison with the developed annual K-factors by 45% and 35%, respectively. In contrast, it was overestimated by around 20% and 28% in the northwest and central portions of the WBW, respectively. In both watersheds, no significant annual dynamics were observed in K-factor between 2002 and 2013, despite the noticeable deviation of annual Kfactors from the KUSLE. The K-factor from the USLE equation, KUSLE, considered spatial changes with each HRU while neglecting the temporal dynamics with each year, while the use of the remotely sensed MODIS data estimated accurate soil moisture content. Thus, they captured the spatial as well as annual changes in the developed K-factor function that led to improved K-factor estimation in the watershed. For the years between 2002 and 2009, a Kfactor range of 0-0.1 was observed for the Old Tampa Bay and Hillsborough Bay watershed (58% of the total watershed), showing the existence of less erodible soils that can elevate the agricultural productivity encompassing the catchments. Meanwhile, elevated K-factors ranging from 0.1-0.7 were found in the remaining portion of the watershed, The K-factor from the USLE equation, K USLE , considered spatial changes with each HRU while neglecting the temporal dynamics with each year, while the use of the remotely sensed MODIS data estimated accurate soil moisture content. Thus, they captured the spatial as well as annual changes in the developed K-factor function that led to improved K-factor estimation in the watershed. For the years between 2002 and 2009, a K-factor range of 0-0.1 was observed for the Old Tampa Bay and Hillsborough Bay watershed (58% of the total watershed), showing the existence of less erodible soils that can elevate the agricultural productivity encompassing the catchments. Meanwhile, elevated K-factors ranging from 0.1-0.7 were found in the remaining portion of the watershed, which indicates the presence of highly erodible lands and potential risks for soil loss due to erosion. On the other hand, only 20% of the total watershed area showed K-factors between 0-0.1, which leaves a significant portion of the WBW with moderate or highly erodible soils [100][101][102]. Hence it is necessary to adopt sustainable water management measures as well as erosion control strategies within or near the coastal zone for further improvement of agricultural productivity. They will result in coastal water resource improvement [103][104][105].
The trends in the C-factor and K-factor time-series from 2002 to 2013 were determined by computing the C-factors and K-factors averaged annually over all the HRUs of the watershed ( Figure 6). The trend analysis results determined a decreasing trend of the watershed C-factors from 2002 to 2013, which confirms the existence of annual management-induced land cover dynamics at the watershed scale [106][107][108][109][110][111]. This highlights the importance of more detailed temporal assessments of crop and cover management for improved investigation of water management in watersheds. watershed C-factors from 2002 to 2013, which confirms the existence of annual management-induced land cover dynamics at the watershed scale [106][107][108][109][110][111]. This highlights the importance of more detailed temporal assessments of crop and cover management for improved investigation of water management in watersheds.  An evenly transitioning trend was observed in WBW, ranging from 0.218-0.224 in the five years. However, the K-factors in the TBW had minimal changes from 0.206 to 0.210 in the eight years. These findings were substantiated by shreds of research which examined the increasingly important influence that urbanization and its associated soil losses have on the individual drainage basins of TBW and WBW [111]. Thus, the trend test results for the K-factors and C-factors adhere to the real hydrogeological conditions. The inclusion of more precise temporal conditions including seasons, months, days, and hourly intervals holds importance in this context, as these conditions can impact time-series trend detection at the watershed scale, at specific monitoring stations, and concerning point and nonpoint sources of contamination within the watersheds.

Spatiotemporal Predictions of Sediment Yields
The study illustrated the spatial and temporal variability in SY assessment using a connective algorithm of the dynamic models of the C-factor and K-factor in the SWAT model. Annual SY was generated for each HRU of the watersheds from 2002 to 2013. The annual HRU wise SY ranged from 23 ton/ha to 70 ton/ha. Overall, high sediment deposition was observed in Lower Tampa Bay and Hillsborough Bay, which indicates the accumulation of the collected sediments towards the outlet of the watershed. About 16%, 44%, and 50% of the watershed was associated with the lower, medium, and higher proportions of SY, respectively. However, keen observation was essential to identify the slight reductions in the SY in some parts of Lower Tampa Bay and Old Tampa Bay. This is supported by studies which pointed out significant increases in urban spaces, developed areas, and water cover with corresponding reductions in agricultural productivity [21]. This led to reduced water quality constituents such as suspended sediments and turbidity in TBW. WBW showed moderate yield assessment in which 70%, 22%, and 8% of the watershed was associated with the lower, medium, and higher proportions of SY, respectively. SY was observed to spatially increase from the northwest portion to the southeast portion of the WBW from 2009 to 2013 [112,113]. This shows that with variable spatial resolution (30 m to 500 m) and relatively coarse temporal resolution (annual) satellite remotely sensed data, a good estimate of SY in catchments of the southeastern United States can be made. According to Figures 7 and 8, the spatial variability of sediment generation in the HRUs of both watersheds was much more significant than the temporal variations. Therefore, the temporal variability of the K-factor, which represents the severity of soil and sediment loss, can be neglected when compared to its spatial variation. Some studies support the use of dynamic models for valid water quality assessments in catchments [114][115][116][117][118].
For the year 2007, geospatial changes in the C-factor led to significant reductions in sediment generation in HRUs of TBW. This indicates the improved correlation of the dynamic C-factor and K-factor estimates on water quality predictions compared to the USLE estimates. Hence, they strongly infer the probable association of the connective soil erosional components C-factor and K-factor in changing the water quality of sediments from year to year, especially on watershed scales [119]. Overall, the results showed that the predictions using the developed models are feasible, which can help users to reduce the uncertainty of the quantification of USLE parameters of the C-factor and K-factor.
The vital contribution of the study is the remotely sensed data-driven models proposed, which present a manageable and user-friendly determination of SY for advanced water resource management and can yield advanced spatiotemporal K, C, and SY mapping for watersheds in the United States. In this study, we used K-factor and C-factor to capture the SY dynamics, although these might not fully represent the connectivity of the contributing factors of SY in the areas. Because the catchments have similar areas, mixed land covers, similar soil types, and relative terrain features, we assumed that the whole of the catchments contributed proportionally to the sediment delivered to the outlet [120,121]. However, even in small catchments, SY at the outlet depends on other factors such as atmospheric deposition and groundwater tables, as well as extreme stormwater events. Therefore, a complete understanding of different factors and dynamic SY generation entails the integration of multifaceted environmental and hydrogeological data both at very small scales and at large watershed scales. Furthermore, the annual geospatial maps of C-factors, K-factors, and SY can represent the interannual variability of these factors on SY assessment and soil erosion modeling.

Study Limitations
The connective framework of the study, which was incorporated in SWAT and tested for the watersheds of the southeastern United States, can be applied to comparable watersheds with mixed land covers. However, extensive spatiotemporal validation of these dynamic models would bolster their global effectiveness in SY predictions and water quality modeling. Regarding the differences between the actual and estimated SY accounting for soil erosion at the spatial and temporal levels in the watersheds, the error is attributable to the lack of alternative dynamic factors, modeling techniques, interpolation methods, and real hydrogeological conditions involved in the study areas. The findings of the study also highlight the necessity of finer spatial and temporal data in real-time water quality modeling of watersheds. Further studies should investigate the connections between unexplored soil erosional factors and water quality, which will broaden knowledge regarding the variability of different water quality constituents and agricultural management in watersheds with changing space and time. Remote sensing-based approaches can be substantially constrained by the coarse resolution of remotely sensed data and the accuracy of its processing, including land cover classification and estimation of various environmental indices. Therefore, justifying a remote sensing-based approach for different study watersheds requires added assessments of how the benefits of the novel models compare to the uncertainties resulting from the fundamental limitations of remote sensing products

Study Limitations
The connective framework of the study, which was incorporated in SWAT and tested for the watersheds of the southeastern United States, can be applied to comparable watersheds with mixed land covers. However, extensive spatiotemporal validation of these dynamic models would bolster their global effectiveness in SY predictions and water quality modeling. Regarding the differences between the actual and estimated SY accounting for soil erosion at the spatial and temporal levels in the watersheds, the error is attributable to the lack of alternative dynamic factors, modeling techniques, interpolation methods, and real hydrogeological conditions involved in the study areas. The findings of the study also highlight the necessity of finer spatial and temporal data in real-time water quality modeling of watersheds. Further studies should investigate the connections between unexplored soil erosional factors and water quality, which will broaden knowledge regarding the variability of different water quality constituents and agricultural management in watersheds with changing space and time. Remote sensing-based approaches can be substantially constrained by the coarse resolution of remotely sensed data and the accuracy of its processing, including land cover classification and estimation of various environmental indices. Therefore, justifying a remote sensing-based approach for different study watersheds requires added assessments of how the benefits of the novel models compare to the uncertainties resulting from the fundamental limitations of remote sensing products and methods of their processing, which is outside the scope of this study. In addition, the temporal extents (2002-2009 and 2009-2013), which are much less than 30 years, exert some limitations for long-term studies on soil loss calculations. Nevertheless, the core emphasis of future works should be the elimination of the uncertainties involved in the present study.

Conclusions
Various conceptual frameworks of sediment yield assessments have been developed in recent years for enhanced soil erosion modeling. The existing frameworks model singular factors affecting SY; however, they overlook the connection between multiple factors affecting SY spatially and temporally and lack a manageable and user-friendly platform. This study incorporated a new framework of two dynamic factors affecting SY (C-factor and K-factor) into the SWAT model using satellite remotely sensed MODIS data of EVI, SR, AWC, and LAI, and evaluated whether dynamic models can advance spatiotemporal SY predictions in watersheds. We developed a connective algorithm that incorporates all the components required to temporally (annual) and spatially (HRU level) simulate dynamic C-factors and K-factors within the SWAT model. The connective algorithm was validated using sensitivity analysis and three performance indices to determine the impact of the dynamic C-factors and K-factors on SY predictions. The rise in the PR 2 (TBW: 0.636 (Case 2) to 0.668 (Case 3); WBW: 0.533 (Case 2) to 0.644 (Case 3)) with the inclusion of the dynamic factors of C and K showed how well the SY was predicted and that the model is not too complicated concerning performance and overfit. The results showed that the model efficiency of SY predictions increased for every dynamic factor added; for all cases, the highest model efficiency occurred when the two dynamic factors were connected (Case 3: PR 2 : 0.60-0.70, R 2 : 0.70-0.80, NS: 0.65-0.75). Minor changes in K-factors and gradual decreases in C-factors in the watersheds reflected potential land use changes that can impact sediments and agricultural water use in these regions. The spatiotemporal analysis of the SY simulations generated using the dynamic C-factors and K-factors displayed moderate to high sediment deposition at the watershed outlets. The spatial variability of sediment generation in the HRUs of both watersheds was much more significant than the temporal variations. Therefore, we recommend the application of the dynamic models at smaller spatial scales than the catchment scale to eliminate uncertainty in spatial predictions. Overall, the connective system improved the reliability of SY predictions in SWAT and created options for applying remotely sensed data unlike the conventional USLE functions existing in SWAT.
The study concluded that satellite remotely sensed data of EVI, SR, LAI, and AWC can effectively be used in spatiotemporal estimation of the dynamic C-factor and K-factor, as well as SY, within basin-scale hydrologic models. The results emphasized the greater feasibility of SY and water demand predictions in SWAT when all the satellite remotely sensed dynamic variables are implemented in the connective algorithm of the C-factor and K-factor. This can help users to reduce the uncertainty of the quantification of USLE parameters of the C-factor and K-factor with the aid of the SWAT model. The vital contribution of the study is the remotely sensed data-driven models proposed that put forth a manageable and user-friendly determination of SY and can yield advanced spatiotemporal K, C, and SY mapping for watersheds in the United States. Regarding the differences between the actual and estimated SY accounting for the soil erosion for the spatial and temporal extents in the watersheds, the error is attributable to the modeling techniques, interpolation methods, and real hydrogeological conditions involved in the study areas. The findings of the study highlight the necessity of better spatial and temporal data in real-time water quality modeling of watersheds in order to ensure hydro-environmental health. Further studies should investigate the underlying uncertainties of remote sensing-based approaches as well as the relationships between unexplored soil erosional factors and water quality in watersheds with changing space and time. Data Availability Statement: Data will be available upon request.