Most farmers in the survey do not use irrigation, but rather rely on rainwater to cultivate cereal crops. Many farmers do view irrigation as an important opportunity to improve their agricultural system but lack the capacity or knowledge to access and utilized surface and groundwater sources. In the irrigated areas surveyed in the Rift Valley (at Arata and around Lake Ziway), farmers produce horticultural crops that include tomato, onion, pepper and cabbage. The farmers that use irrigation typically achieve better productivity and profits, although they also report occasionally losing their crops due to excessive rain, or lack of profit due to price reductions during productive periods. As demonstrated in the previous section on water quality and as shown in Figures 1 and 2, high quality water sources do exist in the region that could be or were already being used for irrigation. For instance, the existing and on-going development of the floriculture sector in the MER (such as Meki and Ziway towns in the MER) is mainly based on surface water (Lake Ziway and its rivers) resources. Most groundwater sources in the region are too poor in quality for agricultural use, however (e.g., Figure 1; red block area), and increased groundwater extraction may also threaten the sustainability of the aquifer, which is essential for many of the region's rural drinking water supplies. With regard to drinking water, these sources have also previously been shown to be contaminated with elevated levels of fluoride and arsenic [38,39,57]. Continued monitoring of the quantity and quality of the groundwater resources is essential to mitigate the associated negative impacts.

#### *5.3. Perceptions of Climate Stresses and Climate Change*

Most farmers reported in the survey that lack of rain is one of the main constraints to agricultural productivity. Many also indicated that they thought rainfall was becoming more erratic. Figure 5 shows farmers' perceptions of changes in the pattern of rainfall and temperature. More than 70% of the farmers stated that rain comes and stops later than expected, compared to ten years ago. Moreover, a majority perceives that the climate is

drier, hotter, less predictable and generally worse for farming than it used to be, especially in lowland areas. A few farmers (6%) report no changes in temperature and rain. Despite perceptions of more erratic precipitation, some farmers report that their farm output is improving due to their use of improved seeds and fertilizers.

**Figure 5.** Farmer perceptions of climate change in the Ziway–Shala basin.

‐ In focus groups, it was noted that farmers were generally aware that changes to the climate are becoming a serious problem. These changes are widely perceived as being related to rainfall variability (rain comes too early and stops too late, and is uneven in amount, whether too high or too low), and farmers in the Rift especially reported experiencing a significant reduction in overall rainfall, alongside extreme and unexpected, rare heavy rainfall episodes. In the highlands, complaints were primarily related to variability, rather than to the quantity of rainfall. Due to rainfall shortage in the Rift, farmers start planting drought-tolerant early maturing maize varieties, though there is limited supply of such seeds. In some cases when there has been no rainfall in April and May, farmers forgo planting entirely. Maize productivity is also decreasing from time to time; however, farmers are achieving enhanced returns by shifting to wheat, barley and teff cultivation.

‐

‐

‐

#### *5.4. Factors Determining Good and Bad Crop Yield*

Farmers pointed out that uninterrupted and sufficient quantities of rain, availability of farming inputs (fertilizer and seeds), and personal efforts are the most important factors determining a good crop season and productivity (Figure 6). In contrast, lack of or excessive rain or unexpected rain, low crop price, lack of available and affordable inputs, and insufficient personal effort all contribute negatively to crop productivity and profit (Figure S2). Rainfall appears most important in these responses. The study revealed that insect pests and diseases are also major problems facing farmers in the region.

**Figure 6.** Relative importance of factors determining a good crop year.

#### *5.5. Climate Change Impacts on Agriculture, Adaptation Practices and Constraints*

‐ ‐ In the 5 years prior to the survey, 66.7% of the farmers reported experiencing crop failure at least once. Crop failure occurred in the Rift floor (<1750 m) an average of 1.6 times in a 5-year period, in contrast to only 0.5 times in the highlands (>2100 m). The majority (80%) of farmers on the Rift floor and 22% in the highlands had experienced at least one crop failure in the 5 years preceding the survey. Table 3 shows the frequency of crop failure at different altitudes. Lack of rain (61% of respondents) and unexpected rain (27%) are the most important determinants of these frequent crop failures. Agriculture in the Rift faces about 76% more crop failures than in the highlands. Crop failures were stated as the cause of food shortages for some households, and a need for food assistance.

‐ **Table 3.** Frequency of reported crop failures in the past 5-year at different elevations.

⸻ ⸻ Practices that farmers are implementing to improve their farming, as reported in the survey, are summarized in Figure 7. About 70% of sample respondents stated that they use high yield crop varieties or early maturing seeds. A large proportion of farmers use fertilizer—–mostly DAP (Diammonium Phosphate) and urea—–and roughly half of them practice crop rotation to improve soil conditions and yield, and to minimize the buildup of weeds, insect pests and plant diseases. The most common crop rotations are between cereal crops. Among irrigators, vegetables and cereals are also commonly rotated. Farmers use terracing and afforestation, conserve seeds and other natural resources, and use irrigation

‐

as key practices to improve productivity. A few wealthy farmers in the highlands use mechanized technology (e.g., renting tractors) to quickly collect crops during harvest. ‐

**Figure 7.** Different measures farmers use to improve productivity.

‐ ‐ ‐ ‐ As a response to rainfall shortages, farmers in the Rift reported increasing planting of drought-resistant early maturing maize varieties (as well as wheat, haricot, and barley), though there is limited supply of such seeds. Farmers sometimes skip planting maize when the rains fail in April and May, and many were shifting away from maize in favor of wheat, barley and teff. In focus groups, many farmers argued that deforestation for firewood and charcoal production were major contributors to environmental or climate changes and expressed concern over the loss of pollinators from overuse of chemical pesticides in the area's floriculture farms. Farmers in several sites indicated they were involved in the protection and planting of new trees (i.e., area closure plans), as part of efforts organized by the regional government, while others mentioned projects to build ponds, enhance water conservation and flood protection, replace chemical fertilizer with manure/compost, and increase crop rotation and diversity to respond to climate stresses. Furthermore, farmers noted that over 350 million trees had been planted in one day (29 July 2019) across the country as part of Ethiopia's Green Legacy Initiative, setting a new record. This effort was conceived to help counter the effects of deforestation and climate change.

‐ ‐ ‐ Lack of access to water, credit or savings, knowledge, information on weather and climate change were identified as major barriers to adaptation to climate and environmental change (Figure S3). Some farmers in the Rift Valley are also reluctant to take available loans when they get access to credit, due to the high uncertainty concerning rainfall, and the associated risk of crop failure and inability to pay back loans. Farmers unable to pay their loans are sometimes forced to sell to livestock or end up mired in bad debts. While market access or transportation was not identified as a major barrier to adaptation, many farmers noted that markets do not function in a way that allows them to maximize earnings. Specifically, more than 50% of farmers sell their crops within 3 months of harvest, which is also the period when prices are lowest. Farmers feel that they have limited bargaining power over prices, which are typically set by local traders. Farmers associations are nonetheless trying to organize to purchase crops from local farmers while the price is low, so as to sell them to the market when prices rise, and allow for increased revenues.

The regression analysis (Table 4) reveals relevant patterns in the determinants of farmer adaptation to climate change. While receipt of information on climate change does not significantly predict adaptation, it does have a positive relationship as would be predicted based on the role of Development Agents in transmitting information. Other studies have similarly found that extension services are positively associated with adaptation, presumably reflecting access to information and other resources that empower farm households to adjust to climate risks [32,50]. Farmer literacy is positive and significant at the *p* < 0.10 level, which is suggestive of a greater ability to receive and process information from external sources, or of correlation with higher socio-economic status and ability to adapt. Farmer productivity (farm revenue per hectare) also has a positive and statistically significant correlation with adaptation at the *p* < 0.10 level, even when controlling for other measures of wealth and economic well-being (though this may also indicate that adaptation leads to greater productivity). While not significant, several other variables have a positive relationship with adaptation, such as experience with crop failure and the number of neighbors. Farmers must take into account a large variety of factors and both individual and community characteristics in order to determine how and when to adapt [58].

**Table 4.** Regression of adaptation behaviors on farmer characteristics.

Standard errors in parentheses; \* *p* < 0.1.

#### **6. Conclusions, Policy Implications, and Future Work**

Climate change impact assessment studies have shown changes in the quantity and variability of rainfall, as well as rising temperatures in Ethiopia—by an average of 1 ◦C since 1960, which is projected to increase by 1.81 ◦C in 2040–2059 [59]. Climate change is threatening rainfed agriculture and rural livelihoods in many parts of the world, especially in highly vulnerable regions and in low-income countries in Africa, such as Ethiopia. The location of this study, the Rift Valley of Ethiopia, is highly vulnerable to climate stressors [5,20,60]. Our study relied on the analysis of primary data on water quality, obtained to evaluate the suitability of different sources for irrigation, as well as survey data pertaining to farmers' perceptions, their adaptation behaviors, and general barriers that challenge adaptation in the MER.

The most salient climate-related shocks in the region are a lack of rainfall and the high variability of rainfall, and crop failures are common. The study also identified several adaptation mechanisms undertaken by households in response to climate-related shocks: the use of new high yield crop varieties, early maturing seeds, increased fertilizer application, and changing crop rotations. They also use terracing, afforestation, seed saving, the conservation of natural resources and irrigation to increase productivity, reduce risk, and control environmental degradation. Irrigation is a common adaptation response to reduced water availability globally, but our assessments of the suitability of irrigation water point to several threats in the MER. First, judging by measures of EC and SAR, most water sources are categorized as having severe or slight to moderate negative implications for

soil permeability and reduced water infiltration. In addition, high bicarbonate and boron concentrations pose water quality threats for cultivation.

It is observed that the role of Extension ("Development") Agents (DAs) is vital to enhancing farmers' knowledge of potential adaptation strategies. This service is currently available in the study basin and is critical for increasing farmers' adaptability in the face of continuing climatic variability. Meanwhile, the main barriers to adaptation include a lack of access to water, credit or savings, lack of appropriate seeds, knowledge and information on weather and climate. Improving access to credit, increased ability to arbitrage sales over seasons with high market prices for local agricultural products, enhancing access and distribution of sufficient seeds (in variety and amounts), early warning on weather during planting and harvesting seasons, and knowledge transfer on effective farming practices (e.g., the use of fertilizers that are suitable and proper amount for a given soil type) are some of the measures that should to be taken to improve farmers' productivity and adaptive capacity to climate change.

Based on the findings of this study, it would be unwise to further develop irrigated agriculture without establishing that water sources being exploited in specific sites are of sufficient quality. Nonetheless, the viability of rain-fed agriculture is clearly threatened by the changing climate, and crop failures are increasingly prevalent in communities, particularly those located in the Rift floor. This situation calls for careful planning of alternative livelihoods strategies for households, by decision-makers who understand these challenges and are willing to design new strategies for adaptation, food security, and enhanced rural well-being.

#### *6.1. Water and Farmers' Survey Data Collection Time and Its Implication in the Study Area*

While the data collection was carried out over a decade ago, it provides useful and rich reference information that can be used as a baseline for examining future changes pertaining to a very data limited region. We use these data to assess the dynamics inherent in farming practices, farmers' perceptions of and adaptation to changing climate and crop productivity, in a context where climatic changes were already being experienced in significant ways. Observing the livelihood of the communities in this region over time, we are confident that no major changes in practices have occurred currently and that agriculture remains unproductive, with the use of antiquated farming practices, very limited or no technology use, and low capacity and resources to invest in improved productivity. Given most groundwater sources in the Rift are not suitable for agriculture due to salinity, which is one of the reasons for the increased pressure on the use of the fresh water sources from Ziway lake and its tributaries. Subsistence agriculture continues to be the main source of income for most farmers in the region. While it may be necessary to monitor specific changes in farmers' behavior and practices over the recent period, we do not anticipate significant changes given that farmers have generally remained in similar conditions marked by a lack of capacity to respond to changing environmental conditions. This includes responses to climate variability and trends, which have increased vulnerability relative to conditions experienced ten years ago.

#### *6.2. Implication of the Study to Current and Future Local Policies*

Our findings on adaptation and barriers to adaptation suggest the need to intensify agricultural productivity by increasing irrigation and technology use in the farming practices in the region. This will improve farmers' food security and income. Considering the farmers' high reliance on unpredictable rainwater for agriculture, policy driven actions that support irrigation use, with particular attention not only to water availability but also to water quality, are critical. Enforcing the regulation of water quality is particularly important to avoid salinization of the soil, which may further exacerbate declining crop productivity in the face of ongoing climate warming. Since freshwater scarcity in the region is prevalent, water use is more determined by its availability or quantity than water quality in the region.

#### *6.3. Future Work*

