Modeling Water Quantity and Quality Nonlinearities for Watershed Adaptability to Hydroclimate Extremes in Agricultural Landscapes

: Changing water supplies and demands, inherent to climate ﬂuctuations and human activities, are pushing for a paradigm shift in water management worldwide. The occurrence of extreme hydrometeorological and climate events such as extended wet periods and droughts, compounded with contaminants, impair the access to water resources, demanding novel designs, construction, and management across multiple hydrologic scales and biogeochemical processes. A constraint to study-ing hydrologic and biogeochemical disturbances and improving best management practices for water quantity and quality at the watershed scale resides in the suitable monitoring, data availability, and the creation of frameworks. We hypothesize that streamﬂow and contaminants, simulated by the hydrologic model Soil and Water Assessment Tool (SWAT) and evaluated during drought and extended wet periods, capture the nonlinearities of contaminants of multiple biogeochemical complexities, indicating the adaptive abilities of watersheds. Our objectives are to (1) use rain gauge and radar data and linear regression to consolidate long-term precipitation data to simulate streamﬂow and water quality using the SWAT model in the Shell Creek (SC) watershed, Nebraska, U.S.; (2) use drought and extended wet events analytics on observed and simulated hydroclimate and water quality variables to identify SWAT’s performance; and (3) identify the temporal attributions of streamﬂow and water quality to complex biogeochemical patterns of variability. We implement a watershed modeling approach using the SWAT model forced with rain gauge and radar to simulate the intraseasonal and interannual variability streamﬂow, sediments, nutrients, and atrazine loads in the SC watershed. SWAT performance uses a calibration period between 2000 and 2005 and a validation period between 2005 and 2007. We examine the model’s ability to simulate hydrologic and biogeochemical variables in response to dry and extended wet ﬂow regimes. The hydrologic model forced by either radar or rain gages performs similarly in the calibration (NSE = 0.6) and validation (NSE = 0.92) periods. It reproduces medium ﬂows closer to the observations, although it overestimates low–ﬂows up to 0.1 m 3 /s while underestimates high ﬂows by 1 m 3 /s. The water quality model shows higher NSE for streamﬂow and sediments followed by nutrients, whereas it poorly reproduces atrazine. We conclude that seasonal changes and hydroclimate conditions led to the emergence of patterns of variability associated to the nonlinearities and coupling between processes of natural and human-origin sources. extremes, the simulation of water quantity and quality nonlinearities—as properties of complex adaptive hydrologic systems—can contribute to improve the predictability of climate-resilient water resources. and biogeochemical and hydrologic processes the of in pattern is on scenario evidences two potential drivers of the temporal patterns of variability that could also be interdependent. The ﬁrst driver emerges from the non-signiﬁcant intraseasonal differences in the atrazine loads, which could be attributed to the ability of SWAT to simulate atrazine based on the parameters used for this study (Supplementary Materials, Table S5—parameters for SWAT). The third driver can be illustrated by the differences between wet and dry periods (see boxes and whiskers), which may transport the excess atrazine used before planting corn in the spring and winter wheat.