Our study highlighted water quality issues that could threaten irrigation-based adaptation in the MER, while demonstrating that the region's farmers perceive changes in climate (temperature and rainfall variability) but are limited in their ability to adapt to them. These results can be useful as baseline information for additional in-depth studies on the impacts of surface and ground water irrigation use on soil salinity and agricultural productivity. Further work should be carried out to identify suitable water management strategies that address water scarcity while paying heed to water quality challenges (e.g., high salinity, SAR, and bicarbonate). In addition, the study showed the dependence of farmers on specific seed varieties (early maturing, drought-tolerant, and high yielding crops) whose supply is limited, and whose uptake may have other unintended long-term consequences. High-yield seed varieties that are suitable for saline soil or viable in arid and semi-arid conditions continue to be needed, and soil stability and recovery must be carefully considered in future work.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/cli9060092/s1.

**Author Contributions:** T.R.G.: Design of study concept, conducted the field work, data analysis and interpretation, wrote the manuscript, and funding acquisition. M.A.J.: Design of study concept, critical revision of the manuscript for important intellectual content, and funding acquisition. C.J.P. and P.G.M.: Data analysis and interpretation, critical revision of the manuscript for important intellectual content. D.L.B.: Field work assistance, revision of the manuscript for important intellectual content and funding acquisition. All authors have read and agreed to the published version of the manuscript.

**Funding:** We are grateful to the grant from START's Global Environmental Change Research in Africa under the theme climate change, agriculture, and food security with an emphasis on the sustainability of ecosystem services.

**Institutional Review Board Statement:** The study protocol was approved under Duke University's Institutional Review Board (protocol number A0741).

**Informed Consent Statement:** All respondents provided informed consent prior to participation in the study. The confidentiality and anonymity of survey respondents has been maintained.

**Data Availability Statement:** Survey and water quality data of this study can be provided upon request.

**Acknowledgments:** We would like to thank the farmers who participated in the study for their interest, cooperation and partnership. We would like to thank the Kulumsa Agricultural Research Center (KARC) employees who participated in the agricultural survey interview (Eshetu Lemma and Tarekegn Etana) and for facilitating and guiding us through the field work. We also give special thanks to Skip Kaufmann for his consistent support throughout the project; he helped ensure smooth execution of the project and achievement of its aims.

**Conflicts of Interest:** All authors confirm that they have no financial or non-financial interests.

#### **References**

## *Article* **A Model to Assess Eastern Cottonwood Water Flow Using Adjusted Vapor Pressure Deficit Associated with a Climate Change Impact Application**

**Ying Ouyang 1,\* , Theodor D. Leininger 2 , Heidi Renninger 3 , Emile S. Gardiner 2 and Lisa Samuelson 4**

**Abstract:** Short-rotation woody crops have maintained global prominence as biomass feedstocks for bioenergy, in part due to their fast growth and coppicing ability. However, the water usage efficiency of some woody biomass crops suggests potential adverse hydrological impacts. Monitoring tree water use in large-scale plantations would be very time-consuming and cost-prohibitive because it would typically require the installation and maintenance of sap flux sensors and dataloggers or other instruments. We developed a model to estimate the sap flux of eastern cottonwood (*Populus deltoides*. Bartr. ex Marsh.)) grown in bioenergy plantations. This model is based on adjusted vapor pressure deficit (VPD) using Structural Thinking and Experiential Learning Laboratory with Animation (STELLA) software (Architect Version 1.8.2), and is validated using the sap flux data collected from a 4-year-old eastern cottonwood biomass production plantation. With R 2 values greater than 0.79 and Nash Sutcliffe coefficients greater than 0.69 and *p* values < 0.001, a strong agreement was obtained between measured and predicted diurnal sap flux patterns and annual sap flux cycles. We further validated the model using eastern cottonwood sap flux data from Aiken, South Carolina, USA with a good agreement between method predictions and field measurements. The model was able to predict a typical diurnal pattern, with sap flux density increasing during the day and decreasing at night for a 5-year-old cottonwood plantation. We found that a 10% increase in VPD due to climate change increased the sap flux of eastern cottonwood by about 5%. Our model also forecasted annual sap flux characteristics of measured cycles that increased in the spring, reached a maximum in the summer, and decreased in the fall. The model developed here can be adapted to estimate sap flux of other trees species in a time- and cost-effective manner.

**Keywords:** cottonwood; climate change; sap flux; STELLA; vapor pressure deficit

#### **1. Introduction**

Fossil fuel consumption has been strongly tied to various sources of environmental degradation. Climate change, in particular, is attributed in part to CO2 emissions, and this issue has motivated continuing efforts to identify sustainable alternatives to fossil fuels [1–4]. Various types of biomass have been identified as renewable bioenergy feedstocks of global significance [4–9]. Algae, agronomic crops, grasses, trees, municipal wastes, and other biological materials can be sourced for biomass, and biomass can be converted to solid, liquid, or gaseous fuels to yield energy for industrial, commercial, and domestic uses. The US Department of Energy is promoting an increased usage of biofuels in the United States by targeting the creation of a biomass supply of 1 billion dry tons per year by 2030 [4].

**Citation:** Ouyang, Y.; Leininger, T.D.; Renninger, H.; Gardiner, E.S.; Samuelson, L. A Model to Assess Eastern Cottonwood Water Flow Using Adjusted Vapor Pressure Deficit Associated with a Climate Change Impact Application. *Climate* **2021**, *9*, 22. https://doi.org/10.3390/ cli9020022

Received: 11 December 2020 Accepted: 19 January 2021 Published: 23 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Tree and shrub cultivation in short-rotation (2–15 years) provides for high-yielding bioenergy feedstock production. Over the past several decades, researchers have refined improvements in short-rotation tree and shrub species such as eastern cottonwood (*Populus deltoides* Bartr. ex Marsh.), willow (*Salix spp*. L.), and eucalyptus (*E. globulus* Labill.) by selecting for rapid growth, pest tolerance, and site suitability [10,11]. Eastern cottonwood is a fast-growing poplar and one of the largest native North American broadleaf trees [12]. Eucalyptus is among the fastest growing hardwood plantation species utilized around the world. Its suitability to plantation culture is driving current research and development towards application to large-scale biomass production in the southern United States [11,13]. Shrub willows have been selected for short-rotation intensive culture in the northeastern, north–central, and mid-Atlantic regions of the United States [14].

There are, however, concerns about the water use of trees in short-rotation bioenergy systems, especially in regard to the potential adverse environmental impacts of excessive water use [15–19]. Scott and Lesch [20] reported a 90 to 100% reduction in stream flow for eucalyptus and a 40 to 60% reduction in stream flow for pines (*Pinus spp.* L.) in the first 8 years or so after planting. Morris et al. [21] studied water use in the eucalyptus plantations of southern China and noted that annual water use by these plantations is about 550 mm, but potential annual water use could reach 865 mm where soil water was not limiting. Flanagan et al. [19] used the eddy covariance measurement of evapotranspiration (ET) to estimate water use in a riparian cottonwood ecosystem in Alberta, Canada. They found a comparatively high water use driven by relatively warm summer temperatures and high vapor pressure deficit ((VPD) which is the difference between the amount of moisture in the air and the maximum (saturated) amount of moisture the air can hold) and supported by sufficient availability of alluvial groundwater. In regard to VPD, Weih [22] reported that fast-growing hybrid willows are generally more sensitive to water stress and have lower water use efficiency than slower-growing wild clones. Though the studies cited above provide valuable insights into water use by trees in short-rotation bioenergy systems, the long-term, continuous, and real-time monitoring of water use by trees in biomass production plantations has rarely been conducted.

Sap flux measurement is currently the most widely used method for assessing individual tree water use. Its advantages include a well-established technology, the ability to continuously monitor trees for lengthy time periods, and the fact that instrumentation does not modify the ambient environment [23,24]. However, for commercially sized biomass production plantations, it would be very time-consuming and costly to install and maintain a sufficient number of sensors and dataloggers to accurately and continuously assess tree water use with real-time resolution. Therefore, there is a need to develop practical and reliable methods for estimating tree water use in intensively managed short-rotation woody crop systems.

This problem may be addressed with an understanding of the relationships between sap flux and its driving environmental variables [24–29]. Hogg et al. [25] compared the sap flux and eddy flux of water vapor in a boreal deciduous forest and found that the sap flux response lagged behind eddy flux measurements by about 1 hour diurnally—presumably because of water storage within the trees. Nadezhdina [26], who studied sap flux in apple trees (*Malus domestica Borkh.*) over one growing season, found correlations to the diurnal patterns of leaf water potential, air temperature, relative humidity, and solar radiation. Wilson et al. [27] compared forest evapotranspiration with sap flux, VPD, soil water budget, and eddy covariance. Reasonably good agreement was observed between sap flux and VPD as well as sap flux and eddy covariance based on annual data. Liu et al. [29] studied the sap flux of sawtooth oak (*Quercus acutissima* Carruth) and china-fir (*Cunninghamia lanceolate* (LAMB) Hook.) relative to photosynthetically active radiation (PAR), air temperature, relative humidity, VPD, precipitation, and soil water content. Significant linear correlations were found between sap flux and PAR, air temperature, and VPD, but not soil water content. These studies establish mechanistic relationships between tree water flux and certain environmental variables based on the response of stomata to environmental drivers,

but most of the reported relationships were obtained using annual data that provided low correlation coefficients. Diurnal sap flux patterns are typically not distinctly parallel to the diurnal patterns of environmental variables. Therefore, it is difficult and inaccurate to simply predict diurnal sap flux directly from environmental variables.

Climate change is a natural phenomenon. However, anthropogenic activities such as fossil fuel burning, industrial pollution, deforestation, and population growth have greatly accelerated greenhouse gas emissions and resulted in extremely abnormal climate change patterns [30]. Climate change over the last several decades has been linked to atmospheric water vapor content increase, precipitation pattern shifts, snow cover reduction and ice melt, and surficial hydrological cycle changes [30]. Lasch et al. [31] investigated the impacts of short-rotation coppice plantation with aspen (*P. tremula* L.) in Eastern Germany under changing climate conditions. These authors found that the aspen plantations contributed to regional CO2 mitigation and carbon sequestration, but had negative on the regional water budget due to climate change. Griffiths et al. [32] reviewed the environmental effects on short-rotation woody crop production. They cited long-term water use, C dynamics, and soil quality studies as being needed to evaluate the potential effects of climate change. Currently, our understanding of climate change effects upon the sap flux of eastern cottonwood in short-rotation plantations is basically unknown.

The goal of this study was to develop a novel model for estimating diurnal and annual sap flux of trees in short-rotation biomass production plantations associated with climate change impacts. To accomplish this goal, we based our model on adjusted VPD using Structural Thinking and Experiential Learning Laboratory with Animation (STELLA) software, and used eastern cottonwood as the plantation tree species. We investigated the relationship between adjusted VPD and eastern cottonwood sap flux to address the following objectives: (1) establish a model using diurnal patterns and annual cycles of adjusted VPD to predict those of measured eastern cottonwood sap flux; (2) validate the model using an independent sap flux dataset collected from sensor measurements; and (3) apply the model to predict eastern cottonwood sap flux in the absence of measured sap flux data under a changing climate. It should be noted that we approached the problem of predicting sap flux with the variable VPD because VPD data are easy to obtain from local weather stations, and sap flux tends to follow VPD more closely than other environmental variables such as air temperature and PAR. The model developed here provides a new paradigm for estimating the sap flux of trees in short-rotation bioenergy plantations in a time- and cost-effective manner.

#### *1.1. Study Site*

This work is based on sap flux measurements collected at a biomass production plantation located near Hollandale, Washington County, Mississippi, US (see Figure 1). The plantation was established in 2012 to study eastern cottonwood and black willow (*S. nigra* L.) biomass production under several planting densities and harvest regimens.

The plantation is situated on a site with Sharkey clay soil that is poorly drained, and the site was previously used for rice (*Oryza sativa*) production. For this study, we instrumented six representative eastern cottonwood trees (each 3 years old) with a planting density of 1 m × 1 m and a plot area of 1.01 hectare, with 2 cm long laboratory made heat dissipation type sensors [33] to collect sap flux measurements from September 2015 to September 2016. One constructed sensor pair consists of a reference and heated sensor, each containing thermocouples within hypodermic needles. Each sensor was installed radially in the xylem of each tree at a height of 120 cm, where the diameter at breast height (DBH) was measured. The sensors were about 10 cm apart vertically, and the immediate area of the stem was covered with insulated aluminum shielding that prevented solar radiation from impacting readings. All sensors were connected to a CR1000 datalogger and AM 16/32/B multiplexer (Campbell Scientific Inc., Logan UT, USA) and the instrumentation was powered by rechargeable marine deep cycle batteries connected to a solar panel (Figure 1b). The system was programmed to measure voltage differences between sensors

corresponding to thermocouple temperature differences every 30 s, and data were averaged every 30 min. As water flows through the xylem, it cools the heated sensor relative to the reference sensor and, therefore, temperature differences between sensors can be used to infer sap flux rates. Collected data were transmitted to an office computer through a wireless carrier (Verizon Inc.). This real-time system enabled periodic monitoring for malfunctions and quick troubleshooting. Meanwhile, time series VPD data with the same collection interval (i.e., every 30 min) were downloaded from an on-site weather station (https://www.hobolink.com). Missing data from the on-site station were replaced with data from the nearest weather station, which was located about 40 km north of the study site (https://www.wunderground.com/history /airport/KGLH/). The VPD data were calculated using the following formula [34]:

$$D\_{vp} = \left[1 - \left(\frac{RH}{100}\right)\right] \ast 610.7 \ast 10^{\frac{7.87\_{air}}{237.3 + I\_{air}}}\tag{1}$$

where *Dvp* is the vapor pressure deficit (Pa), *RH* is the relative humidity (%), and *Tair* is the air temperature (oC). Differences in sap flux sensor temperature were converted to sap flux densities (g/m2/s) with the widely accepted empirical equation developed by Granier [33] using the software program Baseliner version 4. Data were scaled to the tree level (g/s) with DBH measurements on the assumption that the entire cross-sectional area was conductive sapwood.

**Figure 1.** A short-rotation woody crop plantation used in this study near Hollandale, Mississippi, USA (**a**) and real-time monitoring of sap flux from eastern cottonwood using the sensor–datalogger system (**b**).

#### *1.2. Adjusted VPD and its Correlation with Sap Flux*

௩ = 1 − ൬

100൰൨ ∗ 610.7 ∗ 10 .ହ்ೌ

ଶଷ.ଷା்ೌ

The measured hourly sap flux and VPD over the study period illustrated that the highest diurnal peaks for these two variables were not aligned (Figure 2a). That is, the highest peaks for hourly sap flux occurred mostly in the summer, whereas the highest peaks for hourly VPD were observed in the early fall. Analysis of hourly measured sap flux and VPD over the one-year study period (Figure 2b) revealed a very weak relationship between these variables (R2 = 0.263). Further examination of measured sap flux and VPD at diurnal scales in spring, summer, and fall (there are no leaves in winter) revealed the following three distinct patterns (Figure 3): (1) sap flux peaked around 10:00 a.m. under normal conditions; (2) peak VPD typically lagged about 3 hours behind peak sap flux; and (3) sap flux maintained a constant value close to zero from 9:00 p.m. at night to 6:00 a.m. the next morning, while VPD remained above zero during the same period. Based on these observations, unadjusted diurnal VPD would be unsuitable for predicting diurnal sap flux. However, VPD could be adjusted to develop a stronger correlation with sap flux.

**Figure 2.** Time series plots of the measured sap flux and the measured vapor pressure deficit (VPD) (**a**) and their correlation (**b**) for eastern cottonwood in a biomass production plantation near Hollandale, MS, USA.

We employed the following eight steps (Figure 4) to adjust measured VPD and develop a mathematical relationship between sap flux and adjusted VPD.

estimated from our experimental data, and this value can be changed for different tree species or sites.

**Figure 3.** Diurnal sap flux and VPD as measured on selected spring, summer, and fall dates for eastern cottonwood in a biomass production plantation near Hollandale, Mississippi, USA.

4. A correction factor was imposed to adjust the annual cycle of VPD to match the annual sap flux cycle. This correction factor (*F*) was calculated with Equation (2) below.

$$F = 1.1 \ast EXP\left((-0.3) \ast \left(\frac{t-\beta}{\alpha}\right)^2\right) \tag{2}$$

where *t* is the hour of year, *β* is a time constant that specifies the hour of year when maximum sap flux occurs (5100 h in this example), and *α* is a constant characterizing the shape of the annual sap flux cycle (1500 in this example). This constant was obtained from measured annual eastern cottonwood sap flux.

$$\begin{array}{rcl} \text{Q}\_{sapflow} = & 5E-08\*VPD\_{adjusted}^3 - 0.0006\*VPD\_{adjusted}^2\\ &+ 2.5167\*VPD\_{adjusted} - 16.45 \end{array} \tag{3}$$

*α*

*β*

where *Q* is the sap flux (g/m2/h) and *VPDadjusted* is the adjusted VPD (Pa). A detailed discussion of Equation (3) is provided in the Results and Discussion Section. ௦௪ = 5 − 08 ∗ ௗ௨௦௧ௗ ଷ − 0.0006 ∗ ௗ௨௦௧ௗ + 2.5167 ∗ ௗ௨௦௧ௗ − 16.45

8. Because eastern cottonwood sap flux is affected by age (which is correlated to DBH), we accounted for the age effect with Equation (4):

$$F\_{age} = 0.038e^{0.0431DBH} \tag{4}$$

ଶ

where *Fage* is the age factor and DBH (cm) increases as time elapses. Equation (4) was developed from data reported by Schaeffer et al. [35], who measured the sap flux and DHB of eastern cottonwood growing in the San Pedro River system of southeastern Arizona. Combining Equations (3) and (4) yields the following equation: = 0.038.ସଷଵு

$$\mathcal{Q}\_{sapflow}^{adjusted} = \left(1 + F\_{\text{age}}\right) \* \mathcal{Q}\_{sapflow} \tag{5}$$

where *Q adjusted sap f low* is the age-adjusted sap flux (g/m2/h). Equation (5) is used to predict age-adjusted eastern cottonwood sap flux in this study. ௦௪ ௗ௨௦௧ௗ

**Figure 4.** A flowchart showing the procedures used to develop a Structural Thinking and Experiential Learning Laboratory with Animation (STELLA)-based model for predicting sap flux using the adjusted VPD.

**Figure 5.** Measured sap flux and adjusted VPD for selected dates in spring, summer, and fall for eastern cottonwood in a biomass production plantation near Hollandale, MS, USA.

#### *1.3. STELLA Model*

The eight steps listed above were applied to the development of a STELLA model. STELLA is a software package used to develop system dynamic models by creating a pictorial diagram of a system and assigning appropriate values and functions to it (http://www.iseesystems.com). System dynamic models created with STELLA have been widely used [36–39]. Figure 6 provides the STELLA modeling map, which illustrates the conditions and equations (in the dashed boxes) established for this study. The STELLA model (with an hourly time step) we created is further described below. Referring to Figure 6, data for Step 1 is stored as input data in the converter (or a circle) labeled "Measured VPD with shifting 3 hours ahead." Steps 2, 3, and 5 are completed in the converter labeled "Adjusted VPD" according to the conditions in the dashed box.

**Figure 6.** A STELLA modeling map for estimating the sap flux of eastern cottonwood from adjusted VPD.

The routine first determines if time of year (in hours) is less than 2500 h or greater than 7950 h (representing the winter period). If either condition is true, adjusted VPD is set to zero. Then, if time of a day is <7 h or >21 h (representing the nocturnal period), adjusted VPD is set to the minimum value for a site—the value for our study site is 33 Pa. If the time of day is 10 am, the VPD, after having been shifted 3 hours, was adjusted by multiplying it by 2.4 (as explained in Step 3). Step 4 and Equation (1) are processed in the converter labeled "Annual correction factor." It should be noted that we initiated simulation at the beginning of winter (i.e., 12 December 2015) and terminated simulation at the end of the following fall (i.e., 11 December 2016), but start and end times can be changed as appropriate. It should also be noted that our field experiment started in September 2015 and the collection of sap flux data began in December 2015. After adjusted VPD was calculated, predicted sap flux was estimated with Equation (5) in the converter labeled "Predicted sap flow."

#### **2. Results and Discussion**

#### *2.1. Measured Sap Flux, Adjusted VPD, and Model Validation*

The relationship between measured sap flux and adjusted VPD for the one-year study period is illustrated in Figure 7. The values of R2 = 0.823 and *p* < 0.001 (Figure 7a) indicate

that sap flux can be reliably predicted from adjusted VPD. This was further confirmed in Figure 7b, which shows that the diurnal peaks of these two variables shared a relatively similar annual cycle. Though the results above show that measured sap flux correlated well with adjusted VPD, they do not confirm if sap flux predicted by the STELLA model accurately matches measured sap flux. Therefore, model validation is necessary prior to its application. It should be pointed out that sap flux measurements collected from three of the six eastern cottonwood trees instrumented in this study were used to develop Equation (3) and perform initial validation of the model, while sap flux data from the three remaining trees were used as an independent dataset for additional validation of the model. The initial validation of predicted sap flux against measured sap flux from the first three trees is provided in Figure 8. The validation was measured using statistical parameters such as R2 , the Nash Sutcliffe coefficient (NSE), and normalized root mean square error (nRMSE). The nRMSE is calculated as follows [40]:

$$\text{nRMSE} = \frac{1}{\mathcal{O}\_{ave}} \left( \sqrt{\frac{\sum\_{i=1}^{n} \left( \mathcal{O}\_{i} - \mathcal{S}\_{i} \right)^{2}}{n}} \right) \tag{6}$$

where *Oi* is the field observation, *Si* is the model prediction, *Oave* is the average number of field observations, and *n* is the total number of field observations.

**Figure 7.** The relationship between measured sap flux and adjusted VPD (**a**) and the time series plots of measured sap flux and the adjusted VPD (**b**) for eastern cottonwood in a biomass production plantation near Hollandale, Mississippi, USA.

**Figure 8.** A comparison of the measured and predicted sap flux (**a**) and their time series plots (**b**) during the model validation process for eastern cottonwood in a biomass production plantation near Hollandale, Mississippi, USA.

The NSE is given as [41]

$$\text{NSE} = 1 - \frac{\sum\_{i=1}^{n} (\mathbf{O}\_i - \mathbf{S}\_i)^2}{\sum\_{i=1}^{n} \left(\mathbf{O}\_i - \overline{\mathbf{O}}\right)^2} \tag{7}$$

The NSE ranges from -∞ to 1, with the values of 1 representing a perfect fit,>0.75 representing a very good fit, between 0.36 and 0.75 representing a reasonable fit, and <0.36 representing an unsatisfied fit to the model. With R2 = 0.815, NSE = 0.697, nRMSE = 1.202 (g/m2/h), and *p* value < 0.001 (Figure 8a), we suggest that the model developed in this study performed well in predicting eastern cottonwood sap flux. This conclusion was further supported by data presented in Figure 8b, which illustrate diurnal peaks and seasonal trends of predicted sap flux to coincide with those of measured sap flux. For additional validation of the model, we compared predicted sap flux to the measured sap flux from trees in the independent dataset. This second validation confirmed that the developed model predicted sap flux reasonably similar to the measured sap flux based on the relatively high R2 (0.797) and NSE (0.677) as well as the low nRMSE (2.002) and *p* value (<0.001) and its conformation to measured sap flux peaks (Figure 9a).

**Figure 9.** A comparison of the measured and predicted sap flux for eastern cottonwood in a biomass production plantation near Hollandale, MS, USA (**a**) and an experimental plantation near Aiken, SC, USA (**b**).

To develop users' confidence, we further validated the model by comparing the measured sap flux from a different study site to the predicted sap flux from our model. The measured sap flux data were obtained from Samuelson et al. [42]. These authors conducted an experiment to estimate the influence of irrigation and fertilization on transpiration and hydraulic properties of the 3-year-old eastern cottonwood. Their study used marsh clones at a 24.4 ha experimental plantation on the U.S. Department of Energy Savannah River Site located near Aiken, SC. In this validation process, we used the measured sap flux data from Block 2 of their experiment along with the VPD data calculated from local weather variables (i.e., air temperature and relative humidity). A closer look at the measured sap flux data revealed that the maximum sap flux from Samuelson et al. [42] was 5716.8 kg/m2/s in Aiken, SC but was 586.8 kg/m2/s in our experiment in Hollandale, MS. The former was about 10 times greater than the latter. Since the sap flux, weather conditions, and eastern cottonwood clones in Aiken, SC are different from those in Hollandale, MS, we have modified Equation (3) for the Aiken study site as follows:

$$Q\_{sapflow} = -1E - 06 \ast VPD\_{adjusted}^2 + 1.119 \ast \ VPD\_{adjusted} + 12789 \tag{8}$$

It should be noted that other input values from Equations (2), (4), and (5) were not changed during this validation process. A comparison of sap flux densities between model prediction and experimental measurement is given in Figure 9b. The model predicted most of sap flux densities very well, except for several unusually high measured sap flux densities. With R2 = 0.586, NSE = 0.580, nRMSE = 0.897, and *p* < 0.01, we concluded that the model predicted sap flux density reasonably well.

#### *2.2. Model application*

The model was applied to predict eastern cottonwood sap flux at the same study site for a simulation period beginning on 30 November 2016 and ending on 30 October 2017. Input data were the same as for model validation except for measured DBH and VPD. Average DBH was changed to 7.2 cm for this simulation period. This number was based on the measured growth in the plantation. Data to calculate VPD for the simulation period were downloaded from the on-site weather station.

The diurnal pattern of predicted sap flux by eastern cottonwood between 22 April and 28 April 2017 is shown in Figure 10a. This projection conforms to the expectations for the diurnal sap flux volume and the pattern of eastern cottonwood. Relatively low sap flux on the first day (0 to 24 h) was a result of low VPD, normally observed on cloudy or rainy days. When higher VPD prevailed, such as a VPD between 72 and 96 h, we observed predicted sap flux in a pattern similar to previous measurements of sap flux. Maximum sap flux during this time period was 3577 g/m2/h, or 16 g/h for the 7.2 cm DBH tree. This rate was within the range reported by Schaeffer et al. (2000), who found that eastern cottonwood sap flux in the spring averaged about 20 g/h. Results indicated that the model provided for reasonable prediction of diurnal sap flux by eastern cottonwood grown in a biomass plantation.

**Figure 10.** Predicted diurnal (**a**) and annual (**b**) sap flux of eastern cottonwood in a biomass production plantation near Hollandale, Mississippi, USA.

An annual cycle of predicted eastern cottonwood sap flux is presented in Figure 10b. The simulation estimated total annual sap flux at 6,369,309 g/m2/y. Converting this value to be representative of a tree with a DBH of 7.2 cm yields an annual water use of 25,919.52 g/year.

We assumed soil water content did not limit eastern cottonwood sap flux in this study. Soil water dynamics and hydrological processes would need to be included in the STELLA

model if this factor were to be considered important. While we developed this model by studying eastern cottonwood sap flux, it could be adapted to estimate the sap flux of other tree species in short-rotation biomass production plantations. This would be done by revising values outlined in the aforementioned eight steps (Section 1.2) to reflect the tree species of interest and the site.

#### *2.3. Impact of Climate Change on Sap Flux*

To estimate the impact of climate change on sap flux, a simulation scenario was developed by increasing the VPD by 10% from the original VPD. Up to date, no effort has been devoted to estimating the impact of climate change on VPD and sap flux in our study area. The 10% increase in VPD was considered as an extreme case based on the study reported by Yuan et al. [43]. These authors studied the global VPD and reported that the mean VPD in the growing season from 2011 to 2015 was 11.26% higher than that from 1982 to 1986.

Comparisons of the predicted daily and overall average hourly sap flux of eastern cottonwood between the original VPD and the 10% increased VPD are shown in Figure 11. In general, an increase in VPD resulted in an increase in sap flux (Figure 11a). The results indicate that as the VPD increased due to the air temperature increase and the precipitation decrease under the changing climate, the sap flux of eastern cottonwood increased. Overall, the hourly average sap flux for the simulation period from 6 November 2016 to 1 December 2017 was 793 g/m2/h with the original VPD but was 838 g/m2/h with the 10% increased VPD (Figure 11b). The latter was 5.67% greater than the former.

**Figure 11.** Comparisons of the predicted daily (**a**) and overall average hourly (**b**) sap flux of eastern cottonwood between the original VPD and the 10% increased VPD.

#### **3. Conclusions**

This study provides a model to estimate the sap flux of eastern cottonwood grown in short-rotation biomass production plantations. Sap flux was predicted from adjusted VPD, as modeled with STELLA software. The model was validated using two different sap flux datasets collected with sensors between 1 December 2015 and 30 November 2016. Very good agreement was found between predicted and measured sap flux. More specifically, predicted diurnal patterns and annual sap flux cycles matched those of measured sap flux.

The model was then applied to predict eastern cottonwood sap flux for a year-long simulation period (30 November 2016 to 30 October 2017). The model appeared to predict reasonable diurnal patterns and annual cycles of eastern cottonwood sap flux. Overall, a 10% increase in VPD increased the sap flux of eastern cottonwood by about 5%. The model provides a new approach for estimating the sap flux of cottonwood in short-rotation bioenergy plantations in a time-saving and cost-effective manner.

The major limitations of the model were that soil water availability is not a limiting factor for sap flux and the very high sap flux peaks were not easily predicted. Therefore, further study is warranted to include soil water dynamics, hydrological processes, and high sap flux peak accommodations into the STELLA model. Additionally, the response of VPD to sap flux is tree species dependent. Although the model was developed for eastern cottonwood, it can be adapted to estimate the sap flux of other tree species with moderate effort.

**Author Contributions:** Conceptualization, Y.O.; methodology, Y.O. and H.R.; software, Y.O.; validation, Y.O., T.D.L., H.R., E.S.G., L.S.; formal analysis, Y.O.; resources, Y.O. and L.S.; data curation, Y.O. and L.S.; writing—original draft preparation, Y.O.; writing—review and editing, Y.O., T.D.L., H.R., E.S.G., and L.S.; visualization, Y.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data available upon request.

**Acknowledgments:** Matt Moran for help in collecting the sap flux data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

#### **References**

*Article*

## **Climate Drivers and Sources of Sediment and Organic Matter Fluxes in Intermittent Rivers and Ephemeral Streams (IRES) of a Subtropical Watershed, USA**

#### **Janet Dewey 1, \*, Je**ff **Hatten 2 , Byoungkoo Choi 3 , Clay Mangum 4 and Ying Ouyang 5**

Received: 28 August 2020; Accepted: 10 October 2020; Published: 16 October 2020

**Abstract:** Climate-driven hydrological models rarely incorporate intermittent rivers and ephemeral streams (IRES) due to monitoring difficulties and their perceived minor effect on river networks. Worldwide, IRES represent approximately 50% of river networks and up to 60% of annual flow and are recognized as conduits and processors of organic matter (OM). Climate induced changes in precipitation and discharge (Q) may impact OM fluxes from IRES. We assessed storm-driven source and flux of total suspended solids (TSS) and OM from small IRES in Mississippi, USA. We used linear Pearson correlations to evaluate relationships between water and storm characteristics (e.g., discharge). Stepwise regression was used to predict change in flux. Dissolved OM was derived from saturated flow through soil whereas particulate OM was derived from channel extension during storms. A power log relationship between Q and materials flux indicated that Q was the driver for flux. A 5% increase in Q within IRES may result in flux increase of 2% TSS and 1.7–2.8% OM. Climate change projections of increased storm intensity over a shorter water year will increase channel extension and soil water transfer resulting in higher material flux to downstream reaches. Climate-driven hydrological models of OM flux should incorporate IRES.

**Keywords:** climate change; IRES; OM; DOC; POC; ephemeral stream; event sampling; headwaters

## **1. Introduction**

Intermittent rivers and ephemeral streams (IRES) have long been perceived as peripheral in their effect on river networks because they only flow in response to rainfall [1,2]. It has been estimated that headwater IRES account for 60% of totalmean annual flow (for streams and rivers in the northeastern US) [3] and approximately 50% of the world's total river networks [4]. Intermittent rivers and ephemeral streams drain and connect with much of the Earth's critical zone, so they have considerable potential to release, transport, transform, or sequester organic matter (OM) (e.g., [5,6]); this, in turn, can be released into the atmosphere in the form of CO2 [7,8]. Approximately 36% of the CO2 outgassed by river networks can be attributed into headwater streams [8]. Climate change can affect the timing, duration, frequency, and severity of storm events—all of which have the potential to control sediment and OM exports from IRES [9].

Organic matter is a transport vehicle for nutrients, carbon, and pollutants, all of which regulate or impact biological processes [10,11]. Production, transport, degradation, preservation, and sorption of OM are influenced by complex biological, chemical, and hydrological interactions [11–13]. Increases in temperature, precipitation, discharge, and groundwater connectivity associated with climate change may accelerate chemical and physical processes governing OM dynamics [14–19]. Therefore, understanding the source–transport–transform–sink behavior of OM in ephemeral streams (ES) is critical to predicting how OM will affect downstream processes and water quality.

Ephemeral streams function as collection sites for organic inputs and typically flow in response to storm events as infiltration excess [20] or saturated excess [21] overland flow with little or no connectivity to the water table (e.g., Table 1 in Reference [22]). Episodic connectivity during precipitation and flow pulses permits transport of OM, including dissolved and particulate C, to downstream perennial reaches. Storms are important drivers for material flux in headwaters; approximately 90% of dissolved and particulate material in headwater streams is exported during storm events [23]. This flux occurs disproportionately during high flow and is often a power function of the discharge rate [24–26]. Ephemeral streams have been poorly studied in comparison to their intermittent and perennial counterparts with regard to nutrient and energy fluxes due to the fact of their measurement difficulties and rapid discharge (i.e., "flashy behavior").

Recent research suggests that in-channel processing of OC within small inland watersheds is a regionally important source of atmospheric CO2 [7,8,27]. Frequent wetting and drying cycles that occur in ephemeral streams provide "hot moments" [28] for dissolved organic matter (DOM) production, transformation, and transport. A number of studies have shown that oscillation between wet and dry stream conditions facilitates dissolved organic carbon (DOC) processing (e.g., [29–31]) and CO2 evasion [32]; therefore, even very small ephemeral streams may be key components of C cycling in riverine systems [33].

Current climate change scenarios predict an increase in the frequency of droughts and high magnitude storm events in the southeastern United States even if mean annual precipitation (MAP) does not change [34]. Predicted changes in frequency [35] and severity of precipitation [36] may increase effective sediment and OC transport processes [37–39]. Equally, the dry component of the wet–dry cycle of ephemeral streams will increase both spatially and temporally as temperatures rise [32] thereby increasing DOM production [40]. Much of the research on streams that exhibit wet–dry cycles has focused on arid and semi-arid regions [41] and neglected ephemeral streams in humid climates. Therefore, understanding the behavior of C in temporary streams of humid ephemeral source areas is important for understanding the contributions of those source areas to C dynamics in general and under accelerated transport efficiency as predicted by some climate change scenarios. Our research was guided by two questions: (1) what are the sources of sediment and OM in IRES and (2) what are the controls on the flux of sediment and OM from very small (<4 ha) IRES to downstream reaches? We used chemical characteristics to determine the source and behavior of OM during storm events within small ephemeral catchments and tracked the behavior of potential exports during storm events. We discuss the implications of our study in the context of climate change.

#### **2. Materials and Methods**

#### *2.1. Study Design*

We studied a 32 ha watershed (33◦30′49.00′′ N, 89◦25′40.00′′ W) within the Upper Gulf Coastal Plain in Webster County, Mississippi. The site has a humid subtropical climate with a mean temperature ranging from 7 ◦C in winter to 26 ◦C in summer (US National Weather Service Station #222896, Eupora, MS, USA). The 30 year mean annual precipitation is 1451 mm, much of which falls during winter and spring. Annual discharge is 240 × 103 m3 year−1 [42]. Soils are well- to moderately well-drained Sweatman series: US Department of Agriculture (USDA) taxonomy = fine, mixed, semiactive, thermic Typic Hapludult; World Reference Base (WRB) taxonomy = Profondic Alisol. Soils are high in clay with A-horizons of either loam or silt loam; pH ranges from 4.5 to 5.5. Slopes are steep, soils are erodible, and dominant land use is woodland [43]. Forest vegetation is characteristic of the southeastern Mixed Forest Province and consists of even-aged loblolly pine (*Pinus taeda* Linnaeus) with a smaller component of mixed hardwoods; a detailed description is provided in Reference [44]. The site has been in pine timber production for 30 years.

Four small ephemeral catchments, nested within a perennial stream watershed (Figure 1), were monitored for dissolved and particulate contributions to downstream water. We used the following field criteria [22] to classify streams: (1) weakly defined channels; (2) interrupted flow, most of which was in response to storms; (3) a water table below the channel surface for most of the year; (4) lacking in aquatic insects; (5) lacking in obvious material transport; (6) scoured beds; and (7) intermittent organic buildup. Ephemeral streams were mapped as depressions by the US Geological survey (1983). We classified the streams as ephemeral based on the aforementioned criteria, recognizing that the boundary between ephemeral and intermittent classification is approximate. Our major focus in this study was on the ephemeral streams (ES) of the IRES.

**Figure 1.** Schematic project layout depicting monitored ephemeral catchments in the study watershed. Numbers 1–4 correspond to the named ephemeral 1–4 watersheds described in this study. Topographic base excerpt are from the Little Sand Creek, Mississippi quadrangle [45].

For monitoring purposes, we defined the upstream limits of the ES segments by the upslope limit of channel development; downstream limits were defined by evidence of a clearly defined channel, over bank deposits, and seasonal streamflow [44]. All monitored ES lacked evidence of a persistent intermittent segment: the ES conducted surface water directly to the perennial stream and flowed in response to storm events. However, storm events in this region occurred frequently enough that steady-state flow (hereafter referred to as base flow) was periodically observed at the downslope limits of the channels. The ES watershed size ranged from 1.8 to 3.8 ha; the remainder (20.6 ha) comprised a mix of uncut and recently cut forest (Table 1) which had recovered for two years prior to monitoring for this study. Soil disturbance and forest cover were not a focus of this paper but are described as a frame of reference for interpreting OM dynamics in the watersheds. Ephemeral channels differed in their forest cover (Table 1 and described in References [44,46]) and were chosen to represent a continuum of soil disturbance and forest cover that may influence OM dynamics.

**Table 1.** Watershed characteristics and sampling record.

Total Q was calculated from all Q measurements across the entire length of the record. \* Reference site at which no harvesting has occurred since 1978. † Riparian buffer in which logging debris was not permitted in the ephemeral drainage channel. ‡ Riparian buffer in which logging debris was permitted in the ephemeral drainage channel. § Clearcut with all merchantable timber removed and no riparian buffer.

#### *2.2. Field Sampling*

Ephemeral streams (Figure 2) were equipped with stations to monitor stream flow and collect event samples over 17 months (17 February 2010 to 13 July 2011). A 1.8 m long, 254 mm diameter PVC pipe was installed at the outlet of each ephemeral watershed to constrain flow within a measurable cross-sectional area. Water depth, velocity, and discharge within pipes were measured by area velocity sensors controlled by flow loggers (ISCO Inc., Lincoln, NE, USA) programmed to record at 15 min intervals. Discrete water samplers (ISCO Inc., Lincoln, NE, USA) were linked to flow loggers and programmed to start collecting samples on the rising limb of each event. Event sampling was triggered by a 20 mm increase in flow depth in the pipe over a 15 min interval. Discrete samplers were programmed to cease sampling as velocity and/or depth reverted to pre-event levels. Flow was variable among streams and events, so the number and timing of sampled storm events differed among streams. Steady-state base-flow samples were hand-collected from pipe outflow when possible.

**Figure 2. Left**: Gauged reference stream with a 254 mm diameter pipe, ISCO sampler, and flowmeter. **Right**: Ephemeral stream flowing during storm event.

The perennial stream monitoring station was located 200–300 m downstream from ephemeral monitoring stations. Perennial stream stage was measured with a pressure transducer (In-Situ Level TROLL 300) inside a stilling well (91.8 cm tall, 152 mm diameter). Velocity was measured using a portable flowmeter (Hach Company, Loveland, CO, USA). Discharge was determined using a stage-discharge rating curve, developed from 15 measurements of stream profile, stage, and velocity measurements across a range of flows. The perennial stream was equipped with an area velocity sensor and a discrete sampler which was programmed to initiate sampling in the event of increased depth and/or flow. Base-flow samples were hand-collected during field visits when discharge was at steady state.

Groundwater, precipitation, canopy throughfall, and soil pore-water samples were used to characterize potential source contributions to streamflow. Groundwater was sampled in ephemeral watersheds via a grid of wells arranged in five transects perpendicular to the stream channel; each transect contained five wells spaced at 5 m intervals. Wells were made of 5 cm inner diameter PVC pipe and screened the entire subsurface interval. Wells penetrated to a depth of 2.5 m or until refusal. Well water was sampled monthly using a bailer with a check ball after purging wells. Depth of sampling for the groundwater source varied depending upon well and time period [44,46]. Regional precipitation data (15 min interval) for the period 10 February 2010 to 30 May 2010 were acquired from the Eupora 2E, MS, US (33.56◦ N, 89.24◦ W) National Oceanographic and Atmospheric Administration (NOAA) weather station [47]. Onsite precipitation data were collected for the period 30 May 2010 through 13 July 2011 using a tipping bucket rain gauge (ISCO Inc., Lincoln, NE, USA) installed in an open area of the watershed. Precipitation and canopy throughfall were sampled using acid-washed plastic 18.9 L buckets with nylon window screens over the top to prevent contamination by fallen debris. Zero-tension lysimeters collected soil pore water four times during the spring of 2011 at three depths within the soil profile (O, A−10 cm, and A−20 cm). O, A, B, and C soil horizons were sampled from soil pits at representative locations (Figure 1).

#### *2.3. Sample Analysis*

We examined the source and flux of water and OM. The flux of OM from IRES is partially controlled by the concentration of the material in the water, which is strongly controlled by source. To understand how climate change will affect the flux of the materials, we must understand the source and mobilization processes of OM.

Water samples were fractionated in the lab into dissolved and particulate components. Samples were filtered to collect particulate matter (PM) through 0.7 µm glass fiber filters (GFFs) that were pre-combusted for 4 h at 350 ◦C to remove organics. The PM-laden filters were oven-dried overnight at 60 ◦C and weighed for total suspended solid (TSS) concentration. Filtrate was collected and frozen for analysis of dissolved constituents. Soils and suspended sediment were analyzed for total C and N. Soils were dried overnight at 70 ◦C, ground to pass a 60 mesh sieve, and dried again prior to analysis. Dried, sediment-laden GFFs were subsampled using a single-hole-punch tool. Samples were analyzed for total organic carbon (OC) and N by flash combustion elemental analyzer (Costech ECS 4010, Valencia, CA, USA). Molar ratios are used where N:C data are presented. Particulate organic carbon (POC) concentration (mg L−1 ) was calculated as: POC (mg L−1 ) = g of carbon (from NC analysis)/(area of filter subsampled (cm2 )/area of whole filter (cm2 )). The OC concentration of sediment was calculated as: %OC = POC (mg L−1 ) / TSS (mg L−1 ) \* 100.

Filtrate was analyzed for chloride (Cl−), nitrite (NO2 −), nitrate (NO3 −), and ammonium (NH4 +) by ion chromatography (Dionex DX−500, Sunnyvale, CA, USA). We measured ultraviolet absorbance at 254 nm (UVA254) using a BioMate3 spectrophotometer with a 1 cm cell in order to have a low-cost measurement to correlate with DOC and be able to convert dissolved organic carbon to specific UV absorbance (SUVA254) which provides insight into the nature of the organic matter present (e.g., high SUVA254 indicates dominance of humic matter and, thus, a soil origin for DOC). Total nitrogen (as TKN) was converted to ammonium cation by sulfuric acid digestion with mercuric oxide catalysis in a microwave digester (MARS X Press, CEM, Matthews, NC, USA). Samples were digested for one hour with a 15 min ramp time and 45 min hold time at 200 ◦C and analyzed for TKN by salicylate method on a Technicon Autoanalyzer III wet chemistry analyzer (Bran + Luebbe, Norderstedt, Germany). Dissolved inorganic nitrogen (DIN) was determined as sum of nitrite-N (NO2-N), nitrate-N (NO3-N), and ammonium-N (NH4-N). Dissolved organic nitrogen (DON) was calculated as total N minus DIN. We used DON:DOC and PON:POC (rather than the more commonly used DOC:DON and POC:PON) herein because it is statistically more robust due to the higher number (%OC) in the denominator.

Filtered, frozen water, suspended sediment, and soils were sent to the UC Davis Stable Isotope Facility and analyzed within 2 months of sample collection. Dried soils needed no preservation; water samples were preserved by filtering with a 0.2 micron filter and frozen according to UC Davis standard operating procedures. Water was analyzed for DOC and DOC δ 13C PDB using an OI Analytical model 1030 TOC Analyzer (OI analytical, College Station, TX) interfaced to a PDZ Europa 20–20 Isotope Ratio Mass Spectrometer (Sercon Ltd., Cheshire, UK) with a GD−100 gas trap interface (Graden Instruments). Stable isotopic composition of particulate organic carbon (POC; δ 13C) and particulate organic nitrogen (PON; δ 15N) were determined by flash-combustion elemental analysis coupled with isotope-ratio mass spectrometry (EA-IRMS; Elementar Analysensysteme GmbH, Hanau, Germany) interfaced to a PDZ Europa 20–20 isotope ratio mass spectrometer (Sercon Ltd., Cheshire, UK). Isotopic data for C and N were reported in standard delta notation (‰ per mil) relative to Pee Dee Belemnite (PDB) and air standards, respectively. Standard reference materials covering a range of 13C and 15N signatures were used to evaluate precision and accuracy of mass spectrometer response; analytical precision for δ 13C of DOC was 0.4‰; precision for δ 13C and δ 15N was tighter than 0.3‰.

#### *2.4. Data Analysis*

Average concentrations of measured stream water constituents were calculated for each gauged watershed for the length of the study. To determine the effect of storm-specific attributes on composition of dissolved and particulate material, individual flow events were delineated by the change in discharge. The beginning of an event within ES was defined as (1) any discharge >0.002 m3 s −1 or (2) rate of change in discharge >0.0002 m3 s −1 15 min−1 or (3) start of flow. The beginning of an event within perennial streams was defined as any discharge >0.002 m3 s −1 . Where there was a one-hour decline in discharge, followed by an increase in discharge, the event clock was reset, and a new event began. Linear Pearson correlations of log-transformed data were used to evaluate relationships between independent event variables and dependent surface water concentrations and characteristics. We chose event variables that were likely to affect the flux of sediment and OM within our watersheds. Independent event variables tested included: rate change in Q over time; water year date (number of days after 1 October (which was the beginning of the water year); time since flow began in hours; time since event initiation; and time before or after peak flow. Event variables are defined as: Precipitation = the amount of precipitation that fell in the previous 15 min; Q = instantaneous discharge at the time the sample was collected; PeakTime = time before or after peak flow; FlowTime = time since flow began in hours; RateChange = change in Q over time where positive values are on the rising limb and negative values are on the receding limb; EventTime = time since event initiation; WY Day = number of days after 1 October (beginning of water year). Stepwise multiple linear regression was used to explore relationships among event variables and surface water characteristics; data were log-transformed to meet assumptions of normality. The lm function in base R [48] was used to model constituent fluxes from independent event variables; step Akaike Information Criteria (stepAIC function) from the MASS package [49] was used to optimize the size of the model. Herein, we report the standardized coefficients with adjusted *R* 2 and whole-model *p*-values for purposes of examining the relationships among variables as well as the coefficients and equation used to develop predictive calculations. All estimates of mass fluxes were calculated using instantaneous flow-weighted averages across all storms and streams of within-event measurements.

#### **3. Results and Discussion**

#### *3.1. Physical and Chemical Attributes*

Over the 17 month study period, 16 events generated discharge in one or more streams (Table 2). Precipitation was monitored for 12 months (May 2010 through May 2011; Figure 2) and ranged from 21–110 mm per event with a mean of 55 mm. The maximum precipitation rate occurred at the onset of an event and ranged from 10–72 mm h−1 with a mean of 36 mm h−1 . Precipitation between January and May accounted for 53% of the annual total and 94% of the discharge events over 15 months. Event duration ranged from 1.0 h to 22.3 h with a mean of 8.0 h. Precipitation events triggered 233 flow

event samples ranging from a total of 14 to 82 samples per watershed per event (Table 1). Mean event discharge across all recorded events for ephemeral and perennial streams was 749 m3 event−1 and 7324 m3 event−1 , respectively.

**Table 2.** (A) Average instantaneous flow-weighted within-event values of selected dissolved and particulate stream water constituents across all sampling times. (B) Average values from potential sources.

Values are the mean ± standard error. Stream water values are the means of the base flow and event flow. The O-horizon and mineral soil values are from lysimeter solutions. %OC is the % of OC in sediment by mass. The *p*-values are from Mann–Whitney U tests between perennial and ephemeral streams.

Particulate constituent concentrations among ES were highly variable (Table 2-A). The TSS concentration ranged from 54.1 mg L−1 to 1713.0 mg L−1 , and POC concentrations ranged from 2.27 mg L−1 to 21.33 mg L−1 . Average ES TSS (566.1 mg L−1 ; instantaneous flow-weighted average across all streams and storms of within-event measurements; from Table 2) was higher than that of the perennial stream (101.1 mg L−1 ) by a factor of 5.6; thus, it appears that the perennial stream was a sink for suspended solids. Average ES POC concentration (9.30 mg L−1 ) was higher than that of the perennial stream (3.46 mg L−1 ) by a factor of 2.7.

Dissolved constituent values DOC, Cl, UVA254, and specific UV absorbance at 254 nm (SUVA254) were consistent among ES (Table 2-A). The UVA254 of ES was consistently higher than that of the perennial stream, probably due to the average DOC being (14.51 mg L−1 ) 2.5 times higher than that of the perennial stream (5.83 mg L−1 ). The SUVA254, on the other hand, was similar among perennial and ES. The ES Cl concentration was consistently lower than that of the perennial stream indicating that the perennial stream acquired Cl from another source(s).

We quantified concentrations of five potential source contributors to surface waters: precipitation, canopy throughfall, O-horizon soil leachate, mineral soil leachate, and groundwater (Table 2). The DOC concentration increased from precipitation through the canopy as throughfall and O-horizon leachate then declined through the soil profile reaching a minimum in the groundwater samples suggesting that it was removed or diluted with depth. The Cl concentrations generally increased from the canopy through the soil and into the groundwater. The Cl concentration of ephemeral event water was similar to that of canopy throughfall and O-horizon soil, whereas the Cl concentration in perennial event water was most similar to that of mineral soil leachate and groundwater. There was a strong linear correlation (Pearson's *R* = 0.923, *p* < 0.001) between UVA254 absorbance and DOC concentration, explaining the similar pattern between UVA254 and DOC. The UVA254 in precipitation was low but increased as

water passed through the forest canopy. Highest UVA254 levels were observed in O-horizon leachates; UVA254 decreased with depth as water passed through the soil profile. Mean UVA254 in ES water (0.56) was nearly three times that of perennial stream water (0.19) suggesting that C was lost or diluted during transport. The ratio of DON:DOC was highest in precipitation and groundwater and reached minimums in throughfall and mineral soil leachate while being relatively high in the O-horizon. Organic horizon water had the most depleted DOC-δ 13C, and there was a steady enrichment with soil depth (O, A, to groundwater). The DON:DOC had a relatively small standard error among potential sources and, hence, is probably a good tracer of DOM in stream water. On the other hand, the stable isotopic composition of DOC varied <2‰ among stream water and <4‰ among potential sources; therefore, DOC is indistinguishable from potential sources using stable isotopic signature alone (Table 2). Ephemeral stream DON:DOC and %OC was consistently lower than that of the perennial stream and most closely resembled that of O-horizon leachate, except within the disturbed clear-cut site in which DON:DOC resembled that of mineral soil. Perennial stream DON:DOC was similar to that of groundwater and precipitation. The POC:DOC was similar for ES and perennial stream waters; however, %OC of the perennial stream was significantly higher (*p* = 0.002) than that of ES. Mean ES PON-δ 15N was 1.4 (ranged 0.9 to 2.0) as compared to 1.0 in perennial stream water (Table 2). The PON:POC was similar among streams.

We calculated the area-weighted average instantaneous flux for all storms over the water year. Instantaneous flux is the flow-weighted flux at a given time within an event; thus, instantaneous flux could represent a point on the rising limb, peak flow, or the receding limb of an event. Values represent the means for all events in which we were able to collect water samples and quantify Q with rating curves (Table 3). Average Q was higher in the perennial stream than in the ES (Table 1); however, area-weighted discharge (Unit Q) in the ES was more than twice that of the downstream perennial stream (Table 3), even though the overall area was the same for the ES and perennial streams. Event sampling in ES was triggered by a 20 mm increase in flow over 15 min; thus, the Unit Q measurements only represent periods of high flow (storms) in the ES. There was higher Unit Q from the ES than at the downstream perennial stream gauge suggesting that some of the water discharged to the perennial stream was lost to groundwater before exiting the perennial stream gauge.

**Table 3.** Area-weighted average instantaneous flux of selected surface water constituents.

Means of unit Q and constituents were computed from within-event samples collected over an entire water year. Values in parentheses are one standard deviation of the mean. Some data were not normally distributed, so we present un-transformed data (A) and log-transformed data (B).

The ES exports were 2.0–6.5 times that of the perennial stream, depending upon the parameter (Table 3). Unit Q was higher for the ephemeral streams as a result of this value being calculated from the sample set from which the constituents were determined. The ephemeral catchments are very small—from 1.8–3.8 ha; yet during storms, their exports were significant which has implications for changes in storm events due to the climate forcing. Surface water, groundwater, and soil storage water are linked, so we used geochemical signatures from the stream exports and potential source areas in order to understand how the pools and fluxes might interact with changes in climate.

#### *3.2. Source of OM in Surface Waters during Storms: Active and Contributing Areas*

In order to get a better understanding of whether ephemeral and perennial stream waters were derived from similar source areas and were therefore connected during storm events, we compared the Cl and UVA254 of ephemeral and perennial stream event water to that of groundwater and soil leachate. Cl was used as a conservative tracer of potential contributors to stream water and UVA254 was used as a first-order evaluation of DOC inputs. We used UVA254 rather than DOC because of the strong linear correlation between the two in our data and because the larger UVA254 dataset gave us better control on variability. The UVA254 and Cl had a significant negative correlation (*R* = –0.595, *p* < 0.001) such that during an event, UVA254 was elevated and Cl concentration was depleted with respect to base flow conditions. The relationship between Cl and UVA254 amongst pools suggests that ES Cl and DOC are controlled by water routing through the soil profile, whereas the perennial stream is influenced by groundwater inputs.

We did not attempt to quantify the relative contributions of groundwater versus Horton- and Dunne-type overland flow; however, we used OM characteristics of ES and potential source areas to clarify from which pools ES water may have been derived (Figures 3 and 4). We did not attempt to trace N through the system; however, N and C cycling are coupled therefore we incorporated both components in source-attribution, and the ensuing discussion is framed in terms of OM. Ephemeral storm water DOC-δ 13C was intermediate between that of groundwater and soil solution dissolved OM with a DON:DOC composition close to that of soil solution and the *minimum* DON:DOC value (max N:C value) for groundwater (Figure 3).

**Figure 3.** Relationship of N:C and δδ 13C of DOM from groundwater, soil lysimeter, and ephemeral stream event water. Lysimeter samples were collected from below the litter layer, 10 cm into mineral soil, and 20 cm into mineral soil (undifferentiated in this plot).

A comparison between DOM and suspended sediment particulate organic matter (POM) of ES collected during storm events (Figure 4) illustrates the difference in source attributes for DOM and POM. Ephemeral stream DOM and suspended sediment POM clustered tightly in separate-δ 13C versus N:C mixing spaces: storm water plotted in a space that was close to intermediate soil dissolved OM values, whereas suspended sediment particulate DOM plotted closest to values of A- and Bhorizon mineral soils (Figure 4a). Soil leachate values plotted close to that of O-horizon POM (and to a lesser extend A-horizon POM) suggesting that O-horizon soils were an accessible source for labile DOM during storm events. Workers have reported that as much as 66% of DOC in stormflow from

headwater streams originated as saturated flow through the uppermost soil horizons within the riparian zone [50]. Estimated contributions from the upper portions of the soil profile (for lowland riparian soils) represent 64–86% of DOC flux in some watersheds (e.g., [51]). High DOC values were attributed to a larger proportion of wetland soils (i.e., larger potential DOC pool) and larger watershed area [51]. Storm-derived ES POM was more closely related to A-horizon and to a lesser extent B-horizon soils (Figure 4a–c), which suggests that these soils are mobilized or reworked enough to release POM during storm events.

δ δ **Figure 4.** N:C versus (**a**) δ 13C of stream DOM, soil solution DOM, POM from soils (O, A, B, C horizons), suspended sediment (SSPOM); (**b**) δ15N; and (**c**) %OC by mass. Stream DOM and suspended sediment POM values are for ES event samples only. In Figure 5a, purple circles with lower N:C represent DOM from A-horizon lysimeters; purple circles with higher N:C represent DOM from O-horizon lysimeters.

#### *3.3. Storm-Driven C Dynamics*

Samples collected during storm events provided evidence that the ES interacted with the soil profile during rainfall events. The next question we wanted to answer was whether the source or proportion of sources changed within a storm event. We selected storm events for which there was sufficient data to construct plots illustrating the behavior of TSS and OM over the course of an event (Figure 5). The TSS concentrations in ES were highest at the start of an event before peak Q and declined as the event progressed (Figure 5a). Peak perennial stream TSS concentration lagged behind ES peak TSS and coincided with an initial increase in Q. The N:C of particulate material remained constant as ES TSS decreased throughout an event (Table 2 and Figure 5b), thus there is a separation of C-poor material and C-rich material as the event progresses. There is a weak trend of enrichment in δ 13C with event time (Figure 5c) which suggests a within-event shift in C to progressively deeper sources.

δ δ **Figure 5.** Representative behavior of suspended sediment and POC during selected storm events. Blue bars are precipitation amounts. Solid lines indicate within-storm Q. Circles indicate measured values of discrete samples. PER = perennial stream. EPH = ephemeral stream. (**a**) Suspended Sediment (**b**) N:C of POM in suspended sediment (**c**) Change in δ 13C of POM in suspended sediment. Sampling start and stop times were controlled by preset criteria (see methods) and thus do not span the time range from start of precipitation to return to base flow. The single anomalously low δ 13C value may be due to entrainment of a piece of fresh leaf litter.

We attempted to track the responses of UVA254 (as a proxy for DOC) and Cl to storms in the perennial and ephemeral watersheds by constructing hysteresis plots from each of the 5 potential sources. In general, we found that the ES (31 December 2010, 30 April 2010, and 4 April 2011 storm panels in Figure 6) had higher UVA254 and lower Cl than the perennial stream (23 April 2010 storm panel). However, there were different trends among watersheds and storms.

**Figure 6.** Selected (typical) responses of UVA254 versus Cl concentration for several storms. The date of each event is listed above each pair of plots. Major ticks on the date/time axis of storm discharge plots are spaced every 3 h. The top plot of each pair includes a box inset with hydrographs (line representing Q with scale on left) and precipitation where available (bars with scale on right). Hydrographs are color-coded as to the catchment: brown = ephemeral 2; red = ephemeral 4; black = perennial. Precipitation data for April 23 and April 30 do not exist. Below the hydrograph insets, samples for a given event are plotted on a UVA versus Cl plot. Plots of UVA254 versus Cl are in alphabetical order by time stamp of sample collection. Letters represent the time point when a sample was collected (in alphabetical order by time stamp of sample collection); corresponding letters on the hydrograph and the UVA–Cl plot are the same sample. Within the plot sample schema: solid blue triangles are precipitation (precip), unfilled blue circles are throughfall (TF), unfilled purples circles are O-horizon leachate (O), solid purple triangles are A horizon leachate (A), and solid red triangles are groundwater samples (GW) (see legend in the center of the figure). Sources are average ± 95% confidence intervals of samples collected throughout the entire study.

An early season storm (31 December 2010) showed what appears to be flushing of a high Cl and UVA254 source that could be recently fallen litter. Cl concentration becomes diluted as the source of stream water shifted towards throughfall and then a mixture of O-horizon leachate and precipitation. The source then shifted towards throughfall and A-horizon leachate and eventually towards the groundwater source.

This pattern essentially repeats itself in a later season storm (30 April 2010) in the same watershed (ES2) without the flushing of recently fallen litter. For a late spring storm (4 April 2011) in ES4 (clearcut), there appears to be less influence from throughfall on the earliest storm samples, but the source of water eventually shifts through the soil and towards groundwater. The source of water from the perennial stream (23 April 2010) appears to be from deeper within the soil profile at the onset of a typical storm event relative to the ES. Eventually the source of water appears to be entirely derived from a groundwater source late in the event.

Differences in hysteresis at different times of the year suggest that there may be critical thresholds in some storm parameters that control both the source from which the sediment and OM are derived and the flux of sediment and OM within the system. Due to our sampling protocol, it is possible that the timing of our sampling in relation to the storm hydrograph varied with season; however, our sampling procedure was designed to eliminate as much bias as possible. We compared surface water characteristics to seven event variables related to storms: discharge, peak time, flow time, rate of change, event time, water year day, and precipitation (Table 4) to glean which event variables exerted the greatest control over surface water concentrations and fluxes in ES and perennial streams. Where distributions were not normal, we log-transformed the data. Data presented are those that were significant predictors of constituent flux. The larger the value, the greater the control. The sign determines whether the event variable was negatively or positively associated with the dependent variable.

Standardized coefficients (SCs) and standard error (SE) of predictors that were selected through AIC and stepwise procedures. Dependent variables: Precipitation = amount of precipitation that fell in the previous 15 min; Q = instantaneous discharge at the time the sample was collected; PeakTime = time before or after peak flow; FlowTime = time since flow began in hours; RateChange = change in Q over time where positive values are on the rising limb and negative values are on the receding limb; EventTime = time since event initiation; WY Day = # of days after 1 October (beginning of water year). Q, TSS, POC, PN, and DOC were all log-transformed to meet assumptions of normality.

As expected, Q and precipitation exerted the greatest control on TSS and OM for both streams. Precipitation is interesting in that all of the measured exports of ES increased with increasing precipitation, whereas the reverse was true in the perennial stream either as a result of dilution of exports as they reached the perennial stream or reduced carrying capacity in the perennial stream, or both. An increase in DOC with increasing precipitation and Q has been reported in other headwater studies (e.g., [26]). Rate Change, which is tied to the shape of the hydrograph (i.e., flashy behavior equates to a more rapid rate of change), was a significant predictor for TSS, POC, and PN in the perennial stream, but was weakly correlated to TSS in the ES. The other independent variables, PeakTime, FlowTime, EventTime, and WY DAY, were significant only for ES. Dependent variables decreased with increasing FlowTime and EventTime, which suggests that there was a pool of sediment, DOM, and POM which was accessed early in an event and became depleted over time. The positive

correlation with WY DAY for all dependent variables indicates that the time at which the pools were accessed may be significant in terms of interpreting potential responses to climate change. We recognize that WY DAY likely has a sinusoidal relationship and, thus, may not be a good predictor; however, we explored the use of sin transformation and it did not improve *R* 2 .

#### *3.4. The Water Year*

There is a growing body of evidence that the seasonal component to OM pools and fluxes in headwater streams (e.g., [24]) is linked to the wet–dry cycles that are prevalent in IRES and affects the quality and quantity of OM exports from IRES (e.g., [26,30,52,53]). If that is the case, then there should be a relationship between flux and timing within the water year, as the watershed transitions from hot, dry summer conditions to cool, wet winter conditions. We examined correlations between dissolved and particulate constituents and timing within the water year (WY Day; 1 October to 30 September) to see if there were effects associated with wet and/or dry periods (Table 4). There were significant effects between WY Day and flux of TSS, POC, PN, and DOC for ES, whereas these fluxes were not significant for the perennial stream. Flux of TSS, POC, and PN all had positive correlations with WY Day which indicates that flux of these particulates increases over the course of the water year. The most likely mechanism for this would be active downcutting and headward erosion which releases sediment and associated particulate OM into ES with successive storm events.

The negative correlation between water year date versus DOC indicates that ES DOC concentration was highest at the beginning of the "wet" season and gradually decreased over time which is consistent with data from a meta-analysis of 30 forested watersheds in the eastern United States documenting higher DOC exports at the end of a dry period [26]. Flushing of high Cl and UVA254 sources (e.g., throughfall and O-horizon material; Table 2) during an early season storm (Figure 6; 31 December 2010 storm panel) is not as evident later in the season. The late season storm (Figure 6; 30 April 2010 storm panel) has a Cl versus UVA254 signal which is similar to GW (i.e., wet season base flow) at the event onset, followed by throughfall and A-horizon soil. These relationships suggest that early in the season there was a larger pool of labile DOC available within and leached from the upper soil profile (e.g., O-horizon soil solution) and that over multiple flushing events the labile C pool was depleted. If storm events are "hot moments" for the flushing of stored DOC [28,40], then peak DOC efflux should occur at the onset of the wet season (late fall) and decrease throughout the wetter winter months during successive leaching events. The DOC pool will build during drier months such that events that occur after prolonged dry conditions will have high DOC concentrations [26,54]. Our data support belowground storage of a leachable C pool during the dry season and flushing of the soil C pool during storm events throughout the wet season, such that availability of DOC decreases to minimum at the beginning of the next dry season. In contrast, the flushing signal that is evident within the ES is not detectable within the perennial stream, most likely due to the dilution from groundwater inputs (Figure 6; 23 April 2010) or CO2 evasion.

#### *3.5. The Climate Connection*

Climate projections for the southeastern US indicate a slight increase in warm days and especially warm nights, which will lead to a longer duration of frost-free days, and an increase in drought length [55]. Concurrent with the overall increase in temperature, there is "high confidence" that there will be an increase in the number of annual days for which there are heavy (>75 mm) precipitation events [55] even though MAP is expected to be relatively stable [34]. A comparison of several models from the Mississippi River Basin also predict an increase in temperature which leads to a shift in the seasonality of the precipitation such that precipitation will decrease during summer months and increase during winter months [56]. Regional climate models from the central U.S. postulate a localized summertime warming minimum or "warming hole" [57] such that convective storms will result in precipitation that is greater than evapotranspiration, setting the stage for increases in soil moisture at

the end of the summer season. These models pose a conundrum for how climatic drivers will affect fluxes from IRES to downstream rivers.

If MAP does not change but event precipitation increases, then there must be less frequent events of higher intensity with increased drying time between events. This scenario would result in higher fluxes of Sediment and OM (until it is depleted seasonally) from the ephemeral watersheds with increasing precipitation and Q (Table 4) especially during the first flood pulses as has been reported for both particulate and dissolved OM e.g., [5,26]. The four ephemeral watersheds we gauged had a combined area of 11.6 ha, representing 36% of the source area contributing to the perennial stream at our gauge point, yet as small as they were, they discharged an average of 328 g s−1 ha−1 of OM (Table 3; POC + PN + DOC + DON) during storm events. The high standard deviations of these measurements indicate that the fluxes could be much higher or much lower. We used the existing relationship between precipitation/Q and the dependent variables within a regression function to provide an indication of how the watershed might respond to increases in precipitation or Q. The equation used was:

$$\begin{array}{rcl} \text{Ln(X)} = \text{B0} + \text{B1} & \ast \text{Ln(unitQ)} + \text{B2} \ast (\text{Peak time}) + \text{B3} \ast (\text{Floouttime}) + \text{B4} \\ \ast (\text{Rate of Change}) + \text{B5} \ast (\text{TimeEvent}) + \text{B6} \ast (\text{WY Day}) + \text{B7} \\ \ast (\text{Precip}) & \end{array}$$

where *X* is the dependent variable of interest (e.g., TSS, POC, PN, and DOC), B0 is the intercept, and B1–B7 are the coefficients for each of the predictor variables (0 if the predictor variable was not a significant predictor of the dependent). The resultant projected increase or decrease (%) in TSS and OM components given a 5% increase in precipitation or Q is provided in Table 5.

**Table 5.** Projected increase in TSS and OM components with 5% increase in precipitation or Q and coefficients from which they were derived.

Projected increase or decrease in each of the dependent variables are in percentage in the two columns to the right. ES = ephemeral streams. PER = perennial stream. Dependent variables: Precipitation = amount of precipitation that fell in the previous 15 min; Q = instantaneous discharge at the time the sample was collected; PeakTime = time before or after peak flow; FlowTime = time since flow began in hours; RateChange = change in Q over time where positive values are on the rising limb and negative values are on the receding limb; EventTime = time since event initiation; WY Day = # of days after October 1 (beginning of water year).

A century-long study of precipitation trends in forest lands of the Lower Mississippi River Alluvial Valley (LMRAV) reported precipitation increases (estimated to be 7%) at the south coastal area of the LMRAV and constant or a slight decrease at the middle area of the LMRAV [58]. Our study site was between the south and middle parts of the LMRAV and, therefore, a 5% increase in precipitation is a reasonable estimate. A 5% increase in stream discharge is considered as the extreme case because not all of the 5% increase precipitation contributes to discharge.

Where there is a 5% increase in precipitation, within ES we would expect increases across all particulate constituents. However, given the perennial stream sensitivity to precipitation there would most likely be a negative effect of on the flux of most constituents. When and how the precipitation occurs can change the outcome insofar as precipitation and Q are coupled. Higher magnitude storms, which may cause an increase in Q, will be the primary drivers for additional materials flux. Given a 5% increase in Q, there was an increase from 1.6% to 2.8% in average instantaneous flux of materials

from both the ES and the perennial stream (Table 5). For example, an increase of 0.63 × 103 m3 s −1 ha−1 in Q within ES would result in an increase of 69.0 g s−1 ha−1 average instantaneous flux of TSS and 7.8 g s−1 ha−1 of OM (Table 6). While increases appear to be small, it is important to remember these are instantaneous fluxes, so over the course of an event the flux will be much higher depending upon the change in the discharge-rating curve relative to historical data. Changes between ES and the perennial stream are similar, which suggests that materials flux from the ES continues in the perennial stream and that the perennial stream acts as a conduit for materials generated in the ES.

**Table 6.** Modeled increase in average within-event instantaneous flux of materials given a 5% increase in Q.

Instantaneous flux is Q × C, where Q is the instantaneous discharge and C is the concentration of the sample at the time of sample collection.

The hydrologic transition from dry to wet in these IRES could result in pulses of DOC release [30,31,59] followed by mineralization and CO2 release from the perennial stream [60], especially since rates increase with increasing catchment area [61]. Much of the DOC and POC in the ES derived from the soil profile (Figures 2 and 3). If, as a result of climate change, there is an increase in convective storms and greater soil moisture at the end of the summer season, then wetting and drying cycles of soils and streams will be moderated and there may be lower DOC production and flux simultaneously with higher POM due to the increased energy in the system to mobilize particulates and sediment.

If MAP does not change but rainfall shifts to winter months, then presumably the summer drought season would be lengthened. Increased drought during the warmer months followed by short-term high intensity events would result in increased downcutting and headward erosion followed by greater flux of TSS, POC, and PN (Figure 4 and Table 3). A compilation of climate models for an agricultural watershed in Mississippi predicted that an increase in summer temperatures would result in a decrease in runoff, but approximately 12% greater TSS loads, presumably due to the higher intensity precipitation events [62]. Perhaps of greater import is that prolonged drought may increase the areal extent of IRES [1] which means that there would be a greater time and area over which soil OM can accumulate. If Horton infiltration excess is the primary path for DOC to streams, then perhaps there would be no change beyond what is mobilized with sediment due to the erosion. If Dunne saturation excess is the primary source of DOC as is suggested by our data (Figures 2 and 4a), then there will be a greater time over which the pool can build and potentially larger releases to the stream network. Longer-term studies are needed to determine the impact of climate change on the flux of (OM) from IRES in the southeastern US.

#### **4. Conclusions**

Our first question was: what are the sources of sediment and OM in ES? DOM was derived from Dunne-type saturated flow through the upper part of the soil profile, particularly during the first flush of storms and after dry periods. Sediment and POM were derived from down cutting and headward erosion during the rising limb of storm events. Changes in processes that affect erosion are going to lengthen stream channel networks and increase the areal extent of ES.

Our second question was: what are the controls on flux of sediment and OM from very small (<4 ha) ES to downstream reaches? Precipitation, WY Day, and variables associated with the shape of the hydrograph all had linear relationships with sediment, DOM, and POM. No matter what happens to the other variables (e.g., precipitation, length of the water year, shape of the hydrograph), discharge will have an overriding effect on the flux of materials because there is a power log relationship between Q and materials flux. Any changes to Q (i.e., caused by higher-intensity storms) will have a large control over the flux.

The flux of OM in response to climate change is going to be dependent upon how the change affects the source of OM. Projected climate effects are: (1) less frequent, higher-intensity storms and/or (2) shift to heavier precipitation in winter and longer drought in summer. If there is a shorter water year but the same amount of precipitation, there needs to be a longer/larger stream network to accommodate the flow. Precipitation during storms is likely to result in more water transfer through the upper part of the soil profile in conjunction with increased channel incision. Higher intensity storms resulting in higher Q within ES channels will result in greater flux of materials to downstream perennial reaches.

Do small IRES matter in the context of climate change? Given the storm-driven flux of sediment and OM from these small ephemeral watersheds and the potential for that OM to be processed and evolved as CO2 under projected climate scenarios for the region, the answer is yes—IRES matter. There are three reasons that ES should be incorporated into climate-driven hydrological models: (1) ES are intimately connected to sources of OM; (2) ES are indirectly connected to precipitation through Q; (3) there is likely to be an expansion of the areal extent of ES. For these reasons we need to consider ES when modeling the response of hydrological systems to climate change.

**Author Contributions:** Conceptualization, J.D. and J.H.; methodology, J.D., J.H.; B.C., and C.M.; software, J.D. and J.H.; validation, J.D., J.H., B.C., and C.M.; formal analysis, J.H.; investigation, J.D., J.H., B.C., and C.M.; resources, J.D., J.H., and Y.O.; data curation, J.H.; writing—original draft preparation, J.D.; writing—review and editing, J.D., J.H., and Y.O.; visualization, J.D. and J.H.; supervision, J.D. and J.H.; project administration, J.D. and J.H.; funding acquisition, J.D. and J.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the US Geological Survey through Mississippi Water Resources Research Institute, grant number 06HQGR0094 06090820; the US Department of Agriculture Forest Service, grant number 06JV11330127186 07010040, the National Council for Air and Stream Improvement, Inc., grant number 07030270, and Weyerhaeuser Company, grant number 07100923.

**Acknowledgments:** Mississippi State University Forest and Wildlife Research Center provided laboratory support; Weyerhaeuser Company provided site support. Rick Maiers and the students at the Mississippi State University assisted with the project setup.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

# **Using AnnAGNPS to Simulate Runo**ff**, Nutrient, and Sediment Loads in an Agricultural Catchment with an On-Farm Water Storage System**

**Juan D. Pérez-Gutiérrez 1 , Joel O. Paz 2, \* , Mary Love M. Tagert 2 , Lindsey M. W. Yasarer 3 and Ronald L. Bingner 3**

Received: 5 October 2020; Accepted: 10 November 2020; Published: 12 November 2020

**Abstract:** On-farm water storage (OFWS) systems are best management practices that consist of a tailwater recovery (TWR) ditch used with a storage pond to provide irrigation water and improve downstream water quality. These systems have been increasingly implemented in the southeastern US, but the individual and cumulative effects of these systems on a watershed scale are unknown. In this study, the runoff, nutrient, and sediment loads entering a TWR ditch in an agricultural catchment were quantified, and contributing sources were identified using the annualized agricultural non-point source (AnnAGNPS) model. Fields with larger areas and soils with a high runoff potential produced more runoff. The volume of runoff exceeded the TWR ditch storage volume approximately 110 times, mostly during the winter and spring seasons. During years when corn and winter wheat were planted, NO3–N loads increased because these crops need nitrogen fertilization to grow. Planting winter wheat in priority subwatersheds reduced the total phosphorous (TP) and sediment loads by about 19% and 13%, respectively, at the TWR ditch inlet. Planting winter wheat can reduce runoff, TP, and sediment loads but also result in higher NO3–N loads. AnnAGNPS simulations quantified the benefits of an OFWS system to advance the understanding of their impact on water availability and quality at a watershed scale.

**Keywords:** tailwater recovery ditch; AnnAGNPS; BMP

#### **1. Introduction**

The global population has been projected to increase between 9.6 billion and 12.3 billion by 2100 [1]. As a result, agriculture will have to double its productivity to feed people [2,3]. Meeting this increasing food demand has led to the intensification of agriculture, which translates into the greater use of chemical inputs and pressure on soil and water resources. Agricultural activities often transfer excess fertilizers through non-point source (NPS) pollution into aquatic ecosystems causing devastating ecological and economical effects [4–8]. Therefore, a balance between increasing yields and mitigating adverse environmental impacts is crucial for the future of agriculture.

The implementation of best management practices (BMPs) has been recognized to significantly reduce NPS pollution from croplands to downstream waterbodies [9–11]. However, much effort is still devoted to quantifying the effectiveness and benefits of such practices, especially newer BMPs. Measuring BMP performance can be a challenge because their benefits vary spatially and temporally due to the heterogeneity of the landscape and seasonality of hydrological factors [12]. The challenge

is even greater when attempting to evaluate the effects of several BMPs combined [13–15], address the shift in BMP performance over time [16], and understand lag times in water quality response [17]. The use of watershed models can be a feasible alternative to quantify the effectiveness of BMPs while accounting for these challenges [16,18–24]. However, due to the diversity of BMPs that can be implemented over agricultural fields, information describing the performance of each BMP and the combined use of different BMPs is limited.

Structural BMPs such as tail-water recovery (TWR) ditches and agricultural ponds (i.e., on-farm water storage systems—OFWS) can improve downstream water quality by significantly reducing nutrient loads from agricultural watersheds [25–27]. The dual benefit of reducing nutrient pollution and supplying water for irrigation is important in areas such as the Lower Mississippi River Valley, where agricultural production strongly depends on irrigation. Farmers and landowners in this region are tasked with (1) reducing off-site movement of nutrients, which contributes to the hypoxic zone in the northern Gulf of Mexico, and (2) conserving water resources to slow the declining groundwater levels in the Mississippi River Valley Alluvial Aquifer (MRVAA), the primary source of water for irrigation. Consequently, OFWS systems have been implemented across the MRVAA, especially in areas with more severe groundwater declines. Although research has highlighted in-ditch and in-pond nutrient reductions, less attention has been placed on estimating the water and nutrient loads entering and exiting OFWS systems [28,29]. Evaluating the OFWS drainage area will quantify the impacts of these systems on a watershed scale, which will advance our understanding of the cumulative effects of these systems on water quality and groundwater decline as more OFWS systems are implemented throughout a watershed.

The annualized agricultural non-point source (AnnAGNPS) model [30,31] is a watershed model that has been designed to measure the impact of agricultural management practices on hydrological and water quality responses in watersheds. Most studies using the AnnAGNPS model demonstrate its performance after model calibration with field-observed data [22,32–39]. However, AnnAGNPS has also been used for similar purposes without conducting calibration with measured data. A study used the model with no calibrated parameters to estimate the runoff and sediment in an agricultural watershed within the Mississippi Delta Region [21]. They concluded that the model had an adequate ability to simulate monthly and annual runoff and sediment yield with no calibration process. This finding is of significant interest because it is difficult and costly to secure sufficient observed data to conduct calibration processes in watershed-scale modeling studies. To date, efforts to assess the benefits of tail-water recovery ditches using watershed scale models has been limited. This is, at least in part, because these ditches have multiple inlets. Therefore, the collection of measured runoff and water quality data (and especially continuous measured runoff and water quality data) is usually impractical and expensive, making it a challenge to collect a sufficient quantity of measured data on even one OFWS system. However, estimating the effects of OFWS systems at various scales is critical to stakeholders and agencies to better identify where these systems can be implemented to improve water quality and relieve pumping pressure on groundwater.

In this study, the AnnAGNPS model was used without calibration with field data to simulate runoff, sediment, and nutrient loads and identify the main contributing areas draining to a TWR ditch established as part of an OFWS system within Porter Bayou Watershed, Mississippi. The quantification of the water, nutrient, and sediment loading is an essential step in understanding the water quality and quantity benefits of OFWS systems when implemented in agricultural watersheds, as well as how the management of these systems might be altered to improve performance. Predictions on the volume of water captured by OFWS systems can also help determine their potential for alleviating local flooding. Finally, this study advances the incorporation of an OFWS system and the associated watershed impacts into the AnnAGNPS model, with the goal of improving the ability of AnnAGNPS to model these systems on a larger scale.

#### **2. Materials and Methods**

#### *2.1. Study Site*

This investigation was conducted in the Porter Bayou Watershed (PBW; 33◦26′39"–33◦51′38" N, 90◦48′54"–90◦31′34" W) in the Mississippi Delta Region (MDR), an intensively farmed area located in northwest Mississippi, USA (Figure 1). For the AnnAGNPS simulation, a small 214 ha watershed within the PBW located north of Indianola, MS, in Sunflower county, was selected, and it included Metcalf farm and the surrounding area that drains into the OFWS system outlet at M3 (Figure 1). The major crops grown from 2012 to 2016 were soybeans, corn, and rice (Table 1). Winter wheat was planted during 2013 and 2014 in four and two fields, respectively. The watershed is dominated by the following soil types: Forestdale, Tensas, Dundee, Pearson, and Dowling, which are primarily poorly drained soils and prone to produce high runoff. In addition, the watershed is relatively flat with surface elevations ranging from 130 to 135 m. The average temperature ranged from −8.3 ◦C in winter 2014 to 32.2 ◦C in summer 2012 (Figure 2a). The lowest total seasonal precipitation was 99 mm observed in summer 2015, while the highest was 579 mm in spring 2014 (Figure 2b). Average annual precipitation was 1308 mm, and in 2013 and 2014 this average was exceeded by 242 mm and 114 mm, respectively. Meanwhile, 2012 and 2015 were drier years with a total annual rainfall of 1100 and 1133 mm, respectively.

**Figure 1.** Map of the study area showing the agricultural watershed with the on-farm water storage system investigated in the Mississippi Delta Region. The blue arrows represent the runoff flow direction towards the outlet. (M1: tailwater recovery (TWR) inlet; M2: TWR mid-canal; M3: TWR outlet; labels within fields indicate the subwatershed identification used by the annualized agricultural non-point source (AnnAGNPS) model).

**Table 1.**Basic characteristics and land use including the planting dates of the subwatersheds.

TRNAR: turn area; WW: winter wheat. **Figure 2.** Average temperature and total precipitation by season during the period 2012–2016 at the study site. (**a**) Boxplots of the average temperature. (**b**) Bar chart showing the total precipitation. Winter: January–March; Spring: April–June; Summer: July–September; Fall: October–December. Boxplots were set at 90th (the upper whisker), 75th (the upper quartile), 50th (the median), 25th (the lower quartile), and 10th (the lower whisker) percentiles. Outliers were considered those observations 1.5 times beyond the 25th and 75th percentile and are shown as grey circles.

The OFWS system consists of a trapezoidal-shape TWR ditch and an elongated agricultural pond which have a combined storage volume of 128,020 m3 (TWR ditch: 13,320 m3 ; pond: 114,700 m3 ). Water flows from north to south through the ditch, which is 818.8 m long and 1.8 m deep on average; the pond is 2.4 m deep with a surface area of 4.45 ha. Runoff is routed to the single outlet pipe (33◦39′35.6" N, 90◦39′11.9" W) set at 1.2 m above the canal bed (Figure 1). The system was designed according to National Resources Conservation Service (NRCS) [40] guidelines, and more information about its characteristics can be found at Pérez-Gutiérrez, Paz and Tagert [25].

#### *2.2. Model Description*

AnnAGNPS is a physical process model developed to simulate the runoff, sediment, nutrient, and pesticide yields at a daily time step in small watersheds. The model divides the watershed into subwatersheds based on homogeneous physical characteristics such as soil type, land use, and land management. AnnAGNPS is a continuous simulation model and has been primarily developed to evaluate the impacts of different agricultural management practices on watersheds. As with other physical process watershed-scale models, the major input data are climate, land characterization, field operations, chemical characteristics, and feedlot operations. A detailed description of the model can be found in Geter and Theurer [30], Bingner, et al. [41], Bosch, et al. [42], Theurer and Cronshey [43], Cronshey and Theurer [44].

#### *2.3. Model Input*

A detailed field survey was conducted to identify the field boundaries and collect elevation data required by the model. Eighteen fields were identified as subwatersheds (or cells), and associated reaches were defined for routing runoff to the outlet (M3) in the AnnAGNPS model (Figure 1). Because all fields were precision leveled, they were defined as homogeneous drainage areas or subwatersheds. The delineation of the watershed was done manually with the aid of GoogleTM Earth and geographic information system (GIS) technologies. Parameters describing the subwatersheds such as area, average elevation, and average land slope were determined from the field survey (Table 1). Parameters representing the time of concentration and travel time were computed from data provided by the field reconnaissance following USDA-SCS [45] methods, modified by Theurer and Cronshey [43]. Soil data and physical properties were obtained from the Soil Survey Geographic (SSURGO) database [46]. Although 15 soil types were identified throughout the watershed, the dominant soil type was determined for each subwatershed as required by the model, using GIS operations (Table 1). Crop planting dates were obtained from the Sunflower county USDA–NRCS office. A crop management schedule was assigned to each field according to the typical management practices conducted in the MDR (Table 2). Irrigation was included in the crop management schedule, starting in late May or June and ending in August. The SCS curve number is an important model parameter used to estimate runoff. Table 3 shows the curve numbers used in the model, based on different land use categories and hydrologic soil types in the watershed. Weather data were recorded automatically at 15-min intervals from March 2012 to December 2016 by a WatchDog 2700 Weather Station (Spectrum® Technologies, Inc., Aurora, IL) located 9.2 km southeast of the outlet (M3). The weather data were subsequently processed to create daily time scale files as required by the AnnAGNPS model. Data from the Parameter-elevation Regressions on Independent Slopes Model (PRISM) Climate Group (http://www.prism.oregonstate.edu/explorer/) were used to fill any gaps in rainfall records. Other missing weather data such as daily maximum and minimum temperature, dew point, wind velocity, and solar radiation were patched using data from the Moorhead, MS climate station, which is located 22 km south of the study site and managed by the National Oceanic and Atmospheric Administration.

**Table 2.** Typical crop management operation for the crops planted in agricultural watersheds within the Mississippi Delta Region (MDR) used in this study.

**Table 3.** Soil Conservation Service (SCS) curve numbers selected for runoff estimation relative to cropland at the agricultural watershed. Source: USDA-SCS [47].

AnnAGNPS was used to simulate the runoff, sediment, and nutrient loads entering the TWR ditch. The model was set to run for two initialization years (2010–2011) to establish antecedent conditions before beginning the simulation period over the next five years (2012–2016). As described by Bosch, Theurer, Bingner, Felton and Chaubey [42], AnnAGNPS outputs are predefined by the user for the watershed source of interest (subwatersheds, reaches, among others). The model produces event-based output as well as monthly and annual summaries of hydrologic and water quality parameters. The model partitions nitrogen and phosphorus load into sediment-bound and dissolved fractions. This study used dissolved nitrate nitrogen (NO3–N) and total phosphorus (TP) to represent NO3–N and TP, respectively, and employed AnnAGNPS to estimate the runoff, sediment, NO3–N, and TP loads entering a TWR ditch.

#### **3. Results and Discussion**

Table 4 presents the average annual runoff, nutrient and sediment loads entering the TWR ditch per unit area, as well as the contributing fields and total area draining to each point—M1, M2, or main outlet M3—in the TWR ditch. The load values represent the aggregate contributions divided by the area of the fields draining into the TWR ditch at each outlet. These values provide a unique picture of the impact of each source area regardless of size. Annual runoff volume per unit area entering the ditch at mid-point M2 was slightly higher than the runoff entering at the M1 or M3 by about 600 m3 ha−1 yr−1 , even though M2 had the smallest drainage area. This difference is attributed to the runoff potential of the fields draining to M2. Subwatershed C4 is the field covering the third largest area of all subwatersheds, and soil in this field has the highest runoff potential (hydrologic soil group D). In addition, reaches transporting runoff flow directly into the ditch from subwatersheds C4 and C5. This direct flow path minimizes the amount of water allowed to infiltrate or evaporate and results in a greater volume of water draining into the ditch.

**Table 4.** Contributing fields to each channel reach, production of average annual runoff and load of nutrients and sediments entering the TWR ditch.

The NO3–N load per unit area was higher at M1 as most of the fields that were planted with corn and winter wheat during the simulation period are located upstream of M1. The TP load entering the ditch at M3 tended to be higher. Subwatershed C3, one of the three fields that drained to M3, was the only field that was simulated as a turn area. Typically, turn areas have compacted soils and do not have vegetation. This means that the soil in this field is highly susceptible to erosion and the off-site movement of nutrients by rainfall–runoff, which might explain the higher TP load simulated by AnnAGNPS at M3. Sediment loading did not show much variation through the TWR ditch segments, indicating that the sediment load might be equally distributed over the total area contributing to the TWR ditch outlets.

#### *3.1. Spatial Variation*