Introduction
The rising frequency and intensity of extreme hydrometeorological and climate events (EHCEs) attributed to climate change [1][2][3] evidence the need for paradigm shifts in water resources management [4,5]. Changes in hydrologic regimes [6][7][8][9] and water quality [10][11][12][13][14] triggered by climate fluctuations also require more integrative resilience and sustainability frameworks [15][16][17][18]. Our goal is to set the basis for predicting hydrologic resilience to EHCEs such as droughts and extended wet periods, identifying the ability of modeling to simulate water quantity and quality in a watershed dominated by human activities.
Watersheds' rainfall-runoff processes contribute to the mobilization of sediments, nutrients, and organic matter from croplands. Such transport by streamflow and other fluxes within the water cycle is extended to chemicals such as pesticides and emerging contaminants [13,19,20]. The use of watershed models to evaluate hydrologic consequences of volatile climate conditions represents a valuable opportunity to improve our understanding of their functionalities and mitigate monitoring gaps [21]. These models rely on the quality of forcings and parameters for accurate simulations and predictions. Long-term precipitation, a critical forcing in hydrologic modeling, has been monitored using rain gages at specific locations expanding records across multiple locations. Yet only 10.2% out of 5948 stations worldwide recorded at least 80% of daily values in the past century [3]. With constrained periods and resolutions, although geospatially distributed, remote sensing products (radar and satellite) have improved the spatial representation of precipitation gradients. Also, rain gage and radar products have been successfully used for watershed modeling [22][23][24]. Though, recommended long-term records to assess streamflow responses to EHCEs conditions [3,25] and the lack of hydroclimate data remain challenging for hydrologic modeling. The integration of monitoring precipitation techniques (radar and rain gages) represents an alternative to improve the performance of hydrologic simulations and extend the records of hydrologic responses to EHCEs [12,[26][27][28].
On the other hand, agricultural runoff is the largest non-point source of pollution of water bodies as it transports sediments, nutrients, organic matter, and chemicals from croplands [29]. Contamination of surface waters affects aquatic biota and flora across multiple hydrologic scales. Nutrients and organic matter reduce dissolved oxygen and lead to eutrophication [30]. Contaminants of exclusive human origin, such as atrazine, are used as synthetic herbicides in cropping systems, primarily evident in the United States of America's agricultural watersheds [31,32]. The adverse effects on aquatic life include bioaccumulation, reduction in reproduction, tissue abnormalities [33], sex changes [34], and photosynthesis inhibition [35]. The simulation and monitoring of this "cocktail" of contaminants [14] represent a challenge for identifying watersheds' functioning and resilience. As with water quantity, shortages of water quality observations [10,[36][37][38] limit the model's ability to simulate the fate and transport of pollutants in the aquatic environment. This factor in water quality and quantity studies also limits our understanding of the patterns of temporal variability in climate and the hydrologic responses, which can be site-specific and complex [39]. Some of these complexities reside in the nonlinear relationships between hydrometeorological processes and the lagged effects of the water quantity/quality. Thus, integrating theoretical, lab, and numerical approaches [40] becomes a need to diagnose and predict the hydrologic systems' ability to adapt to a changing environment.
As a property of socioenvironmental complex adaptive systems, nonlinearity is described by [41], together with aggregation, hierarchy, and biodiversity. These properties represent a shred of theoretical evidence to identify the temporal adaptive patterns through the underlying water quality and quantity variation and their effects on hydrologic systems in particular and the Earth system in general [13,40,42]. Hence, our scientific question is: Can a hydrologic model such as Soil and Water Assessment Tool (SWAT) capture the temporality of the nonlinear responses of contaminants in streamflow simulated during drought and extended wet periods? We hypothesize that the simulated streamflow and contaminants, parameterized by the hydrologic model SWAT and evaluated during drought and extended wet periods, follow the nonlinearities of contaminants of multiple biogeochemical complexities, indicating the adaptive abilities of watersheds. To test this hypothesis, we formulate three objectives: (1) use rain gauge and radar data and linear regression to consolidate long-term precipitation to simulate streamflow and water quality using the SWAT model in the Shell Creek (SC) Watershed, Nebraska, U.S.; (2) use drought and extended wet events analytics on observed and simulated hydroclimate and water quality variables to identify SWAT's performance; and (3) identify the temporal attributions of streamflow and water quality to complex biogeochemical patterns of variability.
Thus, in rural settings, highly diverted, contaminated, and ungauged watersheds experience poor continuity or lack of water quantity and quality monitoring. This data-deficit scenario precludes the identification of the nonlinearities in streamflow and contaminant, together with the design and implementation of better water management practices. Furthermore, the lack of data and understanding of nonlinearities in changes in water quantity and quality opaque the identification of the compounded attributions of streamflow, sediments, nutrients, and pesticides to EHCEs, land use/land cover changes, and contaminants of human origin [14,[43][44][45].

Study Area
The Shell Creek (SC) watershed drains approximately 1200 km 2 in east-central Nebraska and is located primarily in Boone, Madison, Platte, and Colfax counties of Nebraska. The towns of Schuyler, Platte Center, Lindsay, and Newman Grove are in the watershed, and the city of Columbus is directly to the south. SC has three major tributaries in the watershed: Elm Creek, Loseke Creek, and Taylor Creek (Figure 1a). The five towns in the watershed combined have a population of 1675 people [46]. The daily discharge at the United States Geological Survey (USGS) gage station averaged 1.40 ± 5.13 m 3 /s between 1947 and 2014, and the average annual precipitation was 671 ± 159 mm/year between 1910 and 2014. The watershed is home to heavily agricultural activities. Land cover types include cultivated (78.2%), herbaceous (14.6%), and forest (1.85%); while developed urban areas cover only 4.4% of the watershed. Corn (49%) and soybean (27%) dominates land use, with a high percentage of the cropland for irrigated row crops (46%) (Figure 1b and Sup- The watershed is home to heavily agricultural activities. Land cover types include cultivated (78.2%), herbaceous (14.6%), and forest (1.85%); while developed urban areas cover only 4.4% of the watershed. Corn (49%) and soybean (27%) dominates land use, with a high percentage of the cropland for irrigated row crops (46%) (Figure 1b and Supplementary Materials, Table S1). SC crop season onsets in spring-summer and depends on deep aquifer recharge or rainfed irrigation. The counties comprising the watershed include approximately 1550 farms with over 1,050,000 head of swine, cattle, and poultry [47]. Manure from livestock facilities is typically collected and land applied to agricultural fields as a soil conditioner and fertilizer. In this region, antibiotics are commonly administered prophylactically for disease prevention as well as for disease treatment. Animals are typically raised in confined animal feeding operations (CAFO) where they are confined to production barns (swine and poultry) or feedlots (cattle).

Rain Gauge Precipitation
We used data from rain gages outside the SC watershed since no instruments following the World Meteorological Organization guidelines were found. Radar estimates to characterize the temporal variation and spatial distribution of precipitation were obtained from the National Oceanic and Atmospheric Administration (NOAA) Next Generation Weather Radar (NEXRAD). A total of 12 rain gages were identified with at least 100 year-long records (1910 to 2014). Ten stations belong to the National Climatic Data Center (NCDC) from the Global Historical Climatology Network (GHCN)-Daily database [48] and two to the High Plains Regional Climate Center (HPRCC). After screening quality data, stations with more than 90% of daily observations were kept, or otherwise disregarded. This process led to the selection of eight stations (Supplementary Materials, Figure S1).

Radar Precipitation
The National Centers for Environmental Prediction Stage IV product is a nationwide mosaic of multi-sensor (radar and gages) precipitation analyses (MPEs), produced by the 12 River Forecast Centers (RFC). It has offered real-time, hourly/six-hourly precipitation analyses since 2002 with a 4 km resolution [49]. Daily values are summed up from the six-hourly analysis and the data might have manual quality control.

Streamflows and Pollutant Loads
We accessed the National Water Information System repository of the [50] USGS 2012 using the R package "dataRetrieval" [51] to download the streamflow and water quality records for SC. The availability of records varies from streamflow daily records since 1947. However, water quality variables have less than sixty daily observations constraining the calibration and eventual performance of the model described below. Watershed models cannot be calibrated properly for a wide range of hydrologic conditions, sediments, and nutrients loads when watershed lacks monitoring data [36][37][38].
The disparity in the observations and the need for continuous time series of contaminants for the sensitivity analysis and the discussion on contaminant complexity led to the use of a linear regression to create a long-term synthetic time series. Initially, correlations between the predictor variable (streamflow) and response variables (e.g., sediments load) were implemented when data were available, between 1992 and 1994 and 2008 and 2009. The pollutant loads were calculated as the product of streamflow times its concentration. A logarithmic transformation was identified to describe the correlation (Equation (1)). Secondly, a linear model was fit to the Log10 of sediments, nutrients, and atrazine loads. The performance of the regression was evaluated by the multiple R 2 and the p-value of F-statistics. With multivariate regression, each response variable follows its regression equation, as well as the covariance that exists among each pairwise response variables. For instance, multivariate regression would consider the relationship between the nitrogen and phosphorus levels in the water, if one exists. The flows observed in the monitoring periods of water quality (Table 1) were randomly selected as a representative sample of streamflow in SC, ranging from 0.23 to 200 m 3 /s.

The Shell Creek Model
SWAT is a semi-distributed and physically based model that simulates hydrologic processes (i.e., infiltration, overland flow, evapotranspiration, shallow and groundwater), crop growth, and fate and transport of nutrients and chemicals in soils and surface waters [52][53][54]. We used ArcSWAT, the ArcGIS-ArcView extension, and a graphical user interface for the Soil and Water Assessment Tool-SWAT2012, Revision 637 under the Ar-cGIS version 10.2.2. SC watershed was discretized in Hydrological Response Units (HRU), overlapping soil physical properties, land use, management practices, and topography. After a precipitation event, similar hydrologic response units (HRU) generate the same amount of runoff to route throughout the sub-basin streams.
SWAT can simulate the physicochemical processes before water runoff from landscapes. Nitrogen and phosphorus cycles are represented like soil processes where elements move in and out of pools. These include nitrification mineralization of nitrogen or adsorption, phosphorus fixation, and plant uptake. Inputs can be point and non-point sources. Nutrient cycles are simulated after the Environmental Policy Integrated Climate (EPIC) model [55,56]. On the other hand, soil erosion is estimated using the Modified Universal Soil Loss Equation (MUSLE) [57], which is a modification of the Universal Soil Loss Equation (USLE) [58].
Fate and transport of pesticides involve degradation, infiltration, and volatilization, whose algorithms in SWAT were taken from the EPIC model. Pesticides can be moved in the solid (soil phase) or dissolved in the solution (liquid phase). SWAT model transports soluble pesticides via surface runoff, lateral flow, or infiltration as a sorbed pesticide is moved by surface runoff only to the main channel. SWAT simulates crop rotation and best management practices.
Additionally, SWAT simulates the in-stream processes of nutrients, pesticides, bacteria, and heavy metals using the Enhanced Stream Water Quality Model (QUAL2E) model [59]. The algorithm estimates 1-D longitudinal and temporal variation of water quality in-stream. We use an adaptation of GLEMS FOR SWAT for the simulation of agricultural chemicals, according to [60]. Refer to SWAT's documentation for a more detailed description of the process and equations for water quality [54].

Input Data
The SC model was forced independently with rain gage and radar precipitation estimates. The spatial and temporal attributions of radar-based and rain gauge precipitation estimates, respectively, are used to validate the spatial and the long-term temporal patterns of variability. Since SWAT assigns the weather station closer to the mass centroid of each subbasin for water balance calculations, six rain gages and 40 grid centroids-treated as "virtual rain gages"-were used. A global sensitivity analysis and the calibration-validation of the two models were compared.
Spatially distributed soil types, land use, topography, and climate data were input to the model. Soil types were obtained from the digital US General Soil Map called the STATSGO2 database, which is pre-installed in ArcSWAT. Land use-at a spatial resolution of 30 m cell size raster-was downloaded from the United States Department of Agriculture (USDA)-National Agricultural Statistics Service Information (NASS) [47]. We also used the USGS National Elevation Dataset (NED) at 1 arc-second (30 m) [61,62], the USGS National Hydrography Dataset (NHD), and the USDA-Natural Resources Conservation Service's (NRCS) Watershed Boundary Dataset (WBD), downloaded from The National Map [63]. The NHD and WBD guided the SC streamflow network delineation. The threshold for flow accumulation was set to 2231 ha, which equals 30,000 upslope grid cells. A total of 54 sub-basins were delineated with an area mean of 22.3 km 2 ( Figure 1). Meteorological records-since 1988-including relative humidity, solar radiation, air temperatures, and wind speed, and were obtained from the two HPRCC stations. The Soil Conservation Service (SCS) curve number method [64], an empirical model commonly used in the U.S.A. since the 1950s, estimated the surface runoff. The channel routing and evapotranspiration flow were calculated by the Muskingum method [55] and the Penman-Monteith approach [65], respectively.
In Equation (2), Q sur f is the accumulated runoff in millimeters (mm), R day is the daily rainfall depth, I a is the initial abstraction including surface storage, interception, and infiltration prior runoff, and S is the retention parameter, shown in Equation (3). The S parameter depends on soils, land use, management, and slope and soil water content.
where CN is the curve number.

Crop Management
Terracing was set as a management practice in SC. Irrigated crops and dryland are distributed equally across the watershed (38%) (Supplementary Materials, Table S1). The auto-irrigation was assumed from a deep aquifer after [39] and set randomly to the croplands in half of the subbasins, whereas the other land use was set as drylands. Table 2 presents the management operations set up to preserve the distribution of corn (49%) and soybean (27%) in the SC model, estimated from the 2005 Nebraska land-use patterns map at 30-m cell size resolution according to the Center for Advanced Land Management Information Technologies [66]. A common conservation practice in the US corn-belt is crop rotation that might involve corn-soybean or wheat-corn [67]. Corn was allocated randomly to 24 subbasins, corn-soybean rotation in 30 subbasins (27%), and soybean exclusively rotated to corn in 30 sub-basins (27%). Liquid manure from feeding animal lots (CAFO), commonly applied to cultivated land, was used on any cropland in SC. The application rates of fertilization on corn and soybean ranged between 15 and 90 kg/ha, and the atrazine application rate was 1 kg/ha after previous studies [10,39,67]. The pesticide was assumed to be homogeneously applied across the study area.

Parameter Selection
We performed a global sensitivity analysis (GSA) of the SC model using the sequential Uncertainty Fitting ver. 2 (SUFI-2) algorithm in the SWAT-CUP software version 5.1.6 [68]. A t-test estimates the parameter sensitivity, while the p-value is used for the significance. The most sensitive parameters had larger t-test absolute values and significant sensitivity (p-value ≤ 0.1). SUFI-2 identifies the sources of uncertainty using two factors: the P-factor that quantifies the 95% prediction uncertainty (95PPU); the d-factor that estimates the mean thickness-from 0 to infinity-of the 95PPU band divided by the standard deviation of the observed data. If the P-factor and d-factor tend to be one and zero, respectively, it also indicates that simulations match observations. We examined the sensitivity of 31 parameters associated with streamflow generation (Supplementary Materials, Table S1). A Latin Hypercube Sampling (LHS) generated 150 parameter combinations within initial ranges for a year-long simulation.

Calibration and Validation
Following the GSA, the hydrologic model was calibrated and validated in SUFI-2. We used radar data from 2002 to 2013; the 12 year-long records were split into seven-year (2002-2008) and five-year (2009-2013) periods for calibration and validation, respectively. Studies suggest a proper streamflow calibration to ensure a better water quality simulation [36,37]. The water quality model was calibrated for 2002-2010 using the SUFI-2 procedure in SWAT-CUP. Since runoff influences the fate and transport of sediments, this was first calibrated, followed by nutrients, and atrazine loads. Calibration started with parameter selection reported in the literature [10,29,[36][37][38][39][67][68][69][70] and SWAT documentation [53]. A total of 400 parameter combinations were generated using the LHS for each variable. The use of the SUFI-2 procedure followed two steps. First, an LHS generated a set of 400 parameter combinations within the initial ranges of the most sensitive parameters. After the simulations ran, the SUFI-2 procedure narrowed the range limits, which reduced the uncertainty associated with the value of parameters. Since the new ranges might fall beyond the initial values, these were visually checked and manually adjusted when needed. The target objective function was the NSE, which identified the best simulation. Second, the "best" ranges set another iteration. Next, we generated another 400 parameter combinations and ran the respective simulations. Here, if the NSE improved at least 5%, iterations would continue passing the "best" ranges to the next iteration. Otherwise, iterations should stop, and the ranges kept for the validation without further change. In the validation step, another precipitation time frame was input. We again used 400 parameter combinations within the "best" ranges from the last iteration of the calibration step and ran the model. These steps continued for nutrients and atrazine. No validation was performed for water quality estimates. Once the most sensitive parameters were fitted using the linear regression, the SC model was forced into a long-term simulation (1910 to 2010).
To evaluate the model performance, the Nash-Sutcliffe efficiency (NSE) coefficient [71] compared monthly observed and simulated streamflow. NSE can range from −∞ to 1. A model performs better if it tends to or equals 1, but the performance decreases when it falls below 0. Complementary metrics of efficiency include the root mean square error (RMSE) and the standard deviation estimated for observed streamflow. Also, the percent bias (PBIAS) and the ratio of the root mean square error (RMSE) and the standard deviation (RSR) were also calculated. For the simulation i, time steps t varies from 1 to n; O and S are the observed and simulated streamflow, respectively; and O is the O average. A negative or positive PBIAS (%) evidence an overestimated or underestimated streamflow, and its absolute value quantifies the bias magnitude. RSR can become any positive number from zero, with lower values suggesting a better model.

Hydrologic Model Implementation and Analyses
Once the model was calibrated, including expanding the water quality time series, we ran long-term simulations of streamflow and water quality variables. Then, using the simulations we identified the drought and extended wet spells on water quantity and quality in SC. Next, we provided a synthesis of the model implementation process, the description of the estimations of EHCEs, and the development of the complexity framework.

Hydrologic and Water Quality Simulations
The SC flows have complete records since 1950, suitable to validate long-term simulations. After rain gauge data and the associated simulations were validated for the S.C. watershed, the model was forced into long-term simulation (1950 to 2010). The "best" calibrated ranges were kept similar to the validation step, and the 400 parameters-combinations were reused. Then, the simulation was run and model performance was evaluated.

Drought and Extended Wet Conditions
The three-month Standard Precipitation Index (SPI) identified years with extreme hydrometeorological and climate conditions: drought and extended wet periods in SC [62]. We used rain gauge data since they offer time frame records longer than radar. SPI values below zero indicate dry conditions, ranging from mild (−0.99 to 0) to moderate (−1.49 to −1.00) to severe (−1.99 to −1.50) to extreme (≤−2) drought [72]. Conversely, one can consider SPI values greater than zero from mild (>0 to 0.99) to extreme (>2) wet conditions. This study assumed a very wet condition when SPI was above 1 and with a very dry condition below −1.
On the other hand, a statistical approach was used to analyze the flow conditions of the SC watershed. We fitted the streamflow data into a gamma distribution, and the probability density function was estimated using the observed and the simulated flows. The flow percentiles were also estimated using the inverse gamma cumulative distribution function in MATLAB. Flow conditions were identified as high flows assumed to be equal to or greater than the 90th percentile flow; the medium flows ranging between the 20th and 80th percentile, and low flows those lower than the 10th percentile.

Results and Discussion
The data and information analyses in this section are organized and synthesized to support the thesis that extended wet and drought conditions in the SC watershed expand/diminish the concentrations of contaminants associated with agricultural activities. We started reviewing rainfall's monthly and seasonal climatologies as evidence of rain gauge and radar suitability for the assessment of the temporal and geospatial patterns of variability, respectively. Then, we assessed the consolidation of long-term time series and SWAT's parameter sensitivity and performance analyses simulating long-term streamflow and water quality. In each section is discussed the temporal patterns of streamflow vari-ability and water quality in response to drought and extended wet events and whether the nonlinear biogeochemical attributions unveil SC as a complex adaptive watershed.

Hydroclimate Variability
This study identified 12 rain gages with records at least 100 years long (1910 to 2014). Ten stations belong to the NCDC from the GHCN-Daily database [48] and two to the HPRCC. As a watershed dominated by agriculture (>75%), SC is located within the wettest of four hydroclimate regions in Nebraska [73], where rainfed agriculture relies on the precipitation onset in April, considered the month with the lowest variance, and the peak of precipitation in June (Figure 2). The author in [73] identified a similar monthly hydroclimatology with a monthly high of June, indicating a consistent temporal hydroclimate pattern of variability in Elkhorn River Basin's watersheds such as SC. The regionalization of precipitation along the Platte River basin captures a series of large-scale patterns of variability, but reports unclear evidence of the occurrence and alteration of watershed hydrologic processes due to the low resolution of gridded products or the absence of rain gauges. Figure 1b shows constrained biodiversity, also evidenced by the intense agricultural activity for a drainage area of approximately 1200 km 2 (49 and 27% of the surface used to grow corn and soybean, respectively). Furthermore, Figure 3 shows a homogeneous distribution of spring and summer radar-based precipitation, reflecting the lack of emerging patterns in response to aggregation in SC. Biodiversity and aggregation are properties of complex adaptive systems, as [41] suggests. Yet, there was a clear difference between wet and dry years (Table 3). Wet-summer precipitation is almost 40% larger than the dry-summer precipitation, representing an incentive to improve sub-seal to seasonal forecasts of drought and extended wet conditions. Further, the estimation of SPI for the SC indicates that summers are wetter and winters drier based on the SPI's trend analysis (Figure 4). The regionalization of precipitation along the Platte River basin captures a series of large-scale patterns of variability, but reports unclear evidence of the occurrence and alteration of watershed hydrologic processes due to the low resolution of gridded products or the absence of rain gauges. Figure 1b shows constrained biodiversity, also evidenced by the intense agricultural activity for a drainage area of approximately 1200 km 2 (49 and 27% of the surface used to grow corn and soybean, respectively). Furthermore, Figure 3 shows a homogeneous distribution of spring and summer radar-based precipitation, reflecting the lack of emerging patterns in response to aggregation in SC. Biodiversity and aggregation are properties of complex adaptive systems, as [41] suggests. Yet, there was a clear difference between wet and dry years (Table 3). Wet-summer precipitation is almost 40% larger than the dry-summer precipitation, representing an incentive to improve sub-seal to seasonal forecasts of drought and extended wet conditions. Further, the estimation of SPI for the SC indicates that summers are wetter and winters drier based on the SPI's trend analysis (Figure 4).

Sensitivity Analysis and Model Calibration
We contrasted radar-based and rain gage-based precipitation data to validate the long-term time series of the latter (Supplementary Materials, Figure S1). The implementa-

Sensitivity Analysis and Model Calibration
We contrasted radar-based and rain gage-based precipitation data to validate the long-term time series of the latter (Supplementary Materials, Figure S1). The implementa-

Sensitivity Analysis and Model Calibration
We contrasted radar-based and rain gage-based precipitation data to validate the longterm time series of the latter (Supplementary Materials, Figure S1). The implementation of this approach led to two findings. The first emerged from SWAT's global sensitivity analysis of both digital resources. In Figure 5, the most sensitive parameters (p < 0.10) were the runoff curve number (CN2), the fraction of field capacity water content (FCB), effective hydraulic conductivity in the main channel alluvium (CH K2), and the average slope of tributary channels (CH S1) for both, radar and rain gauge-based forcings. Since the runoff volume strongly depends on the selected method for calculation and the SCS curve number, one can expect the CN2 to be a very sensitive parameter. The authors in [12,26] reported similar SWAT sensitivities to CN2 for both radar and rain gauge-based forcings. Further, the sensitivity of the radar-based forced model to ESCO led to an annual ET lower than that simulated by the rain gage-forced model (Figure 5c,d). While the average ET was similar for both, some hydrologic units for rain gauge-based forcings were 10 to 20 mm higher than those for radar. Also, more radar-based hydrologic units experience low ETs, which might be due to the homogeneous distribution of rainfall and irrigation practices across the watershed (Figure 3b,c and Figure 5d). Infiltration and subsurface storage have been found to be sensitive parameters to adjust water budgets [76]. Thus, calibrating the ESCO and FFCB parameters balanced the evaporative demands in SC from deeper soil layers to match the simulated and observed streamflow. The second finding is how the consistency in the sensitive parameters captures the nonlinearities in lumped hydrologic processes, whether spatially distributed or punctual forcings force them. Yet, the disaggregated expression of radar-based forcings might indicate multiple geospatial attributions of the underlying hydrologic processes [75]. In agriculturally dominated landscapes, refs. [73,77] identified shorter recovery times of remote sensing-based leaf area index values after the incidence of droughts and extended wet events.
were the runoff curve number (CN2), the fraction of field capacity water content (FCB), effective hydraulic conductivity in the main channel alluvium (CH K2), and the average slope of tributary channels (CH S1) for both, radar and rain gauge-based forcings. Since the runoff volume strongly depends on the s e l e c t e d method for calculation and the SCS curve number, one can expect the CN2 to be a very sensitive parameter. The authors in [12,26] reported similar SWAT sensitivities to CN2 for both radar and rain gauge-based forcings. Further, the sensitivity of the radar-based forced model to ESCO led to an annual ET lower than that simulated by the rain gage-forced model (Figure 5c,d). While the average ET was similar for both, some hydrologic units for rain gauge-based forcings were 10 to 20 mm higher than those for radar. Also, more radar-based hydrologic units experience low ETs, which might be due to the homogeneous distribution of rainfall and irrigation practices across the watershed (Figures 3b,c and 5d). Infiltration and subsurface storage have been found to be sensitive parameters to adjust water budgets [76]. Thus, calibrating the ESCO and FFCB parameters balanced the evaporative demands in SC from deeper soil layers to match the simulated and observed streamflow. The second finding is how the consistency in the sensitive parameters captures the nonlinearities in lumped hydrologic processes, whether spatially distributed or punctual forcings force them. Yet, the disaggregated expression of radar-based forcings might indicate multiple geospatial attributions of the underlying hydrologic processes [75]. In agriculturally dominated landscapes, refs. [73,77] identified shorter recovery times of remote sensing-based leaf area index values after the incidence of droughts and extended wet events.  The performance of the SWAT model to the forcings confirms the scenario of parameter consistency in the sensitivity analyses. The efficiency metrics for the calibration/validation for radar and rain gauge-based forcings were similar. The ability of SWAT to simulate streamflow was practically the same (Supplementary Materials, Table S4), which also indicates the reliability of the streamflow simulations. Calibration efficiency indices were similar to efforts using monthly rain gauge-based data for SWAT [24,28,78] and radar-based simulations [22,23,26,79,80]; however, the validation efficiency indices were higher than those for calibration. NSE and RSR metrics for calibration were 0.58 and 0.65, respectively; for validation, the values were 0.92 and 0.28, respectively. This response was observed in the reduced streamflow temporal variability during the validation period (2009 to 2013), the model overestimated low-flows below the tenth percentile, while underestimated high-flows, especially above the 99th percentile ( Figure 6). The SC model improved the simulation of different flow regimes in a forest-agricultural area (89%) where a rain gageforced model constantly underestimated low, medium, and high flows at the 5th, 50th, and 95th percentiles, respectively [12]. The improved model performance during the validation might be aligned with shifts in precipitation like those reported by [81]. Further analyses may indicate whether this watershed is transitioning to another hydrologic regime attributed to the intensification of agriculture or the intensification of droughts or floods in the region [44]. Thus, these performance metrics make the SWAT model suitable to capture the nonlinearities between hydrologic and biogeochemical processes and their response to hydroclimate spells such as drought and extended wet events. ydrology 2022, 9, x FOR PEER REVIEW 13 of 2 also indicates the reliability of the streamflow simulations. Calibration efficiency indice were similar to efforts using monthly rain gauge-based data for SWAT [24,28,78] and ra dar-based simulations [22,23,26,79,80]; however, the validation efficiency indices wer higher than those for calibration. NSE and RSR metrics for calibration were 0.58 and 0.65 respectively; for validation, the values were 0.92 and 0.28, respectively. This response wa observed in the reduced streamflow temporal variability during the validation period (2009 to 2013), the model overestimated low-flows below the tenth percentile, while un derestimated high-flows, especially above the 99th percentile ( Figure 6). The SC mode improved the simulation of different flow regimes in a forest-agricultural area (89% where a rain gage-forced model constantly underestimated low, medium, and high flow at the 5th, 50th, and 95th percentiles, respectively [12]. The improved model performanc during the validation might be aligned with shifts in precipitation like those reported by [81]. Further analyses may indicate whether this watershed is transitioning to another hy drologic regime attributed to the intensification of agriculture or the intensification o droughts or floods in the region [44]. Thus, these performance metrics make the SWAT model suitable to capture the nonlinearities between hydrologic and biogeochemical pro cesses and their response to hydroclimate spells such as drought and extended wet events

The Statistical Model and the Construction of Long-Term Time Series
The authors in [13,14] identified temporal patterns of variability in water quality and quantity attributed to contrasting hydroclimate conditions. Their work also evidenced how droughts and extended wet events regulate the intraseasonal and inter-annual wate quality variations. Such variations follow the periodicity of the underlying hydrologica and biogeochemical processes in some cases. In the SC, as in many ungauged watershed worldwide, the robustness of those records grapples with the lack, discontinuity, or short ness of time series, constraining the development of integrated water management plans We mitigated the conspicuous deficit of water quality data by using a linear regression model (Equation (1)) that enabled the enhancement of the observed sediment, nitrate phosphate, and atrazine records in SC's streamflow. The R 2 (p < 0.01), observed in Tabl 4, suggests that the logarithmic transformation of the streamflow and the water quality variables explained between 75 to 89% of the total variability of phosphorous, nitrogen phosphorous, and sediments; for atrazine, the model explained 56% of the total variability (Supplementary Materials, Figure S2). The authors in [14] proposed the term "chemica cocktails" to describe the hydrological connection of emerging contaminants. This con

The Statistical Model and the Construction of Long-Term Time Series
The authors in [13,14] identified temporal patterns of variability in water quality and quantity attributed to contrasting hydroclimate conditions. Their work also evidenced how droughts and extended wet events regulate the intraseasonal and inter-annual water quality variations. Such variations follow the periodicity of the underlying hydrological and biogeochemical processes in some cases. In the SC, as in many ungauged watersheds worldwide, the robustness of those records grapples with the lack, discontinuity, or shortness of time series, constraining the development of integrated water management plans. We mitigated the conspicuous deficit of water quality data by using a linear regression model (Equation (1)) that enabled the enhancement of the observed sediment, nitrate, phosphate, and atrazine records in SC's streamflow. The R 2 (p < 0.01), observed in Table 4, suggests that the logarithmic transformation of the streamflow and the water quality variables explained between 75 to 89% of the total variability of phosphorous, nitrogen, phosphorous, and sediments; for atrazine, the model explained 56% of the total variability (Supplementary Materials, Figure S2). The authors in [14] proposed the term "chemical cocktails" to describe the hydrological connection of emerging contaminants. This conceptualization emulates another property of complex adaptive systems proposed by [41]: aggregation. The aggregation in ecosystems and "chemical cocktails" can lead to patterns unseen in compounds individually due to interactions between biogeochemical or resource management (i.e., fertilization in cropping systems vs. natural occurrence of nutrients). Following the discussion on SWAT's ability to capture the nonlinearities of hydrological and biogeochemical processes, we suggest that the logarithmic regression model also evidences the gradients of biogeochemical complexity in the SC streamflow. Thus, [13,14] suggest that the compounding of chemicals and hydroclimate events supports emerging patterns such as increasing contaminant loads. On the other hand, the construction, calibration, and analysis of a long-term time series of streamflow and water quality variables simulated by the SWAT model can provide some insights into SWAT's ability to simulate water quality, the biogeochemical responses to drought, and extended wet events, and ultimately, provide some insights to whether SC is a climate-resilient complex system. SWAT's model calibration parameters (Supplementary Materials, Tables S5 and S6) led to NSE values for streamflow, sediments, nitrogen, phosphorus, and atrazine (0.78, 0.82, 0.61, 0.72, and 0.01, respectively). These values coincide with those reported by [10,34,39] for streamflow and sediment simulations in watersheds across the Northern High Plains. Once again, the gradient of biogeochemical complexity based on the NSE values may illustrate the interdependency of the contaminants with the natural and the human origin sources and the alterations of the water and biogeochemical cycles in SC waters. Further, the reduced performance for contaminants of a single human origin source can be attributed to unknown processes (e.g., chemical transport), socioenvironmental processes (e.g., landslides, irrigation, and water diversion), or constrained model parameterizations, which lead to the propagation of uncertainties in the model outputs [36,[67][68][69]. Hence, droughts and extended wet events can make more evident the contrasting hydrological and biogeochemical patterns of temporal variability and their origin.

Extreme Conditions
Section 3.1 used the SPI to identify the interannual trends of drought and extended wet events during summer and winter (Figure 4). In this section, we characterize the drought and extended wet months by fitting the long-term precipitation records into a gamma distribution function (Figure 7). Thus, months with percentiles above 90% and below 10% were considered under drought and extended wet months, and the criteria used to select water quality values in any given month. We applied the Tukey-Kramer honest significant difference to test differences between summer and winter and whether wet and dry conditions were different for streamflow and loading pollutants (Figure 8). Consistent with the sensitivity and calibration/validation results, three patterns emerged from the aggregation and nonlinearity associated with streamflow and the biogeochemical gradient (see the letters at the bottom of each plot in Figure 8). In the first pattern, precipitation is the driver of intraseasonal variability of streamflow. Given the size of the watershed and the strong irrigation activity, streamflow presents the largest variability and magnitude. In the second pattern, sediments, nitrogen, and phosphorous present similar patterns differentiating wet summers from the rest of the conditions. Yet, the conspicuous high ni-trogen and phosphorous loads during the wet summers are likely attributed to fertilization practices and other biogeochemical and hydrologic processes that foster the transport of contaminants in the streams. The third pattern is based on atrazine's lack of response to intraseasonal climate variability (i.e., the wet and dry or the winter and summer periods). This scenario evidences two potential drivers of the temporal patterns of variability that could also be interdependent. The first driver emerges from the non-significant intraseasonal differences in the atrazine loads, which could be attributed to the ability of SWAT to simulate atrazine based on the parameters used for this study (Supplementary Materials, Table S5-parameters for SWAT). The third driver can be illustrated by the differences between wet and dry periods (see boxes and whiskers), which may transport the excess atrazine used before planting corn in the spring and winter wheat.
x FOR PEER REVIEW 15 of 21 to fertilization practices and other biogeochemical and hydrologic processes that foster the transport of contaminants in the streams. The third pattern is based on atrazine's lack of response to intraseasonal climate variability (i.e., the wet and dry or the winter and summer periods). This scenario evidences two potential drivers of the temporal patterns of variability that could also be interdependent. The first driver emerges from the nonsignificant intraseasonal differences in the atrazine loads, which could be attributed to the ability of SWAT to simulate atrazine based on the parameters used for this study (Supplementary Materials, Table S5-parameters for SWAT). The third driver can be illustrated by the differences between wet and dry periods (see boxes and whiskers), which may transport the excess atrazine used before planting corn in the spring and winter wheat.

Conclusions
The thesis presented here states that hydrologic and biogeochemical simulations using the SWAT model evidence the abilities of the Shell Creek' watershed's adaptability to extended wet and droughts. We supported this thesis on three premises. First, we suggest that the temporal patterns of variability in the hydroclimate based on long-term time series of rain gauge-based precipitation aggregate the distribution of radar-based rainfall. Second, parameter sensitivity, calibration, and validation analytics are used to improve SWAT's model performance to simulate temporal changes in hydrologic and biogeochemical variables (e.g., sediments, nitrogen, phosphorous, and atrazine loads) those responses from extended wet and drought events. Third, the approaches above and the resultant hydrologic and biogeochemical simulations evidence the nonlinearity and aggregation properties of complex adaptive systems defined by [41].
The hydroclimate, sensitivity analyses, and calibration/validation resultant from the comparison of rain-gauge and radar-based time series data indicate that both digital resources are comparable. Since the latter spans a significantly larger time series, it was selected to calibrate the SWAT model and used as forcing for long-term hydrologic and biogeochemical simulations.
We used a log-transformed regression model to enhance a short-term time series of biogeochemical variables to calibrate the SWAT model and create long-term simulations of sediments, nitrogen, phosphorous, and atrazine loads in SC. The results were encouraging for all but atrazine.
The hydrologic and biogeochemical simulations were aggregated for wet and dry conditions during the summer and winter. The simulations identified three patterns attributed to agriculture's intensification and extent in the area: water availability for irrigation (during dry periods), the absence of a growing season, and the decrease of agronomic activities during the winter. The first pattern evidenced the influence of the rainy season and irrigation on the summer streamflow, with large flows in both, the wet and dry summers. During the winter, the lowest flows were observed. The second pattern is illustrated by the difference between wet summers and the rest of the seasons and dry conditions,

Conclusions
The thesis presented here states that hydrologic and biogeochemical simulations using the SWAT model evidence the abilities of the Shell Creek' watershed's adaptability to extended wet and droughts. We supported this thesis on three premises. First, we suggest that the temporal patterns of variability in the hydroclimate based on long-term time series of rain gauge-based precipitation aggregate the distribution of radar-based rainfall. Second, parameter sensitivity, calibration, and validation analytics are used to improve SWAT's model performance to simulate temporal changes in hydrologic and biogeochemical variables (e.g., sediments, nitrogen, phosphorous, and atrazine loads) those responses from extended wet and drought events. Third, the approaches above and the resultant hydrologic and biogeochemical simulations evidence the nonlinearity and aggregation properties of complex adaptive systems defined by [41].
The hydroclimate, sensitivity analyses, and calibration/validation resultant from the comparison of rain-gauge and radar-based time series data indicate that both digital resources are comparable. Since the latter spans a significantly larger time series, it was selected to calibrate the SWAT model and used as forcing for long-term hydrologic and biogeochemical simulations.
We used a log-transformed regression model to enhance a short-term time series of biogeochemical variables to calibrate the SWAT model and create long-term simulations of sediments, nitrogen, phosphorous, and atrazine loads in SC. The results were encouraging for all but atrazine.
The hydrologic and biogeochemical simulations were aggregated for wet and dry conditions during the summer and winter. The simulations identified three patterns attributed to agriculture's intensification and extent in the area: water availability for irrigation (during dry periods), the absence of a growing season, and the decrease of agronomic activities during the winter. The first pattern evidenced the influence of the rainy season and irrigation on the summer streamflow, with large flows in both, the wet and dry summers. During the winter, the lowest flows were observed. The second pattern is illustrated by the difference between wet summers and the rest of the seasons and dry conditions, likely produced by decreased cropping activities and the constrained use of fertilizers during droughts. The third pattern evidences the low SWAT model performance simulating atrazine based on the parameter estimation and selection proposed in this study (not the architecture and parameterizations of the model per se).
Finally, we used the shreds of evidence above in Section 3 to support the conceptualization of SC as a complex adaptive system based on [41] properties of aggregation and nonlinearity. The consistency in the sensitivity analysis, indicating that both radar and rain-gauge-based forcings led to the similar temporal variability of streamflow, but the different distribution of ETs indicates the relevance of seasonal changes and aggregation of spatial distribution patterns. Furthermore, the biogeochemical complexity gradient built upon the differences among sediments, nitrogen, phosphorous, and atrazine simulations and model performances indicates the gradual increase of human activities from sediment load to nitrogen and phosphorous to atrazine.
Future work on properties of complex adaptive systems such as hierarchy and biodiversity my lead us to identify potential sources of predictability for hydrologic shifts [82] and management practices and policies that favor hydrologic resilience [17] across multiple scales [73].
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/hydrology9050080/s1. Figure S1. Rain gage stations and radar grid centroids used as forcings of the Shell Creek model in SWAT; Figure S2: Logarithmic regression between sediments, nutrients, and atrazine, and streamflow in Shell Creek. The black line is the fitted regression. Red dashed lines display the 95% confidence band. Sed: sediments; Atrazine; TN: total nitrogen; TP: total phosphorus; Q: streamflow; Table S1. Distribution of irrigation and land use in Shell Creek; Table S2: Periods considered for trend analysis in summer (S) and winter (W) at precipitation stations near Shell Creek; Table S3: Precipitation trend analysis in summer (S) and winter (W) at rain gage stations near Shell Creek; Table S4: Performance of calibrated and validated models using rain gage and radar forcings; Table S5. Daily and monthly performances of the water quality model for Shell Creek. Table S6. SWAT parameter values for atrazine.