PROTECTING WATER QUALITY AND INCREASING RESILIENCY OF CROP PRODUCTION WITH CLIMATE-SMART DRAINAGE STRATEGIES By Babak Dialameh A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Biosystems Engineering – Doctor of Philosophy 2024 ABSTRACT Accurate phosphorus (P) load estimation in subsurface drainage water is critical for assessing the field-scale efficacy of conservation practices and minimizing the environmental impact of P loss from subsurface-drained fields to freshwater bodies like the Great Lakes. Also, to build a more resilient crop production system in a changing climate, it is crucial to understand how future weather patterns affect subsurface drainage design and whether subirrigation will be needed in the future for crop production. In this study, we used high-frequency P concentration measurements to investigate P transport dynamics and evaluate the effects of water sampling strategies on the uncertainty of P load estimation. We also evaluated the change in the optimum drain spacing from using historical (1994-2023) and future (2030- 2059) weather data and assessed the efficacy of subirrigation to alleviate yield reduction due to drought stress in southeast Michigan, USA. We used the HydroCycle-PO4 instrument to measure total reactive P (TRP) concentration at high resolution from a subsurface-drained field, evaluating four hypothetical water sampling strategies: time-proportional discrete and composite sampling, and flow-proportional discrete and composite sampling. All strategies underestimated TRP load compared to the reference dataset, with underestimation ranging from 0.2% to 51% depending on the interval and method. Flow- proportional strategies provided more accurate estimates of cumulative P load with fewer samples, capturing more samples at higher flow rates. TRP concentration showed a transport-limited chemodynamic pattern, increasing with flow during events. A 1% increase in drainage discharge led to a 1.36% increase in TRP load, with flow events contributing 89% of P loss and the highest 7.7% of flow transporting 75% of the TRP load. Most flow events (30 out of 36) displayed a flushing effect, indicating preferential flow as a pathway for TRP loss. High-frequency P sampling suggested that management and conservation practices should target flow events to reduce P loss. A total of 27 general circulation models with a moderate greenhouse gas emission scenario (shared socioeconomic pathway 2-4.5) were used for climate projections. Simulations were performed using the DRAINMOD model, and the optimum drain spacing was determined based on the maximum average annual return on investment. Results showed that the projected 30-year average annual precipitation is not expected to change significantly while that of the temperature will increase by 2.5°C in the future. Future optimum drain spacings for depths of 75 cm and 125 cm were found to be 300 cm and 600 cm wider than historical spacings, respectively. On average, there was a 23% decrease in the 30-year average annual drainage discharge, attributed to an average 17% increase in evapotranspiration. Drought stress is projected to be the primary cause of yield loss in the future, due to increased temperatures and an average 8% deeper water-table depth. Subirrigation shows high potential in reducing year-to-year crop yield variability in the future (decreasing the coefficient of variation for the yield from 0.26 to 0.06, on average) and increasing yield by up to 31%. In the past, subirrigation initiation was feasible in late June with a weir depth of 70 cm. However, in the future, subirrigation may have to start sooner in early to mid-June, and weir depth may have to be shallower ranging from 65 cm to 55 cm to be more advantageous. In conclusion, a wider drain spacing, providing reduced drainage intensity, along with subirrigation may be needed in the future to mitigate crop yield loss from drought stress. In conclusion, using future climate projections may allow wider drain spacing to retain more water and reduce drought stress. Without such data, implementing subirrigation could help mitigate potential yield losses from drought. Copyright by BABAK DIALAMEH 2024 To my parents v ACKNOWLEDGEMENTS I would like to express my deepest gratitude to my advisor, Dr. Ehsan Ghane, for his invaluable guidance, unwavering support, and insightful feedback throughout the journey of completing this dissertation. His mentorship has been instrumental in shaping my research and academic growth. I am also immensely thankful to my esteemed committee members, Dr. Mohamed Youssef, Dr. Timothy Harrigan, and Dr. Wei Zhang, for their expertise, constructive criticism, and encouragement, which have greatly enriched the quality of this work. My heartfelt appreciation extends to my dedicated colleagues, Jason Piwarski, Samantha Smith, Youssef AbdalAal, Sami Shokrana, Manal Askar, Josue Kpodo, and Joseph Letavis, for their collaboration, camaraderie, and support. Their diverse perspectives and shared experiences have enriched my academic journey immeasurably. I extend my sincere thanks to the staff at BAE for providing a professional and friendly environment that facilitated my research and academic endeavors. I am also grateful to Michigan State University, Michigan Department of Agriculture and Rural Development, and Environment, Great Lakes, and Energy – State of Michigan for their financial support. Last but certainly not least, I would like to express my deepest gratitude to my parents for their unwavering love, encouragement, and sacrifices. Their endless support and belief in my abilities have been the driving force behind my academic pursuits. To all those who have contributed to this endeavor in ways big and small, I am profoundly grateful. vi TABLE OF CONTENTS LIST OF ABBREVIATIONS ........................................................................................................ ix CHAPTER 1 EFFECT OF WATER SAMPLING STRATEGIES ON THE UNCERTAINTY OF PHOSPHORUS LOAD ESTIMATION IN SUBSURFACE DRAINAGE DISCHARGE ............1 1.1. Abstract ................................................................................................................................ 2 1.2. Introduction and Literature Review ..................................................................................... 3 1.3. Materials and Methods ......................................................................................................... 5 1.3.1. Study site ....................................................................................................................... 5 1.3.2. Precipitation data .......................................................................................................... 6 1.3.3. Measurement of drainage discharge and phosphorus concentration ............................ 6 1.3.4. Calculating the reference hourly TRP load ................................................................... 8 1.3.5. Description of sampling strategies for TRP load estimation ........................................ 8 1.3.6. Calculating uncertainty of TRP load estimation ......................................................... 11 1.4. Results and Discussions ..................................................................................................... 12 1.4.1. Relationship between hourly drainage discharge and hourly TRP concentration ...... 12 1.4.2. Time-proportional discrete sampling scenarios .......................................................... 13 1.4.3. Time-proportional composite sampling scenarios ...................................................... 15 1.4.4. Flow-proportional discrete sampling scenarios .......................................................... 17 1.4.5. Flow-proportional compositing sampling scenarios ................................................... 18 1.4.6. Practical application of the findings ........................................................................... 19 1.5. Conclusion ......................................................................................................................... 22 REFERENCES ......................................................................................................................... 24 CHAPTER 2 INVESTIGATION OF PHOSPHORUS TRANSPORT DYNAMICS USING HIGH-FREQUENCY MONITORING AT A SUBSURFACE-DRAINED FIELD IN THE WESTERN LAKE ERIE BASIN ..................................................................................................28 2.1. Abstract .............................................................................................................................. 29 2.2. Introduction ........................................................................................................................ 30 2.3. Materials and Methods ....................................................................................................... 32 2.3.1. Study site ..................................................................................................................... 32 2.3.2. Precipitation data ........................................................................................................ 33 2.3.3. High-frequency measurement of drainage discharge ................................................. 33 2.3.4. High-frequency measurement of phosphorus concentration ...................................... 34 2.3.5. Estimating the TRP load ............................................................................................. 35 2.3.6. Determining phosphorus transport dynamics ............................................................. 35 2.3.7. Description of hypothetical sampling strategies ......................................................... 36 2.3.8. Delineating events ....................................................................................................... 37 2.3.9. Quantifying hysteresis and flushing indices ............................................................... 37 2.4. Results ................................................................................................................................ 38 2.4.1. Relationship between TRP concentration and drainage discharge ............................. 38 2.4.2. Relationship between TRP load and drainage discharge ............................................ 40 2.4.3. Effect of sampling strategy on the accuracy of P transport dynamic analysis ............ 42 2.4.4. Summary of flow events ............................................................................................. 43 2.4.5. Contribution of event flow to P loss ........................................................................... 43 vii 2.4.6. Hysteresis and flushing patterns ................................................................................. 45 2.5. Discussion .......................................................................................................................... 48 2.5.1. Relationship between TRP concentration and drainage discharge ............................. 48 2.5.2. Contribution of baseflow and event flow to P transport ............................................. 49 2.5.3. C-Q Relationship ........................................................................................................ 50 2.5.4. L-Q Relationship ......................................................................................................... 51 2.5.5. Hysteresis and flushing indices ................................................................................... 52 2.5.6. Effect of sampling strategy on the accuracy of P transport dynamic analysis ............ 54 2.6. Conclusions ........................................................................................................................ 55 REFERENCES ......................................................................................................................... 57 CHAPTER 3 SUBSURFACE DRAINAGE DESIGN AND SUBIRRIGATION AS CLIMATE- SMART STRATEGIES FOR RESILIENT CROP PRODUCTION ............................................63 3.1. Abstract .............................................................................................................................. 64 3.2. Introduction ........................................................................................................................ 65 3.3. Materials and Methods ....................................................................................................... 69 3.3.1. Site description............................................................................................................ 69 3.3.2. Measurement of drainage discharge and water-table depth ........................................ 70 3.3.3. The DRAINMOD model setup to simulate hydrology, relative crop yield, and subirrigation .......................................................................................................................... 71 3.3.4. Evaluating the effect of climate change on drainage design ....................................... 78 3.3.5. Economic analysis of various drainage designs.......................................................... 78 3.3.6. Selecting the optimum drain spacing .......................................................................... 80 3.3.7. Suitability of subirrigation .......................................................................................... 81 3.3.8. Statistical analysis ....................................................................................................... 82 3.4. Results and Discussions ..................................................................................................... 83 3.4.1. Comparing temperature and precipitation the past and future .................................... 83 3.4.2. Results of calibration and validation of DRAINMOD ............................................... 87 3.4.3. Optimum drainage design for past and future............................................................. 89 3.4.4. Effect of climate change on hydrology and crop yield ............................................... 91 3.4.5. Efficiency of subirrigation practice in the past and future .......................................... 98 3.5. Conclusions ...................................................................................................................... 102 REFERENCES ....................................................................................................................... 105 APPENDIX A: CHAPTER 1.......................................................................................................111 APPENDIX B: CHAPTER 2 .......................................................................................................121 APPENDIX C: CHAPTER 3 .......................................................................................................127 viii LIST OF ABBREVIATIONS AACPI Average Annual Crop Production Income AADSC Average Annual Drainage System Cost AAR Average Annual Return on Investment ACPC Annual Corn Production Cost CMIP5 Coupled Model Intercomparison Project phase 5 CMIP6 Coupled Model Intercomparison Project phase 6 CP DEP DRP Corn Price Depreciation Per Year Dissolved Reactive Phosphorus DWR Drainage Water Recycling ET FI Evapotranspiration Flushing Index GCM General Circulation Model GS HI ICDS KGE MAE MME MR Growing Season Hysteresis Index Initial Cost of the Drainage System Kling-Gupta Efficiency Mean Average Error Multi-model Ensemble Maintenance Rate NCCS NASA Center for Climate Simulation ix NEX-GDDP NASA Earth Exchange Global Daily Downscaled Projections NGS Non-Growing Season NOAA National Oceanic and Atmospheric Administration NSE P Nash-Sutcliff Efficiency Phosphorus PBIAS Percent Bias PY RH RIAI RY Potential Yield Relative Humidity Interest Rate on Average Investment Relative Yield RZWQM Root Zone Water Quality Model SSP TRP Shared Socioeconomic Pathway Total Reactive Phosphorus WTD Water Table Depth x CHAPTER 1 EFFECT OF WATER SAMPLING STRATEGIES ON THE UNCERTAINTY OF PHOSPHORUS LOAD ESTIMATION IN SUBSURFACE DRAINAGE DISCHARGE 1 1.1. Abstract Accurate phosphorus (P) load estimation in subsurface drainage water is critical to assess the field-scale efficacy of conservation practices. The HydroCycle-PO4 instrument measures real- time total reactive P (TRP) concentration without the need for sample filtration, thereby enabling comparative evaluation of different sampling strategies. The main objective of this study was to evaluate the effects of water sampling strategies on the uncertainty of P load estimation. Hourly TRP concentration and hourly drainage discharge measurements formed the reference P load dataset. Four hypothetical water sampling strategies were evaluated: (a) time-proportional discrete sampling, (b) time-proportional composite sampling, (c) flow-proportional discrete sampling, and (d) flow-proportional composite sampling. All sampling strategies underestimated TRP load compared with the reference dataset, regardless of whether the underestimation was statistically significant. Total reactive P load underestimation changed from 0.2 to 51% as time-proportional discrete sampling intervals increased from 3 h to 14 d. Total reactive P load underestimation changed from 12 to 43% as the time-proportional compositing scenario increased from 1 to 7 d, each with one aliquot per day. In the case of flow-proportional discrete sampling scenario, the lowest (0.6%) and the highest (–5.1%) uncertainties were observed when 1- and 5-mm flow intervals were used. The relative error based on the results provided by the flow-proportional composite sampling ranged from 0.2% when using 1-mm flow interval to –6.7% when using 5- mm flow interval. In conclusion, the flow-proportional sampling strategies provided a more accurate estimate of cumulative P load with fewer number of samples because a greater portion of samples were taken at higher flow rates compared with time-proportional sampling strategies. 2 1.2. Introduction and Literature Review Subsurface drainage has enhanced economical crop production in poorly drained soils (Evans & Fausey, 1999). However, drain pipes can rapidly transfer nutrients to surface water bodies (Fausey et al., 1995; Ghane & Askar, 2021). Phosphorus (P) is a key nutrient transported by subsurface drainage that causes eutrophication and overstimulates the growth of harmful algal blooms (Chen et al., 2018; Ding et al., 2018; Wilson et al., 2019). There is an urgent need to reduce P transport from subsurface-drained fields to downstream water bodies. Edge-of-field monitoring of conservation practices is the first step to develop strategies that reduce P load from subsurface drainage discharge. However, edge-of-field monitoring requires efficient water sampling strategies to accurately calculate P load (Daniels et al., 2018; Dantas Mendes, 2020; Harmel et al., 2018). The water sampling strategy affects uncertainty in P load estimates (Birgand et al., 2011; Moatar et al., 2013; Shih et al., 1994). There are different water sampling strategies for drainage discharge (e.g., water sampling on a daily, weekly, fortnightly, and monthly basis), which are selected depending on budget and equipment limitations (Table A1.1). Long sampling intervals might be sufficient for some indicators, contaminants, or purposes, but they can fail to account for short-term P variability (Bowes et al., 2015; Fones et al., 2020). Long sampling intervals fail to capture P concentration fluctuations, particularly during precipitation events when drainage discharge rapidly rises and recedes (Luo et al., 2012). Selecting the most appropriate sampling strategy is essential for accurate P load calculation in drainage discharge with temporal P variations (Villa et al., 2019). Therefore, there is a need to evaluate the uncertainty associated with different water sampling strategies in estimating P load in drainage discharge. For subsurface drainage application, we only found Williams et al. (2015) that had investigated the effect of sampling strategies on the accuracy of dissolved reactive phosphorus 3 (DRP) load estimation in drainage discharge. The authors used sub-daily (2-4 h) intervals during storm events at some sites, and other sites were samples daily and the daily concentrations were interpolated to estimate hourly DRP concentrations at other sites. The authors recommended a minimum time interval to accurately estimate the annual DRP load when using a time-proportional sampling strategy. However, to our knowledge, there has been no study that has evaluated the accuracy of P load estimation in drainage discharge using flow-proportional sampling strategies. Flow-proportional sampling strategies are commonly used at edge-of-field monitoring. Therefore, there is a need to evaluate the accuracy of P load estimation of flow-proportional sampling strategies. In situ sensors have been developed for measuring high-frequency nutrient concentration while not needing filtration of the water sample. Liu et al. (2020) used such sensors for measuring high-resolution nitrate concentration, whereas we did not find any study that had used high- frequency P sensors to evaluate subsurface drainage application. This study aimed to evaluate the effects of different sampling strategies on P load estimation uncertainty: time-proportional discrete sampling, time-proportional composite sampling, flow- proportional discrete sampling and flow-proportional composite sampling. The value of this study is that it informs decision-making about the suitable sampling strategies needed to evaluate the P- reducing benefits of conservation practices. 4 1.3. Materials and Methods 1.3.1. Study site This study was conducted on a private farm located in Lenawee County, Michigan, USA, from January 2019 to July 2020 (Figure 1.1). The drainage area was 22.5 ha. The dominant soil type was Blount loam (fine, illitic, mesic Aeric Epiaqualfs), which was classified as a poorly drained soil. Subsurface drainage pipes were installed at an average 0.75-m depth with 12-m drain spacing. This system was under conventional free drainage during the study period. Figure 1.1. Geographic location and drainage layout of the study site. The cropping system during the study period was continuous corn with cereal rye as a winter cover crop. The conservation tillage practice was no-till. Manure was surface broadcasted at a rate of 43.2 kg-P/ha in January 2019 and 5.3 kg-P/ha in December 2019. Commercial fertilizer (product formulation 9-18-9 with sulfur) containing 2.9 kg-P/ha was applied in May 2020. 5 1.3.2. Precipitation data Precipitation data were collected using the microclimate sensor suite ATMOS-41 (METER Group, Inc., USA). This device is built-in with a high-resolution precipitation sensor (0.017-mm resolution). This sensor only measures rainfall; therefore, we used the snow water equivalent data from the National Oceanic and Atmospheric Administration (NOAA) weather station located 13.5 km away at the Adrian Lenawee County Airport to complete the precipitation data. 1.3.3. Measurement of drainage discharge and phosphorus concentration 1.3.3.1. Hourly drainage discharge measurement We combined two methods to measure hourly drainage discharge from January 2019 to July 2020. The first method measured drainage discharge with a metal-edge sharp-crest 45° V-notch weir (Agri Drain Corp., USA), which was installed in a 25-cm Agri-Drain water-level control structure. This method was used only when two conditions were satisfied: water flowed inside of the V-notch (i.e., did not exceed the weir capacity); and water-level in the structure's downstream chamber did not exceed the V-notch apex height. A HYDROS-21 water-depth sensor (METER Group) hourly measured the head of water inside the V-notch weir. Then, a calibrated V-notch equation was used to calculate the hourly drainage discharge based on the head. The second method measured hourly drainage discharge using a TIENET-350 area-velocity sensor (Teledyne ISCO, USA), which was installed inside a pipe located downstream of the control structure. This method was used only when either of these two conditions was met: water flowed over the V-notch weir (exceeding the weir capacity); and water level in the downstream chamber of the structure exceeded the height of the V-notch apex. In both methods, the drainage discharge rate was estimated using hourly area-velocity and water-depth measurements. The area-velocity sensor provided high flow rates and the V-notch weir provided low flow rates. 6 1.3.3.2. High-frequency phosphorus concentration measurement using the HydroCycle-PO4 We used an in situ HydroCycle-PO4 (Sea-Bird Scientific, USA) to measure high-resolution P concentration at a 2-hour interval (Figure A1.1). The HydroCycle-PO4 measures P concentration based on a heteropoly molybdenum-blue complex with phosphate that can be detected colorimetrically (Murphy & Riley, 1962). As the P concentration measurement is conducted on an unfiltered sample, the values represent total reactive phosphorus (Rice et al., 2017). The term reactive refers to the inorganic form of P that is readily bioavailable for uptake. The HydroCycle-PO4 has an average accuracy of -5 to 20% at the lower and upper ends of the operating range, respectively (Snazelle, 2018). Johengen et al. (2017) performed an assessment of precision of the HydroCycle-PO4 by computing the standard deviations and coefficients of variation of the five replicate measurements for five different concentration trials (from 0.01 to 0.4 mg/L). The standard deviation of the mean ranged from 0.0005 to 0.0020 mg/L across the five trials, and the coefficient of variation ranged from 0.14 to 5.78 percent. Also, the sensor has a maximum and minimum P detection limit of 1.2 and 0.002 mg/L, respectively. The TRP concentration in drainage discharge was measured from January 2019 to July 2020. TRP concentration measurements were retrieved from the data logger every week and processed using CycleHost software (Sea-Bird Scientific). Finally, the 2-hour interval TRP concentrations were linearly interpolated to obtain hourly TRP concentrations. Although the sensor used a single copper screen mesh with 7.5-µm onboard filters to reduce sediment intake, we observed sediment accumulation inside the unit. Thus, we performed weekly cleaning of the instrument to maintain high-quality data. New cartridges were installed every 3 to 4 months to maintain proper functioning of the HydroCycle-PO4. 7 1.3.4. Calculating the reference hourly TRP load Hourly TRP concentrations (section 1.3.3.2) and hourly drainage discharge measurements (section 1.3.3.1) were used to calculate the reference hourly TRP load as Loadref = ∑ QiCi n i=1 (1) where Loadref is the reference TRP load (kg/ha), Qi is the hourly drainage discharge (m3/h), and Ci is the hourly TRP concentration (mg/L). To adjust for units, a conversion factor of 4.44 × 10−5 was used. Any day that had missing data for any number of hours was eliminated from the analysis. 1.3.5. Description of sampling strategies for TRP load estimation 1.3.5.1. Time-proportional discrete sampling The reference hourly TRP concentration was subsampled to create eight hypothetical time- proportional discrete frequencies: 3-h, 6-h, 12-h, 1-day, 2-day, 3-day, 7-day, and 14-day intervals (Figure 1.2). The load for each sampling frequency was estimated by multiplying the TRP concentration of the selected sample and the cumulative drainage discharge during the relevant sampling interval as follows: Loadest = 4.44 × 10−5 × (∑ QjCj n j=1 ) (2) where Loadest is the estimated load (kg/ha), Qj is the cumulative drainage discharge during the sampling interval (m3/h), and Cj is the TRP concentration in the middle of the sampling interval (mg/L). The starting point for our artificial subsampling was January 11, 2019, the day that we started collecting TRP concentration data. The TRP load for each sampling strategy was estimated with 24 iterations due to the possibility of different starting times during a day (from 00:00 to 23:00). The calculation uncertainty for each load is minimized by performing several iterations, as outlined in Williams et al. (2015). 8 There were some gaps in the TRP concentration dataset. Malfunctioning and removal of the instrument for regular maintenance accounted for 211 days of no concentration data. No-flow condition also occurred for 10 days due to freezing of water inside the control structure during winter and drying of water inside the control structure. As a result, 305 days of flow data out of 566 days were used in the analysis. 1.3.5.2. Time-proportional composite sampling Time-proportional composite sampling involves the collection of numerous aliquots collected at regular intervals over a specified period. These aliquots form a composite sample that is representative of the sampling period. To investigate the effect of time-proportional compositing scenarios on the accuracy of TRP load estimation, we estimated the daily TRP load with various compositing scenarios based on the reference dataset described in section 13.4. The reference hourly TRP concentration dataset was subsampled to create 20 hypothetical time-proportional composite sampling scenarios: 1-, 2-, 3-, and 7-day composites, each with 1, 2, 4, 6, and 8 aliquots per day (Figure 1.2). The TRP concentration in the composite sample was considered equal to the average of the TRP concentrations over the sampling interval. Then, we calculated the TRP load by multiplying the average TRP concentration with its associated drainage discharge during the sampling interval. We did not account for any potential effect of bottles remaining in the automated sampler until retrieval and any effect that delayed filtering may have on P concentration. 9 1.3.5.3. Flow-proportional discrete sampling Flow-proportional sampling is a strategy widely used in water quality monitoring programs (Kladivko et al., 1991; Schleppi et al., 2006; Stone et al., 2000; Ulén & Persson, 1999; Wang et al., 2003). In this strategy, the water sample is collected when a specified volume or depth of water passes the monitoring point. Four volumetric depths of 1.0, 2.0, 3.0 and 5.0-mm were selected as flow intervals (Figure 1.2). These hypothetical manual flow-proportional discrete sampling scenarios were implemented on the reference hourly TRP concentration dataset. 1.3.5.4. Flow-proportional composite sampling The reference hourly TRP concentration was the basis to create hypothetical flow- proportional composites. Since 1-L bottles are commonly used in automated samplers, each bottle can hold six 150-ml aliquots without risk of bottle overflow for a total volume of 900 ml per bottle. To estimate the TRP concentration for each sample bottle, the reference hourly TRP concentration was subsampled for each aliquot based on four flow depth intervals of 1.0, 2.0, 3.0 and 5.0-mm. Then, the six aliquots per bottle were averaged to generate the TRP concentration for each sample bottle (Figure 1.2). 10 Figure 1.2. Sampling strategies and sampling scenarios investigated in this study. 1.3.6. Calculating uncertainty of TRP load estimation The relative error quantifies the uncertainty in load estimates (Harmel et al., 2006; Williams et al., 2015). The relative error was calculated as e(%) = ( Loadest−Loadref Loadref ) × 100 (3) where e is the uncertainty (%), Loadest is the estimated TRP load based on the specific sampling interval or compositing scenario (kg/ha), and Loadref is the reference hourly TRP load (kg/ha). The difference between the estimated load and reference load was calculated after each iteration of the TRP load estimation, which resulted in a distribution of uncertainty values. The bias and precision of the P load estimate were determined from this distribution. We used the distribution (e50) median as a measure of bias and computed precision as the difference between the 95th (e95) and 5th (e5) percentiles of the distribution (Williams et al., 2015). 11 1.4. Results and Discussions 1.4.1. Relationship between hourly drainage discharge and hourly TRP concentration The hourly TRP concentrations in our study varied from 0.007 to 1.161 mg/L with an average of 0.136 mg/L. These hourly TRP concentrations were generally higher than those reported by previous studies, which used larger sampling intervals ranging from daily to monthly (Daigh et al., 2015; Daly et al., 2017; Tiemeyer et al., 2009). Our higher TRP concentrations may be explained by the finer hourly dataset that did not require immediate sample filtration (Harmel et al., 2006, 2018; Massri et al., 2021), as opposed to the coarser temporal resolution used in previous studies. The previous studies collected samples with longer sampling intervals that could have missed higher P concentrations occurring at peak flows. Hourly drainage discharge rates in our study varied between 0.016 and 0.062 mm/h. Differences in soil, climate, drainage system, and agronomic practices also may explain the difference in P concentration between our study and previous studies (King et al., 2015). Temporal variations in hourly drainage discharge and hourly TRP concentration indicated that TRP concentration tended to increase during high flow events, whereas it remained steady during baseflow (Figure 1.3). A similar relationship was reported by other studies (Bende-Michl et al., 2013; Stamm et al., 1998; Vidon & Cuadra, 2010). Our combined results show the importance of high-frequency P sampling to accurately measure P concentration variation during high flow events. Subsequently, an accurate measure of P concentration variation leads to an accurate evaluation of the P transport dynamics (i.e., rapid P changed during storm events, event hysteresis pattern, flushing and flashiness). The hourly TRP concentration and drainage discharge measurements during two other periods are provided in Figure A1.2 and A1.3. 12 Figure 1.3. Hourly TRP concentration and drainage discharge in a subset of reference dataset from March 18, 2020, to May 17, 2020. 1.4.2. Time-proportional discrete sampling scenarios The accuracy of TRP load estimation was considerably affected by sampling frequency (Table 1.1) (Figure A1.4). The estimated TRP load decreased as sampling interval increased from 3 h to 14 days, which led to underestimation of TRP load compared to the reference hourly load (Table 1.1). This underestimation became considerable when the sample collection frequency was longer than one day. For example, the estimated TRP load based on 14-day sampling interval was 0.67 kg/ha, which was 51% less than the reference TRP load. Sampling intervals longer than one day often miss sharp increases in P concentration during the rising limb of the event flow hydrograph and the peak P concentration during storm events, thus, they do not accurately represent P variation (section 1.4.1). Therefore, the shorter the sampling interval, the better representation of the rapid variation of P concentration during storm events (Figure A1.5). Bias increased from –0.6 at 3-h sampling interval to –47.8 at 14-day sampling interval (Figure 1.4 and Table 1.1). The precision of P load estimation generally decreased as the sampling 13 interval increased, which was due to high temporal variations in P concentration in drainage discharge. As P concentration dramatically varies from baseflow to event flows, different sampling intervals might produce considerably different P load estimates, especially when using long sampling intervals. The total number of collected samples varied from only 22 samples, when using a 14-day interval, to 2440 samples, when using a 3-h hour interval (Table 1.1). The decrease in the total number of collected samples from 2440 to 22 resulted in an increase of error from –0.2% to –51.0%. Therefore, the total number of collected samples in the time-proportional discrete sampling strategy considerably affected the P load estimates. Table 1.1. Uncertainty indicators (relative error, bias, and precision) for TRP loads estimated using different time-proportional discrete sampling intervals for the entire period of the study. Common Sampling method Automated Automated Manual Sampling interval 1 h (reference) 3 h 6 h 12 h 24 h 48 h 72 h 7 days 14 days Total number of samples that needs to be analyzed 7320 TRP load (kg/ha) Relative error in TRP load (%) 2440 1220 610 305 154 100 44 22 1.37 1.36 1.34 1.29 1.20 1.10 0.99 0.78 0.67 – –0.2 –1.4 –5.3 –12.2 –19.2 –27.2 –42.7 –51.0 Bias (e50) – –0.6 –1.5 –4.8 –12.0 –20.4 –28.3 –41.1 –47.8 Precision (e5 to e95) – –1/–0.3 –2.8/–1.3 –7.8/–2.6 –18.8/–5.6 –34.5/–2 –40.4/–13.7 –59/–27.7 –72/–36 14 Figure 1.4. Bias and precision of TRP load estimation using different sampling intervals (e50, median representing bias; e5 and e95, 5th and 95th percentiles representing precision, respectively). 1.4.3. Time-proportional composite sampling scenarios The number of aliquots did not considerably affect precision and bias in P load estimation (Table 1.2), which was consistent with the study of Harmel & King (2005). Compositing scenarios using 1 and 8 aliquots had the lowest and highest precision, respectively. Precision did not considerably differ for composite samples from 1 to 7-day intervals. By contrast, bias was considerably affected by compositing interval. The 1-day composite with 1 aliquot is analogous to discrete sampling with 1-day sampling interval. The 1-day composite had the lowest average bias of –12.5, whereas the 7-day composite had the highest average bias of –42.8 across varying numbers of aliquots. These results indicate that the 1-day composite with 1 aliquot per day is a reliable and cost-effective strategy for P load estimation. The 1-day composite had the lowest uncertainty in TRP load estimation compared to 2, 3, and 7-day composites (Table 1.2). The highest relative error (43.0%) was observed for the 7-day composite, suggesting that longer compositing intervals considerably underestimated TRP load. The total number of composite samples varied from 44 to 305 when implementing 1-day and 7-day 15 composite scenario, respectively (Table 1.2). The decrease in the total number of collected samples from 305 to 44 resulted in an increase of error from –12.4% to –43%. Therefore, the total number of collected samples in the time-proportional composite sampling strategy considerably affected the P load estimates. It is important to note that results of this section assume there is no error from delayed filtering of sample. We used HydroCycle-PO4 instrument that provides real-time P concentration, while under field conditions, samples remain in the automated sampler until they are retrieved. Thus, the delay in sample filtration after sample collection generates uncertainty in DRP concentration measurements (Harmel et al., 2006, 2018; Massri et al., 2021). Therefore, the actual error from the time-proportional composite sampling strategy is expected to be slightly higher than the ones reported in this study because P concentration has been shown to decrease over time if not filtered. 16 Table 1.2. Uncertainty indicators (bias and precision) for the TRP loads estimated by different time- proportional composite sampling scenarios for the entire period of the study. Relative error was calculated for the TRP loads were averaged across different number of aliquots (1, 2, 4, 6, and 8 aliquots per day). Compositing scenario Bias (e50) Total Number of samples That Needs to Be Analyzed 7320 TRP Load (kg/ha) 1.36 Relative Error in TRP load (%) – Precision (e5 to e95) Reference 1-day composite 2-day composite 3-day composite 7-day composite Number of aliquots per day – 1 2 4 6 8 1 2 4 6 8 1 2 4 6 8 1 2 4 6 8 305 1.20 –12.4 154 1.10 –19.7 100 0.99 –27.7 44 0.78 –43.0 – –12.0 –12.4 –12.6 –12.7 –12.6 –19.2 –20.0 –19.9 –20.0 –19.7 –28.1 –27.8 –27.7 –27.8 –27.4 –42.7 –42.7 –43.1 –42.7 –42.8 – –18.8/–5.6 –14.3/-11.1 –12.9/-12.2 –12.9/-12 –12.7/-12.5 –23/-16.1 –20.3/-18.9 –20.3/-18.9 –20.3/-18.9 –20.1/-19.5 –30/-23.7 –28.0/-26.9 –28.2/-26.6 –28.0/-26.9 –27.9/-27.3 –45.6/-33.6 –43.0/-42.5 –43.6/-42.8 –43.0/-42.5 –42.8/-42.6 1.4.4. Flow-proportional discrete sampling scenarios The analysis of flow-proportional discrete sampling scenarios showed that the flow interval can affect the accuracy of TRP estimations. However, this effect may not be substantial. The shorter 1-mm flow interval scenario underestimated the TRP load for 0.2% while there was a 5.1% underestimation in TRP load estimation while using the longer 5-mm flow interval scenario (Table 1.3). All scenarios underestimated the TRP load for the study period. Overall, the accuracy of flow- proportional discrete sampling in TRP load estimation was acceptable with any sampling interval. High accuracy is obtained by using this sampling strategy because a greater portion of samples are taken at higher flow rates (Ulén & Persson, 1999) (Figure A1.6). Therefore, the shorter the flow interval, the better representation of the variation of P concentration during a storm event. 17 Even though there were no major differences among relative errors of the flow-proportion discrete sampling strategies, the number of samples that needed to be analyzed were considerably different, from 389 samples when using 1-mm flow interval to 85 samples when using 5-mm flow interval. The six-aliquot compositing also followed the same trend of requiring a smaller number of samples for the 1-mm interval compared to the 5-mm flow interval. Therefore, the larger sampling intervals can provide reasonably well estimates of P load while costing much less for water analysis. Table 1.3. Uncertainty indicators (relative error, bias, and precision) for TRP loads estimated using different flow-proportional sampling scenarios (discrete sampling and six-aliquot compositing). Sampling Scenario Reference Discrete sampling Six-aliquot compositing Flow interval 1 mm 2 mm 3 mm 5 mm 1 mm 2 mm 3 mm 5 mm Total Number of Samples That Needs to Be Analyzed 7320 389 190 127 85 66 34 24 16 TRP load (kg/ha) 1.37 1.37 1.33 1.32 1.30 1.36 1.27 1.30 1.34 Relative error in TRP load estimation (%) - –0.2 –3.1 –2.9 –5.4 –0.5 –7.2 –4.8 –1.9 1.4.5. Flow-proportional compositing sampling scenarios The accuracy TRP load estimation was not sensitive to the flow interval when using flow- proportional compositing sampling. Generally, no trend was observed between the relative error in TRP load estimation and flow interval (Table 1.3). The highest error (-7.2% underestimation) was produced by 2-mm flow interval with 34 composite sample. The least error also was observed when using 1-mm flow interval with 66 analyzed samples. 18 1.4.6. Practical application of the findings According to the results, total number of samples that needs to be collected and analyzed was substantially different when selecting various sampling intervals (Tables 1 to 3). Although a fewer number of samples were collected and analyzed in flow-proportional sampling strategies than time-proportional sampling strategies, lower uncertainty in TRP load estimation was observed when using flow-proportional sampling strategies. For example, if 5% error in P load estimation is the target, 85 and 610 samples need to be analyzed when using flow-proportional and time- proportional sampling strategies, respectively. Although time-proportional sampling strategies are simpler to be implemented, the flow-proportional sampling strategy best represents the cumulative P load in drainage discharge because a greater portion of samples are taken at higher flow rates (Figure 1.5) (Harmel et al., 2003). Therefore, the flow-proportional sampling strategies provided a more accurate estimate of cumulative P load at a lower analysis cost compared to time-proportional sampling strategies. Figure 1.5. An example of the difference in sample timing between time-proportional and flow- proportional sampling. 19 During the decision-making process, the suitable sampling strategy is the one that provides a balance between the purpose of the study and the budget. Flow-proportional sampling strategies can estimate cumulative P loss during a certain period with a fewer number of water samples compared to time-proportional sampling strategies. However, they fail to record the P concentration in drainage discharge at each time step especially during storm events when P concentration fluctuates rapidly. Generally, a shorter flow interval is suitable for smaller drainage areas, when P load is the objective, because a smaller drainage area conveys a less volume of water. Similarly, a longer flow interval is suitable for larger drainage areas. High-frequency time-proportional discrete sampling strategy is needed, if P transport dynamics (rapid P changed during storm events, event hysteresis pattern, flushing and flashiness) are of interest to the monitoring program. If other low-frequency sampling strategies are used to assess P transport dynamics, they cannot capture rapid changes in P concentration during peak flow, thereby creating bias in the results. If accurate P transport dynamics is required, automated sampler or real-time sensors like HydroCycle-PO4 are needed. Due to the limited capacity of automated samplers, they cannot hold many of bottles. Therefore, water samples need to be retrieved from the field frequently, especially during storm events, or when time or flow interval is short. Besides, real-time sensors eliminate any error related to delayed filtering. Under field conditions, depending on the selected sampling frequency, time-proportional discrete sampling can be performed either with manual grab samplers (i.e., sampling intervals longer than a day) or using automated samplers (i.e., sub-daily frequencies). An automated sampler can be used for all four sampling strategies. However, flow-proportional strategies rely on an external flow/depth sensor to send a signal to the data logger, and then the data logger sends another signal to the automated sampler to trigger sampling. This means that there is a greater chance of 20 losing concentration data due to failure or malfunction of one of the sensors because there are more parts involved in the sampling process. A time-proportional strategy has a lower chance of losing the concentration data because sampling does not rely on the flow/depth sensor. A simple diagram showing the cost and relative error of the four sampling strategies for the study period is illustrated in Figure 1.6. For all sampling strategies, except for the 1 to 14-day discrete time-proportional sampling, an automated sampler is needed, so we included the cost of one automated sampler (assuming $5000) for each strategy. For all strategies, a standard flow sensor is required to estimate the P load (assuming $5000) for each strategy. We assumed $10 per sample as the cost of chemical analysis. The flow-proportional sampling strategies produced almost the same accuracy for estimating P load as the high-resolution time-proportional sampling strategies (3h, 6h, and 12h) with lower cost (Figure 1.6). The cost analysis is presented in Table A1.2. 21 Figure 1.6. A simple diagram of the differences in cost and relative error of the four sampling strategies over 305 days. All sampling strategies included a cost for an automated sampler, except for the 1 to 14-day discrete time-proportional sampling. The cost of $10 per sample was included in the analysis. The cost of a flow sensor was included for all strategies. 1.5. Conclusion It is critical to obtain accurate P load estimates from subsurface-drained fields to comparatively evaluate conservation practices. In this study, we evaluated the accuracy of the following sampling strategies in P load estimation in drainage discharge: 1) time-proportional discrete sampling, 2) time-proportional composite sampling, 3) flow-proportional discrete sampling and 4) flow-proportional composite sampling. Our study resulted in eight key conclusions: • All sampling strategies underestimated TRP load compared to the reference dataset. This underestimation should be considered in P budget calculations. • As time-proportional discrete sampling intervals increased from 1 day to 14 days, 22 underestimation of TRP load changed from 12% to 51%. The uncertainty in TRP load estimation declined as the sampling interval decreased. • The underestimation of TRP load changed from 12% to 43% as the time-proportional compositing scenario increased from 1-day to 7-day composite, each with one aliquot per day. • The number of aliquots (1, 2, 4, 6, and 8 aliquots per day) collected for the 1- to 7-day time- proportional composite did not considerably affect the accuracy of load estimations. • In the case of flow-proportional discrete sampling strategies, both discrete and composite sampling produced accurate results (the relative error from –0.2% to –7.2%). • The flow-proportional sampling strategies provided a more accurate estimate of cumulative P load at a lower analysis cost compared to time-proportional sampling strategies. The purpose of the monitoring project should dictate the sampling strategy. If calculating the cumulative P load during a certain period is the main purpose of a monitoring program, flow- proportional sampling strategies (either discrete or composite) can be used to provide fairly accurate results with a smaller P budget underestimation. If P transport dynamics (rapid P changed during storm events, event hysteresis pattern, flushing and flashiness) are of interest, high-frequency time- proportional sampling strategies are recommended. The high-frequency time-proportional sampling strategies can be performed with automated samplers or real-time sensors. This study provides new insight about the accuracy of each sampling strategy as stand-alone method, thereby helping the user make the best decision for choosing a sampling strategy. Even though this study deals with TRP, the findings apply to DRP because both are comprised of dissolved form of P. 23 REFERENCES Bende-Michl, U., Verburg, K., & Cresswell, H. P. (2013). High-frequency nutrient monitoring to infer seasonal patterns in catchment source availability, mobilisation and delivery. Environmental Monitoring and Assessment, 185(11), 9191–9219. https://doi.org/10.1007/s10661-013-3246-8 Birgand, F., Faucheux, C., Gruau, G., Moatar, F., & Meybeck, M. (2011). Uncertainties in assessing annual nitrate loads and concentration indicators: Part 2. Deriving sampling frequency charts in Brittany, France. Transactions of the ASABE, 54(1), 93–104. https://doi.org/10.13031/2013.36263 Bowes, M. J., Jarvie, H. P., Halliday, S. J., Skeffington, R. A., Wade, A. J., Loewenthal, M., Gozzard, E., Newman, J. R., & Palmer-Felgate, E. J. (2015). Characterising phosphorus and nitrate inputs to a rural river using high-frequency concentration-flow relationships. Science of the Total Environment, 511, 608–620. https://doi.org/10.1016/j.scitotenv.2014.12.086 Chen, M., Ding, S., Chen, X., Sun, Q., Fan, X., Lin, J., Ren, M., Yang, L., & Zhang, C. (2018). Mechanisms driving phosphorus release during algal blooms based on hourly changes in iron and phosphorus concentrations in sediments. Water Research, 133, 153–164. https://doi.org/10.1016/j.watres.2018.01.040 Daigh, A. L. M., Zhou, X., Helmers, M. J., Pederson, C. H., Horton, R., Jarchow, M., & Liebman, M. (2015). Subsurface Drainage Nitrate and Total Reactive Phosphorus Losses in Bioenergy-Based Prairies and Corn Systems. Journal of Environmental Quality, 44(5), 1638–1646. https://doi.org/10.2134/jeq2015.02.0080 Daly, K., Tuohy, P., Peyton, D., Wall, D. P., & Fenton, O. (2017). Field soil and ditch sediment phosphorus dynamics from two artificially drained fields on poorly drained soils. Agricultural Water Management, 192, 115–125. https://doi.org/10.1016/j.agwat.2017.07.005 Daniels, M. B., Sharpley, A., Harmel, R. D., & Anderson, K. (2018). The utilization of edge-of- field monitoring of agricultural runoff in addressing nonpoint source pollution. Journal of Soil and Water Conservation, 73(1), 1–8. https://doi.org/10.2489/jswc.73.1.1 Dantas Mendes, L. R. (2020). Edge-of-field technologies for phosphorus retention from agricultural drainage discharge. Applied Sciences, 10(2), 634. https://doi.org/10.3390/app10020634 Ding, S., Chen, M., Gong, M., Fan, X., Qin, B., Xu, H., Gao, S. S., Jin, Z., Tsang, D. C. W., & Zhang, C. (2018). Internal phosphorus loading from sediments causes seasonal nitrogen limitation for harmful algal blooms. Science of the Total Environment, 625, 872–884. https://doi.org/10.1016/j.scitotenv.2017.12.348 24 Evans, R. O., & Fausey, N. R. (1999). Effects of Inadequate Drainage on Crop Growth and Yield. In Agricultural Drainage (Vol. 38, pp. 13–54). https://doi.org/10.2134/agronmonogr38.c2 Fausey, N. R., Brown, L. C., Belcher, H. W., & Kanwar, R. S. (1995). Drainage and water quality in Great Lakes and cornbelt states. Journal of Irrigation and Drainage Engineering, 121(4), 283–288. https://doi.org/10.1061/(asce)0733-9437(1995)121:4(283) Fones, G. R., Bakir, A., Gray, J., Mattingley, L., Measham, N., Knight, P., Bowes, M. J., Greenwood, R., & Mills, G. A. (2020). Using high-frequency phosphorus monitoring for water quality management: a case study of the upper River Itchen, UK. Environmental Monitoring and Assessment, 192(3). https://doi.org/10.1007/s10661-020-8138-0 Ghane, E., & Askar, M. H. (2021). Predicting the effect of drain depth on profitability and hydrology of subsurface drainage systems across the eastern USA. Agricultural Water Management, 258, 107072. https://doi.org/10.1016/j.agwat.2021.107072 Harmel, R. D., Cooper, R. J., Slade, R. M., Haney, R. L., & Arnold, J. G. (2006). Cumulative uncertainty in measured streamflow and water quality data for small watersheds. Transactions of the ASABE, 49(3), 689–701. https://doi.org/10.13031/2013.20488 Harmel, R. D., King, K., Busch, D., Smith, D., Birgand, F., & Haggard, B. (2018). Measuring edge-of-field water quality: Where we have been and the path forward. Journal of Soil and Water Conservation, 73(1), 86–96. https://doi.org/10.2489/jswc.73.1.86 Harmel, R. D., & King, K. W. (2005). Uncertainty in measured sediment and nutrient flux in runoff from small agricultural watersheds. Transactions of the ASAE, 48(5), 1713–1721. https://doi.org/10.13031/2013.20005 Harmel, R. D., King, K. W., & Slade, R. M. (2003). Automated storm water sampling on small watersheds. Applied Engineering in Agriculture, 19(6), 667–674. https://doi.org/10.13031/2013.15662 Johengen, T., Purcell, H., Tamburri, M., Loewensteiner, D., Smith, G. J., Schar, D., McManus, M., & Walker, G. (2017). Performance verification statement for Sea-bird Scientific HydroCycle phosphate analyzer. In Allicance for Coastal Technologies. https://doi.org/http://dx.doi.org/10.25607/OBP-289 King, K. W., Williams, M. R., Macrae, M. L., Fausey, N. R., Frankenberger, J., Smith, D. R., Kleinman, P. J. A., & Brown, L. C. (2015). Phosphorus transport in agricultural subsurface drainage: A review. Journal of Environmental Quality, 44(2), 467–485. https://doi.org/10.2134/jeq2014.04.0163 Kladivko, E. J., Van Scoyoc, G. E., Monke, E. J., Oates, K. M., & Pask, W. (1991). Pesticide and Nutrient Movement into Subsurface Tile Drains on a Silt Loam Soil in Indiana. Journal of Environmental Quality, 20(1), 264–270. https://doi.org/10.2134/jeq1991.00472425002000010043x 25 Liu, W., Youssef, M. A., Birgand, F. P., Chescheir, G. M., Tian, S., & Maxwell, B. M. (2020). Processes and mechanisms controlling nitrate dynamics in an artificially drained field: Insights from high-frequency water quality measurements. Agricultural Water Management, 232, 106032. https://doi.org/10.1016/j.agwat.2020.106032 Luo, Y., Arnold, J., Allen, P., & Chen, X. (2012). Baseflow simulation using SWAT model in an inland river basin in Tianshan Mountains, Northwest China. Hydrology and Earth System Sciences, 16(4), 1259–1267. https://doi.org/10.5194/hess-16-1259-2012 Massri, Z., Nunn, A. N., Ghane, E., Seybold, A., Piwarski, J. R., Adams, J., & Asher, J. (2021). Water sample holding times at appropriate temperature, filtering, and storage conditions for accurate analysis of soluble reactive phosphorus. ASABE Annual International Virtual Meeting 2101116. https://doi.org/https://doi.org/10.13031/aim.21 Moatar, F., Meybeck, M., Raymond, S., Birgand, F., & Curie, F. (2013). River flux uncertainties predicted by hydrological variability and riverine material behaviour. Hydrological Processes, 27(25), 3535–3546. https://doi.org/10.1002/hyp.9464 Murphy, J., & Riley, J. P. (1962). A modified single solution method for the determination of phosphate in natural waters. Analytica Chimica Acta, 27, 31–36. https://doi.org/10.1016/S0003-2670(00)88444-5 Rice, E. W., Baird, R. B., & Eaton, A. D. (2017). Methods for the examination of water and wastewater. American Public Health Association, American Water Works, 23. Schleppi, P., Waldner, P. A., & Stähli, M. (2006). Errors of flux integration methods for solutes in grab samples of runoff water, as compared to flow-proportional sampling. Journal of Hydrology, 319(1–4), 266–281. https://doi.org/10.1016/j.jhydrol.2005.06.034 Shih, G., Abtew, W., & Obeysekera, J. (1994). Accuracy of nutrient runoff load calculations using time-composite sampling. Transactions of the ASAE, 37(2), 419–429. https://doi.org/10.13031/2013.28093 Snazelle, T. (2018). Laboratory evaluation of the Sea-Bird Scientific HydroCycle-PO4 phosphate pensor. https://pubs.usgs.gov/of/2018/1120/ofr20181120.pdf Stamm, C., Flühler, H., Gächter, R., Leuenberger, J., & Wunderli, H. (1998). Preferential Transport of Phosphorus in Drained Grassland Soils. Journal of Environmental Quality, 27(3), 515–522. https://doi.org/10.2134/jeq1998.00472425002700030006x Stone, K. C., Hunt, P. G., Novak, J. M., Johnson, M. H., & Watts, D. W. (2000). Flow- proportional, time-composited, and grab sample estimation of nitrogen export from an eastern coastal plain watershed. Transactions of the ASAE, 43(2), 281–290. https://doi.org/10.13031/2013.2703 26 Tiemeyer, B., Kahle, P., & Lennartz, B. (2009). Phosphorus losses from an artificially drained rural lowland catchment in North-Eastern Germany. Agricultural Water Management, 96(4), 677–690. https://doi.org/10.1016/j.agwat.2008.10.004 Ulén, B., & Persson, K. (1999). Field-scale phosphorus losses from a drained clay soil in Sweden. Hydrological Processes, 13(17), 2801–2812. https://doi.org/10.1002/(SICI)1099- 1085(19991215)13:17<2801::AID-HYP900>3.0.CO;2-G Vidon, P., & Cuadra, P. E. (2010). Impact of precipitation characteristics on soil hydrology in tile- drained landscapes. Hydrological Processes, 24(13), 1821–1833. https://doi.org/10.1002/hyp.7627 Villa, A., Fölster, J., & Kyllmar, K. (2019). Determining suspended solids and total phosphorus from turbidity: comparison of high-frequency sampling with conventional monitoring methods. Environmental Monitoring and Assessment, 191(10), 605. https://doi.org/10.1007/s10661-019-7775-7 Williams, M. R., King, K. W., Macrae, M. L., Ford, W., van Esbroeck, C., Brunke, R. I., English, M. C., & Schiff, S. L. (2015). Uncertainty in nutrient loads from tile-drained landscapes: Effect of sampling frequency, calculation algorithm, and compositing strategy. Journal of Hydrology, 530, 306–316. https://doi.org/10.1016/j.jhydrol.2015.09.060 Wilson, R. S., Beetstra, M. A., Reutter, J. M., Hesse, G., Fussell, K. M. D. V., Johnson, L. T., King, K. W., LaBarge, G. A., Martin, J. F., & Winslow, C. (2019). Commentary: Achieving phosphorus reduction targets for Lake Erie. Journal of Great Lakes Research, 45(1), 4–11. https://doi.org/10.1016/j.jglr.2018.11.004 Wang, J. R. Frankenberger, & E. J. Kladivko. (2003). Estimating nitrate-N losses from subsurface drains using variable water sampling frequencies. Transactions of the ASAE, 46(4). https://doi.org/10.13031/2013.13965 27 CHAPTER 2 INVESTIGATION OF PHOSPHORUS TRANSPORT DYNAMICS USING HIGH- FREQUENCY MONITORING AT A SUBSURFACE-DRAINED FIELD IN THE WESTERN LAKE ERIE BASIN 28 2.1. Abstract To minimize the environmental impact of phosphorus (P) loss from subsurface-drained fields to freshwater water bodies like the Great Lakes, a detailed understanding of P transport dynamics is vital. The main objective of this study was to investigate P transport dynamics using high-frequency monitoring. We used the HydroCycle-PO4 instrument to measure total reactive P (TRP) concentration at a high resolution from a subsurface-drained farm with continuous no-till and Blount loam soil. We used a dataset containing hourly TRP concentration and hourly drainage discharge measurements in the analysis. Results showed that there was a good relationship between TRP concentration and drainage discharge (R-squared = 0.60) such that TRP had a transport- limited chemodynamic pattern, that is TRP concentration tended to increase with increase in flow during events. A 1% increase in drainage discharge resulted in a 1.36% increase in TRP load, indicating a significant increase in P concentration during high flows. We found that flow events substantially contributed to P loss (89%) because of capturing the rapid increase in P concentration during high flows. The rate of increase in P concentration during the rising limb ranged from 0.02 to 0.66 mg/L per hour. The highest 7.7% of drainage flow transported 75% of the TRP load during the monitoring period. The hysteresis pattern tended to be positive (clockwise) during the study period, indicating that preferential flow was a pathway for TRP loss. Most flow events (30 out of 36) displayed a flushing effect in which P concentration increased with rise in drainage discharge. In conclusion, high-frequency P sampling showed that management and conservation practices should target flow events to reduce P loss. 29 2.2. Introduction Nutrient transport, especially phosphorus (P), through agricultural subsurface drainage systems has been recognized as a problem for large lakes, such as Lake Erie (Ghane et al., 2016; Makarewicz et al., 2007; Pease et al., 2018; Pease et al., 2018; Williams et al., 2016). To reduce phosphorus loss in drainage discharge, management and conservation drainage practices have been recommended by environmental specialists, including 4R approaches, controlled drainage, and phosphorus removing structures (Carstensen et al., 2019; Penn et al., 2007; Hoffmann et al., 2020; King et al., 2018). If informed decisions for P reduction practices at the field scale are to be effectively implemented, a better understanding of P transport dynamics and P movement hydrochemistry from subsurface-drained farms will be needed. High flows are responsible for the majority of P loss at the watershed scale because there is a strong correlation between P concentration and drainage discharge (Bende-Michl et al., 2013; Dialameh & Ghane, 2022; Vidon & Cuadra, 2011). At the field scale, Williams et al. (2015) reported that 2% of the highest flows transported more than 50% of the annual dissolved reactive phosphorus (DRP) load based on sub-daily to daily sampling intervals. Reinhardt et al. (2005) studied a wetland receiving drainage water from agricultural fields and reported that most of the annual P load was related to high discharge events based on daily composites. Dialameh & Ghane (2022) found that P concentration increases rapidly with increase in drainage discharge based on hourly P data. However, to our knowledge, we found no study that had investigated P transport dynamics with hourly data at the field scale. Therefore, there is a need for further investigations of the patterns and controls on P transport at the field scale using fine-resolution hourly data to assess the efficacy of conservation drainage practices aimed at mitigating P loss. With the emergence of in-situ high-tech instruments, researchers have access to real-time 30 high-frequency water quality monitoring to gain a better insight of nutrient transport processes (Ouyang, 2021). Liu et al. (2020) measured drainage discharge and concentration of nitrate at high resolutions of 15 and 45 min to show the use of high-frequency monitoring. Their results showed that high-frequency monitoring can capture the rapidly changing nitrate concentration with varying drainage discharge. Speir et al. (2021) used high-frequency nitrate and drainage discharge measurements to determine the physicochemical controls of nitrate export across time scales at the watershed scale. Their results showed high-resolution nitrate concentration data was able to precisely document the magnitude of the deleterious effects of nitrate transport by flow events in agricultural watersheds. Although some studies have investigated the role of flow events on nitrate transport using high-frequency measurements, the P transport mechanism remains unknown due to technical challenges and the lack of advanced instruments. Although high-frequency data have enhanced our understanding of nitrate transport dynamics and its hydrochemical processes in drainage discharge (Liu et al., 2020), our understanding of P transport dynamics and its hydrochemical processes is still lacking due to the high cost of high-frequency sampling. To our knowledge, no study has investigated P transport dynamics and its hydrochemical processes in drainage discharge using sub-daily high-frequency data. The objectives of this field-scale study were using high-frequency P concentration measurements to: (1) analyze P transport pattern at the event scale, (2) infer likely processes involved in the patterns observed at the event scale, and (3) evaluate the effect of sampling strategy on the accuracy of P transport analysis. The value of this study is that it enhances our understanding of P transport dynamics so that informed decisions for P reduction strategies can be effectively implemented. 31 2.3. Materials and Methods 2.3.1. Study site This investigation was performed from January 2019 to July 2020 on a privately owned subsurface-drained farm (22.4 ha) in Lenawee County, Michigan, USA (Figure 2.1). The Blount loam (fine, illitic, mesic Aeric Epiaqualfs) was the field’s dominant soil type. This soil type is known as a poorly drained soil under natural conditions (USDA-NRCS, 2022). The average drain depth and spacing were 0.75 m and 12 m, respectively. During the study period, the drainage system was under conventional free drainage. Figure 2.1. Drainage layout of the study site (CL Site) in Lenawee County, Michigan, USA. 32 The cropping system was continuous corn during the study period. Cereal rye was planted as the winter cover crop. No-till was adopted as the conservation tillage practice. The farmer surface broadcasted manure at a rate of 43.2 kg-P/ha and 5.3 kg-P/ha in January and December 2019, respectively. In May 2020, pop-up commercial fertilizer was applied at planting in contact with the seed in furrows. The chemical fertilizer had the formulation of 9-18-9 with sulfur and contained 2.9 kg-P/ha. 2.3.2. Precipitation data We used an ATMOS-41 (METER Group, Inc., USA) that is a microclimate sensor suite to collect rainfall data. This built-in device measures rainfall with a resolution of 0.017 mm. Because the ATMOS-41 does not measure snow, we used the snow water equivalent data provided by National Oceanic and Atmospheric Administration (NOAA) weather station. The NOAA weather station was located at the Adrian Lenawee County Airport which was 13.5 km away from the on- farm research site. 2.3.3. High-frequency measurement of drainage discharge During the study period, hourly drainage discharge was measured by two different methods. In the first method, a sharp-crest metal-edge 45-degree V-notch weir (Agri Drain Corp., USA) was used to measure drainage discharge. The weir was inserted inside a 25-cm water-level control structure (Agri Drain Corp., USA). Two conditions should have been satisfied to use this method: (1) the flow should be less than the weir capacity (i.e., maximum weir flow); and (2) the level of water in the downstream chamber of the control structure should not rise above the V-notch apex. The water level inside the control structure was measured hourly using a HYDROS-21 water-level sensor. Then, the measured water levels and a calibrated equation developed by Shokrana and Ghane (2021) were used to estimate the hourly drainage discharge. However, when neither of 33 abovementioned two conditions were met, the second method was used. In this method, a TIENET- 350 area-velocity sensor (Teledyne ISCO, USA) installed at the downstream of the control structure was used to measure hourly drainage discharge. Overall, the area-velocity sensor and V-notch weir provided the high and low flow rates, respectively. 2.3.4. High-frequency measurement of phosphorus concentration We measured P concentration in drainage discharge at the high-resolution of two hours using the HydroCycle-PO4 instrument (Sea-Bird Scientific, USA) (Figure A2.1). The HydroCycle- PO4 is an in-situ instrument that measures real-time P concentration colorimetrically using a heteropoly molybdenum-blue complex (Murphy & Riley, 1962). The HydroCycle-PO4 conducts the P concentration measurements on an unfiltered sample. Therefore, the measurements represent TRP (Rice et al., 2017). The inorganic form of P is known as reactive because it is readily bioavailable. The Sea-Bird Scientific company does not specify an accuracy and precision range for the HydroCycle. However, recent studies have reported an accuracy ranging from –5 to 20%, standard deviation ranging from 0.0005 to 0.0020 mg/L, and coefficient of variation ranging from 0.14 to 5.78% (Johengen et al., 2017; Snazelle, 2018). The maximum P detection limit of the sensor is 1.2 mg/L while it has a minimum P detection limit of 0.002 mg/L. The concentration of TRP in drainage discharge was monitored from January 2019 to July 2020. Every week, we downloaded the TRP concentration data from the data logger and used the CycleHost software (Sea-Bird Scientific) for processing. Ultimately, to estimate hourly TRP concentrations, we linearly interpolated the 2-hour interval TRP concentrations. This linear interpolation was justified because the 2-hour concentrations captured the peak of the concentration during peak flows. Thus, linear interpolation would still capture the peak concentration. 34 The HydroCycle-PO4 was equipped with a 7.5-µm screen mesh to reduce sediment intake. The instrument was cleaned every week to maintain high-quality monitoring. Every 3 to 4 months, new cartridges were installed to maintain proper performance of the instrument. 2.3.5. Estimating the TRP load We used hourly TRP concentration and drainage discharge to estimate the hourly TRP load as Load = ∑ QiCi n i=1 (1) where Load is the TRP load (kg/ha), Qi is the hourly drainage discharge (m3/h), and Ci is the hourly TRP concentration (mg/L). A conversion factor of 4.44 × 10−5 was used to adjust the units. When any number of hours within a day was missing, we disregarded that day from the analysis. 2.3.6. Determining phosphorus transport dynamics To assess P transport dynamic, we evaluated the relationship between TRP concentration and drainage discharge (C − Q relationship) and relationship between TRP load and drainage discharge (L − Q relationship). For the C − Q relationship analysis, a linear regression was developed for the plot of natural logarithm of hourly TRP concentration, against natural logarithm of hourly drainage discharge: LnC = bLnQ + Lna (2) where C is the hourly TRP concentration (mg/L), Q is the hourly drainage discharge (mm), b is the slope, and a is a constant. The slope of the regression line provides useful information about how nutrient concentration changes with varying drainage discharge. Using thresholds defined by Bieroza et al. (2018), we classified the slopes of the regression line into three categories: 1) chemostatic indicating no significant change in P concentration when drainage discharge changes ( −0.1 < slope < 0.1 ), 2) transport-limited chemodynamic indicating an enrichment response 35 across the range of flows observed (e.g., increasing drainage discharge results in increasing TRP concentrations) (0.1 < slope), or 3) source-limited chemodynamic indicating a dilution response across the range of flows observed (e.g., increasing drainage discharge results in decreasing TRP concentrations) (slope < −0.1). For L − Q analysis, we fitted a linear regression to the plot of the natural logarithm of hourly TRP load, L (kg/ha), versus the natural logarithm of drainage discharge, which is written as LnL = bLnQ + Lna (3) a and b are a constant and slope, respectively. When the slope is equal to one, nutrient concentration remains constant under varying drainage discharge. When the slope is greater than one, high drainage discharge leads to high nutrient concentration, whereas low drainage discharge leads to low nutrient concentration. The regression slope is the percent increase in nutrient load that is induced by 1% increase in drainage discharge (Alexander et al., 1996; Tomer et al., 2003). We used the 95% confidence interval to determine whether the slope was significantly different from one. 2.3.7. Description of hypothetical sampling strategies The resolution of the data is important in the outcome of analyzing the P transport dynamics. To investigate the effect of sampling strategy on the L − Q relationship, we artificially subsampled the reference hourly dataset to create eight hypothetical time-proportional discrete sampling strategies. The sampling intervals evaluated were: 3-h, 6-h, 12-h, 1-day, 2-day, 3-day, 7-day, and 14-day intervals. For each dataset, the daily TRP load was estimated and was plotted against the daily drainage discharge. 36 2.3.8. Delineating events The flow events were delineated for event-scale analysis. The beginning of an event was the time at which there was at least 20% higher drainage discharge compared to baseflow (Liu et al., 2020) (Figure A2.2). The selection of the end point of an event is more arbitrary (Linsley et al., 1975). Compound events with more than one spike were removed from the from the dataset for this analysis. The amount of P loss during each flow event was estimated. We also determined the contributions of rising and falling limbs to P transport for each event. 2.3.9. Quantifying hysteresis and flushing indices Hysteresis analysis provides insights into spatial and temporal dynamics of nutrient source availability, nutrient storage, and hydrological pathways (Duvert et al., 2010; Williams, 1989). Normalized TRP concentrations (Eq. 4) were used for each flow event to calculate the hysteresis index (HI), which quantifies the strength and direction of the hysteresis effect. We calculated HI as follows (Lloyd et al., 2016b, 2016a): Ci,norm = Ci−Cmin Cmax−Cmin (4) HI = n i=1 ∑ (Ci,rising−Ci,falling) n (5) where Ci,norm is the normalized TRP concentration; Ci is the ith event value of measured TRP concentration; Cmin and Cmax are the event minimum and maximum TRP concentrations, respectively; HI is the event hysteresis index; Ci,rising is the normalized concentration of ith segment rising limb; Ci,falling is the normalized concentration of ith segment on the rising limb; and n is the scaling factor. A positive HI value means that the TRP concentration peak occurs before the drainage discharge peak, whereas a negative HI value means that the TRP concentration peak occurs after the drainage discharge peak. HI values between –0.1 and 0.1 display insignificant 37 hysteresis patterns (Liu et al., 2020). The flushing index quantifies event flushing/dilution effects and is calculated as described by Vaughan et al. (2017): FI = CQpeak − C0 (6) where FI is the flushing index, CQpeak is the normalized TRP concentration at the peak of the event hydrograph, and C0 is the normalized TRP concentration at the beginning of the event hydrograph. A positive FI indicates a flushing effect (chemodynamic transport-limited pattern) when TRP concentration in drainage discharge is higher during the event than in baseflow. By contrast, an event has a dilution effect (chemodynamic source-limited pattern) when FI is negative (Liu et al., 2020). The FI increases the understanding of P behavior at peak flows (Speir et al., 2021). FI values between –0.1 and 0.1 display insignificant dilution/flushing effects or chemostatic pattern (Liu et al., 2020). 2.4. Results 2.4.1. Relationship between TRP concentration and drainage discharge The hourly drainage discharge rate differed from 0.016 to 0.062 mm/h. The average hourly TRP concentration was 0.136 mg/L with minimum and maximum values of 0.007 mg/L and 1.161 mg/L, respectively. The high-frequency measurements of TRP concentration and drainage discharge showed that TRP concentration rapidly increased during flow events, while it did not fluctuate during baseflow (Figure 2.2). The linear regression between hourly TRP concentration and hourly drainage discharge was statistically significant (P-value  0.001) (Figure 2.3). 38 Figure 2.2. Precipitation, hourly TRP concentration, and drainage discharge during the study period. The green lines show the events that were investigated. Figure 2.3. Relationship between hourly TRP concentration and hourly drainage discharge during the study period (January 2019 to July 2020). 39 The TRP concentration had a transport-limited chemodynamic behavior because the slope was greater than 0.1 (b = 0.36, p < 0.001; Figure 4). According to the seasonal C − Q analysis, all seasons showed transport-limited chemodynamic pattern (spring: b = 0.58, p < 0.001; fall b = 0.43, p < 0.001; winter: b = 0.98, p < 0.001; Figure A2.3), except summer that had a chemostatic pattern (summer: b = 0.06, p > 0.05; Figure A2.3). Figure 2.4. Relationship between hourly TRP concentration and hourly drainage discharge throughout the study period. 2.4.2. Relationship between TRP load and drainage discharge The linear regression of the natural log of hourly TRP load against the natural log of hourly drainage discharge indicated that the slope of the line significantly differed from one (p < 0.001) (Figure 2.5). The slope of 1.36 means that with every 1% increase in drainage discharge TRP load increased by 1.36%. The seasonal L − Q relationship analysis also showed that the slope of the regression line was significantly different from one during spring (b = 1.58, p < 0.001), fall (b = 1.43, p < 0.001), and winter (b = 1.98, p < 0.001) while for summer, insignificant slope was observed (b = 1.06, p < 0.001) (Figure 2.6). The contribution of each flow event to P loss during the study period varied from 0.1% to 14.6%. Flow events in spring and fall had the highest 40 (4.4%) and lowest (1.1%) contribution to P loss, respectively (Figure A2.3). Figure 2.5. Relationship between hourly TRP load and hourly drainage discharge throughout the study period. Figure 2.6. Seasonal relationship between hourly TRP load and hourly drainage discharge. 41 2.4.3. Effect of sampling strategy on the accuracy of P transport dynamic analysis Our analysis showed that sampling strategy significantly affected the slope of the regression line and the coefficient of determination (Figure 2.7). For sampling strategies longer than a day, the slopes were significantly smaller than one meaning that the correlation between daily TRP load/concentration and daily drainage discharge was not strong. The coefficient of determination also significantly decreased when using low-resolution sampling strategies. Figure 2.7. Relationship between daily TRP load derived from different sampling intervals and daily drainage discharge throughout the study period. 42 2.4.4. Summary of flow events A total of 46 flow events occurred during the study period, and 36 events were selected for further analysis. The average event durations and rising limbs were 29.9 and 4.1 h, respectively (Table 2.1). The average of flow during events varied from 0.05 to 0.36 mm/h, with an average of 0.2 mm/h for all events. The hydrological and hydrogeochemical characteristics considerably differed between individual events. The average of TRP concentration during events varied from 0.14 to 0.91 mg/L. These variations in average drainage discharge and TRP concentration led to a wide range of event TRP loads, which differed from 0.001 to 0.198 kg/ha, with an average of 0.022 kg/ha for all events (Table 2.1). 2.4.5. Contribution of event flow to P loss Approximately 76% of P load was transported during the falling limb, whereas 24% of the total P was transported during the rising limb. This was due to the duration of the falling and rising limbs; the longer duration of the falling limb generates more flow, which in turn contributes to more P load (Table 2.1). The results showed that baseflow and event flow contributed approximately equally (50% each) to the cumulative total flow, although 11% and 89% of the TRP load was transported during baseflow and event flows, respectively (Figure 2.8). The higher TRP load during event flows can be explained by the rapid increase in TRP concentration with increasing drainage discharge. The top 7.7% of drainage flow transported 75% of the TRP load during the monitoring period (Figure 2.9). 43 Table 2.1. Key metrics characterizing the 38 flow events investigated in this study. Metric Flushing index, FI Hysteresis index, HI Average of flow during events (mm/h) Peak flow (mm/h) Event total water (mm) 6.5 Average of TRP concentration during events (mg/L) 0.42 Event TRP load (kg/ha) 0.036 Average Standard deviation Min Median Max 0.53 0.09 0.20 0.39 0.43 0.37 0.09 0.16 5.2 0.20 –0.69 –0.72 0.05 0.09 0.8 0.14 0.66 0.05 0.18 0.36 4.3 0.38 1.00 0.94 0.36 0.67 22.3 0.91 0.038 0.001 0.022 0.198 Event duration (h) Duration of rising limb (h) Duration of falling limb (h) 30.0 5.2 25.0 13.8 3.2 12.2 12.0 2.0 6.0 25.0 5.0 23.5 66.0 16.0 54.0 Figure 2.8. Contribution of baseflow and event flow to total flow and phosphorus load. 44 Figure 2.9. Percent of total TRP load across the range of flow conditions (as percentile). 2.4.6. Hysteresis and flushing patterns We observed 14 events with significant positive hysteresis patterns (clockwise with average HI = 0.48), 13 events with significant negative hysteresis patterns (anticlockwise with average HI = −0.30), and 9 events with no hysteresis patterns (Figure 2.10, Figure 2.11). The average and median HI values were 0.1 and 0.05, respectively (Table 2.1). The average of event duration for events with positive HI was 29.9 h, for events with negative HI was 31.2 h, and for events with no HI pattern was 28.3 h. Therefore, event duration only had a minor effect on the hysteresis pattern. Generally, events with negative HI had higher TRP concentration, TRP load, peak flow, and total water on average in comparison with events with zero or positive HI. If classifying events based on the magnitude of average peak flow, it reveals that events with negative hysteresis index had larger peak flow (0.42 mm) than events with positive hysteresis index (0.36 mm). The results for seasonal analysis showed that during the study period (Figure 2.11), events occurred in spring tended to have a positive hysteresis effect (mean HI = 0.26), events occurred in summer and fall tended to have no hysteresis with a slight tendency to positive hysteresis effect 45 (mean HI = 0.06 for summer, 0.22 for fall), and events occurred in winter tended to have a negative hysteresis effect (mean HI = -0.1). However, in all seasons, different HI values were observed suggesting a combination of distal and proximal source of P. During the study period, four farm operations (planting on 06/16/2019 and 06/02/2020 and harvesting on 06/04/2019 and 10/26/2019) were applied to the field, however, no logical relationship was observed between these operations and the hysteresis pattern of events. More data maybe need to investigate the effect of farm operation on the hysteresis and flushing indices. The FI also ranged from –0.69 to 1, with an average of 0.53. We observed 30 events (84%) with significant positive FI (FI < 0.1) or flushing effect, three events with insignificant FI value (−0.1 < FI < 0.1), and three events with significant negative FI (FI > 0.1) or dilution effect. No seasonal pattern was observed in FI values. Figure 2.10. Three different events with three different hysteresis patterns of clockwise pattern (HI>0), counter-clockwise pattern (HI<0), and no pattern (-0.1 0.1). Based on the C − Q analysis for the study period, the TRP concentration had a transport-limited chemodynamic behavior (section 2.4.1). Therefore, the transport-limited chemodynamic C − Q 50 response for individual events came to similar outcome as in the situation when aggregated data were used for the C − Q analysis. The analysis of C– Q relationship for individual events also reflects a combination of near- term controls on P concentrations, such as quantity and intensity of rainfall, which interacts with management practices, such as fertilizer application time and tillage. The distinct C– Q slopes at short versus long times scales is of imperative for accurate estimation of watershed-scale P loads (Minaudo et al., 2019), and has substantial implications for downstream water quality, given that the majority of P loss from these agricultural fields occurred during high flows. 2.5.4. L-Q Relationship According to the slope of the regression line, a 1% increase in drainage discharge resulted in a 1.36% increase in TRP load. This means that high flows led to increased TRP concentrations, which agrees with the results of our linear regression in Section 2.4.1. This result agrees with the general positive FI. The high value of the slope for high flows shows that high flows are a dominant factor for P transport. Therefore, peak flow should be targeted to reduce P loss. To our knowledge, our slope of 1.36 is the highest reported for a subsurface-drained field for any form of P. Madison et al. (2014) reported a slope of 0.96 for DRP for a subsurface-drained farm under corn-soybean rotation based on a flow-proportional composite sampling strategy, where the authors represented a multiday event with a single flow-proportional composite DRP concentration. Their results of finding a reduction in P concentration with increased flow may be due to their high-frequency sampling. Ghane et al. (2016) reported a slope of 1.11 for total phosphorus for a subsurface-drained farm under corn-soybean rotation using a flow-proportional compositing strategy, where a multiday event was represented by a single concentration. The higher slope found in our study can be explained by the need for high-frequency sampling to capture rapid 51 changes in concentration. Dialameh and Ghane (2022) reported that high-frequency time- proportional discrete P sampling is required to accurately assess the departure of the slope from one because it captures rapid variation in P concentration during high flows. Our hourly discrete TRP concentrations led to an accurate estimate of the slope of 1.36 for TRP, which reveals the importance of high flows in transporting P from subsurface-drained fields. The coefficient of determinations (R-square = 0.83) indicates that hydrology has a dominant role in TRP transport (Figure 2.4). A previous study reported that TRP transport was primarily due to drainage discharge (Tiemeyer et al., 2009). Therefore, the implementation of conservation drainage practices that reduce drainage discharge (e.g., controlled drainage) may diminish TRP load transport in drainage discharge (Evans et al., 1995). 2.5.5. Hysteresis and flushing indices The variation in hysteresis index suggests that multiple flow pathways and transport mechanisms contribute to P loss to drain pipes (Williams et al., 2018). The general hysteresis pattern during our study period tended to be positive, suggesting dominance of proximal source of P, or rapid source mobilization over the flow event progress. Liu et al. (2020) found an average negative HI value for nitrate (based on a 45-minute time-proportional sampling strategy), which means that nitrate concentration peaked after the drainage discharge peak. Ulén and Persson (1999) reported equal frequencies of clockwise and anticlockwise HI patterns for DRP in a subsurface-drained field without quantifying HI values based on a flow-proportional composite sampling strategy of 10 consecutive subsamples corresponding to 0.2 mm of drainage discharge. To our knowledge, no study has quantified the HI values for the relationship between P and drainage discharge using high- frequency hourly data. In 14 events with a significant positive hysteresis (clockwise pattern), the TRP concentration 52 peak occurred during the rising limb of the hydrograph and P was transported rapidly to the sample collection point. Positive hysteresis indicates that P rapidly moves from the soil surface to the drainpipe via preferential (macropore) flow pathways without passing through the soil matrix (Williams et al., 2016). Manure was surface broadcasted in 2019 at the study site. The combination of surface-broadcast manure and no-till promotes preferential loss of P through macropores (Jarvie et al., 2017). Overall, the positive hysteresis observed at our study site indicated that preferential flow pathways contributed to P loss in a no-till soil with loam soil. In fields with a positive hysteresis index, macropores are likely the primary pathway for rapid movement of P from the soil surface to the drain pipe (Figure 2.12). Heavy clay soil can promote macropores due to the swelling and shrinking of the clay minerals (Jarvis, 2020; King et al., 2015). Cullum (2009) found that no-till soil developed more macropores than tilled soil. Macrae et al. (2021) recommended that no-till was not combined with surface-broadcast fertilization. Instead, fertilizer should be subsurface injected in a no-till field (Williams et al., 2018). If surface broadcast cannot be avoided, the fertilizer should be incorporated into the soil, such as subsurface banding (Macrae et al., 2021; Williams et al., 2016). Figure 2.12. Schematic view of preferential and matrix flows from the surface to the drainpipes. 53 In the 13 events with a significant negative hysteresis (anticlockwise), the TRP concentration peak occurred after the falling limb of the hydrograph. This may happen when the P source is distant, or P moves through paths that are activated later during the flow event (i.e., slow source mobilization) (Speir et al., 2021). The general positive FI value agrees with the slope of the regression line (1.36) discussed in Section 2.4.2, which showed a general P concentration increase during high-flow events. During the study period, 30 of 36 events displayed flushing effects (FI > 0.1) (Figure 2.11). A positive FI value shows that the P concentration during the flow peak was higher than that in the baseflow immediately before the start of the flow event, indicating that P concentration increases with drainage discharge. Three events displayed dilution effects (FI < 0), indicating that P concentration decreases with increasing discharge. The flushing or dilution effect of an event flow depends on P availability in the soil profile compared to the P concentration in pre-event baseflow (Bieroza et al., 2018; Zhang, 2018). The combined results of FI and HI show that events had varying hysteretic behaviors, with 35.3%, 31.6%, and 8.8% exhibiting anticlockwise with flushing, clockwise with flushing, and clockwise with dilution, respectively. These different concentration-discharge responses indicate that, apart from drainage discharge, the status of soil P, and rapid exchanges between P pools, the event magnitude and duration also regulated solute delivery (Welikhe et al., 2021). 2.5.6. Effect of sampling strategy on the accuracy of P transport dynamic analysis The analysis of sampling strategies showed that the longer the sampling interval, the lower the slope of the regression line and the coefficient of determination (Figure 2.7). Since low- resolution sampling strategies (generally longer than one day) are not able to capture the rapid fluctuations of P concentration in the drainage discharge, especially during high flows when P 54 concentration rapidly increases and recedes, thereby likely underestimating the slope of the regression line and the coefficient of determination (Figure A2.4). Underestimation in estimating these values leads to undermining the role of hydrology especially high flows in P transport. For example, if adopting a time-proportional 7-day discrete sampling strategy (grab sampling once every week), the slope of the regression line would have reduced from 1.37 (hourly) to 0.89, and the coefficient of determination would have reduced from 0.83 (hourly) to 0.48 (Figure 2.7). 2.6. Conclusions We assessed the P transport dynamics and investigated the hysteresis patterns of P concentration based on high-frequency monitoring of TRP concentrations. Our study resulted in the following key conclusions. • TRP concentration increased rapidly as drainage discharge increased. There was a good linear relationship between TRP concentration and drainage discharge (R-square = 0.72). • A 1% increase in drainage discharge resulted in a 1.36% increase in TRP load, indicating a significant increase in TRP concentration at high flows. This shows the importance of peak flow in transporting P from subsurface-drained fields. • High-resolution sampling strategies are important to obtain accurate results due to dynamic behavior of P concentration in drainage discharge. When evaluating the contribution of event flow to P loss, high-frequency monitoring of P in drainage discharge is needed. • Baseflow and event flow contributed to 22% and 78% of the TRP load, respectively, showing that event flows predominantly contribute to P loss. • The top 7.8% of drainage flow transported 75% of the TRP load during the monitoring period. • Hysteresis tended to be positive during the study period, indicating that P was transported 55 rapidly through preferential flow pathways and reached the sample collection point before peak flow. In conclusion, our high-frequency monitoring of P provided new insights into the considerable contribution of event flows to P loss (78%), which had not been previously reported to such a large extent. Our study captured the large contribution of event flows because of the presence of high-frequency hourly TRP sampling, which accurately measured the rapid increase in P concentration during high flows. Generally, high resolution strategies shorter than 6 h are required to capture rapid variation in P concentration during high flow. Given the substantial contribution of event flows to P loss, management and conservation practices should target flow events and reduce the peak discharge as one of the strategies for reducing P loss. 56 REFERENCES Alexander, R. B., Murdoch, P. S., & Smith, R. A. (1996). Streamflow-induced variations in nitrate flux in tributaries to the Atlantic Coastal Zone. Biogeochemistry, 33(3), 149–177. https://doi.org/10.1007/BF02181070 Bende-Michl, U., Verburg, K., & Cresswell, H. P. (2013). High-frequency nutrient monitoring to infer seasonal patterns in catchment source availability, mobilisation and delivery. Environmental Monitoring and Assessment, 185(11), 9191–9219. https://doi.org/10.1007/s10661-013-3246-8 Bieroza, M. Z., Heathwaite, A. L., Bechmann, M., Kyllmar, K., & Jordan, P. (2018a). The concentration-discharge slope as a tool for water quality management. Science of The Total Environment, 630, 738–749. https://doi.org/10.1016/j.scitotenv.2018.02.256 Bieroza, M. Z., Heathwaite, A. L., Bechmann, M., Kyllmar, K., & Jordan, P. (2018b). The concentration-discharge slope as a tool for water quality management. Science of The Total Environment, 630, 738–749. https://doi.org/10.1016/j.scitotenv.2018.02.256 Carstensen, M. v., Børgesen, C. D., Ovesen, N. B., Poulsen, J. R., Hvid, S. K., & Kronvang, B. (2019). Controlled Drainage as a Targeted Mitigation Measure for Nitrogen and Phosphorus. Journal of Environmental Quality, 48(3), 677–685. https://doi.org/10.2134/jeq2018.11.0393 Chad J. Penn, Ray B. Bryant, Peter J. A. Kleinman, & Arthur L. Allen. (2007). Removing dissolved phosphorus from drainage ditch water with phosphorus sorbing materials. Journal of Soil and Water Conservation, 62(4), 269–276. Cullum, R. F. (2009). Macropore flow estimations under no-till and till systems. CATENA, 78(1), 87–91. https://doi.org/10.1016/j.catena.2009.03.004 Daigh, A. L. M., Zhou, X., Helmers, M. J., Pederson, C. H., Horton, R., Jarchow, M., & Liebman, M. (2015). Subsurface Drainage Nitrate and Total Reactive Phosphorus Losses in Bioenergy-Based Prairies and Corn Systems. Journal of Environmental Quality, 44(5), 1638–1646. https://doi.org/10.2134/jeq2015.02.0080 Daly, K., Tuohy, P., Peyton, D., Wall, D. P., & Fenton, O. (2017). Field soil and ditch sediment phosphorus dynamics from two artificially drained fields on poorly drained soils. Agricultural Water Management, 192, 115–125. https://doi.org/10.1016/j.agwat.2017.07.005 Dialameh, B., & Ghane, E. (2022). Effect of water sampling strategies on the uncertainty of phosphorus load estimation in subsurface drainage discharge. Journal of Environmental Quality. 57 Djodjic, F., Ulen, B., & Bergstrom, L. (2000). Temporal and spatial variations of phosphorus losses and drainage in a structured clay soil. Water Research, 34(5), 1687–1695. https://doi.org/10.1016/S0043-1354(99)00312-7 Djodjic, F., Ulén, B., & Bergström, L. (2000). Temporal and spatial variations of phosphorus losses and drainage in a structured clay soil. Water Research, 34(5), 1687–1695. https://doi.org/10.1016/S0043-1354(99)00312-7 Duvert, C., Gratiot, N., Evrard, O., Navratil, O., Némery, J., Prat, C., & Esteves, M. (2010). Drivers of erosion and suspended sediment transport in three headwater catchments of the Mexican Central Highlands. Geomorphology, 123(3–4), 243–256. https://doi.org/10.1016/j.geomorph.2010.07.016 Evans, R. O., Skaggs, R. W., & Giiliam, J. W. (1995). Controlled versus conventional drainage effects on water quality. Journal of Irrigation and Drainage Engineering, 121(4), 271– 276. Ghane, E., Ranaivoson, A. Z., Feyereisen, G. W., Rosen, C. J., & Moncrief, J. F. (2016a). Comparison of contaminant transport in agricultural drainage water and urban stormwater runoff. Plos One, 11(12), e0167834. https://doi.org/10.1371/journal.pone.0167834 Ghane, E., Ranaivoson, A. Z., Feyereisen, G. W., Rosen, C. J., & Moncrief, J. F. (2016b). Comparison of Contaminant Transport in Agricultural Drainage Water and Urban Stormwater Runoff. PLOS ONE, 11(12), e0167834. https://doi.org/10.1371/journal.pone.0167834 Hoffmann, C. C., Zak, D., Kronvang, B., Kjaergaard, C., Carstensen, M. V., & Audet, J. (2020). An overview of nutrient transport mitigation measures for improvement of water quality in Denmark. Ecological Engineering, 155(April), 105863. https://doi.org/10.1016/j.ecoleng.2020.105863 Jarvie, H. P., Johnson, L. T., Sharpley, A. N., Smith, D. R., Baker, D. B., Bruulsema, T. W., & Confesor, R. (2017). Increased Soluble Phosphorus Loads to Lake Erie: Unintended Consequences of Conservation Practices? Journal of Environmental Quality, 46(1), 123– 132. https://doi.org/10.2134/jeq2016.07.0248 Jarvis, N. J. (2020). A review of non-equilibrium water flow and solute transport in soil macropores: principles, controlling factors and consequences for water quality. European Journal of Soil Science, 71(3), 279–302. https://doi.org/10.1111/ejss.12973 Johengen, T., Purcell, H., Tamburri, M., Loewensteiner, D., Smith, G. J., Schar, D., McManus, M., & Walker, G. (2017). Performance verification statement for Sea-bird Scientific HydroCycle phosphate analyzer. In Allicance for Coastal Technologies. https://doi.org/http://dx.doi.org/10.25607/OBP-289 58 King, K. W., Williams, M. R., LaBarge, G. A., Smith, D. R., Reutter, J. M., Duncan, E. W., & Pease, L. A. (2018). Addressing agricultural phosphorus loss in artificially drained landscapes with 4R nutrient management practices. Journal of Soil and Water Conservation, 73(1), 35–47. https://doi.org/10.2489/jswc.73.1.35 King, K. W., Williams, M. R., Macrae, M. L., Fausey, N. R., Frankenberger, J., Smith, D. R., Kleinman, P. J. A., & Brown, L. C. (2015). Phosphorus transport in agricultural subsurface drainage: A review. Journal of Environmental Quality, 44(2), 467–485. https://doi.org/10.2134/jeq2014.04.0163 Linsley, R. K., Kohler, M. A., & Paulhus, J. L. H. (1975). Hydrology for engineers. McGraw-Hill. Liu, W., Youssef, M. A., Birgand, F. P., Chescheir, G. M., Tian, S., & Maxwell, B. M. (2020). Processes and mechanisms controlling nitrate dynamics in an artificially drained field: Insights from high-frequency water quality measurements. Agricultural Water Management, 232, 106032. https://doi.org/10.1016/j.agwat.2020.106032 Lloyd, C. E. M., Freer, J. E., Johnes, P. J., & Collins, A. L. (2016a). Technical Note: Testing an improved index for analysing storm discharge-concentration hysteresis. Hydrology and Earth System Sciences, 20(2), 625–632. https://doi.org/10.5194/hess-20-625-2016 Lloyd, C. E. M., Freer, J. E., Johnes, P. J., & Collins, A. L. (2016b). Using hysteresis analysis of high-resolution water quality monitoring data, including uncertainty, to infer controls on nutrient and sediment transfer in catchments. Science of the Total Environment, 543, 388– 404. https://doi.org/10.1016/j.scitotenv.2015.11.028 Macrae, M., Jarvie, H., Brouwer, R., Gunn, G., Reid, K., Joosse, P., King, K., Kleinman, P., Smith, D., Williams, M., & Zwonitzer, M. (2021). One size does not fit all: Toward regional conservation practice guidance to reduce phosphorus loss risk in the Lake Erie watershed. Journal of Environmental Quality, 50(3), 529–546. https://doi.org/10.1002/jeq2.20218 Macrae, M. L., English, M. C., Schiff, S. L., & Stone, M. (2007a). Intra-annual variability in the contribution of tile drains to basin discharge and phosphorus export in a first-order agricultural catchment. Agricultural Water Management, 92(3), 171–182. https://doi.org/10.1016/j.agwat.2007.05.015 Macrae, M. L., English, M. C., Schiff, S. L., & Stone, M. (2007b). Intra-annual variability in the contribution of tile drains to basin discharge and phosphorus export in a first-order agricultural catchment. Agricultural Water Management, 92(3), 171–182. https://doi.org/10.1016/j.agwat.2007.05.015 Madison, A. M., Ruark, M. D., Stuntebeck, T. D., Komiskey, M. J., Good, L. W., Drummy, N., & Cooley, E. T. (2014). Characterizing phosphorus dynamics in tile-drained agricultural fields of eastern Wisconsin. Journal of Hydrology, 519, 892–901. https://doi.org/10.1016/j.jhydrol.2014.08.016 59 Makarewicz, J., D’Aiuto, P., & Bosch, I. (2007). Elevated Nutrient Levels from Agriculturally Dominated Watersheds Stimulate Metaphyton Growth. Journal of Great Lakes Research, 33(2), 437–448. Minaudo, C., Dupas, R., Gascuel-Odoux, C., Roubeix, V., Danis, P.-A., & Moatar, F. (2019). Seasonal and event-based concentration-discharge relationships to identify catchment controls on nutrient export regimes. Advances in Water Resources, 131, 103379. https://doi.org/10.1016/j.advwatres.2019.103379 Murphy, J., & Riley, J. P. (1962). A modified single solution method for the determination of phosphate in natural waters. Analytica Chimica Acta, 27, 31–36. https://doi.org/10.1016/S0003-2670(00)88444-5 Nazari, S., Ford, W. I., & King, K. W. (2021). Quantifying hydrologic pathway and source connectivity dynamics in tile drainage: Implications for phosphorus concentrations. Vadose Zone Journal, 20(5). https://doi.org/10.1002/vzj2.20154 Ouyang, Y. (2021). A flow-weighted approach to generate daily total phosphorus loads in streams based on seasonal loads. Environmental Monitoring and Assessment, 193(7), 422. https://doi.org/10.1007/s10661-021-09199-4 Pease, L. A., Fausey, N. R., Martin, J. F., & Brown, L. C. (2018). Weather, Landscape, and Management Effects on Nitrate and Soluble Phosphorus Concentrations in Subsurface Drainage in the Western Lake Erie Basin. Transactions of the ASABE, 61(1), 223–232. https://doi.org/10.13031/trans.12287 Pease, L. A., King, K. W., Williams, M. R., LaBarge, G. A., Duncan, E. W., & Fausey, N. R. (2018). Phosphorus export from artificially drained fields across the Eastern Corn Belt. Journal of Great Lakes Research, 44(1), 43–53. https://doi.org/10.1016/j.jglr.2017.11.009 Reinhardt, M., Gächter, R., Wehrli, B., & Müller, B. (2005). Phosphorus Retention in Small Constructed Wetlands Treating Agricultural Drainage Water. Journal of Environmental Quality, 34(4), 1251–1259. https://doi.org/10.2134/jeq2004.0325 Rice, E. W., Baird, R. B., & Eaton, A. D. (2017). Methods for the examination of water and wastewater. American Public Health Association, American Water Works, 23. Shokrana, M. B. S., & Ghane, E. (2021). An empirical V-notch weir equation and standard procedure to accurately estimate drainage discharge. Applied Engineering in Agriculture. Applied Engineering in Agriculture, 37(6), 1097–1105. https://doi.org/https://doi.org/10.13031/aea.14617 Snazelle, T. (2018). Laboratory evaluation of the Sea-Bird Scientific HydroCycle-PO4 phosphate pensor. https://pubs.usgs.gov/of/2018/1120/ofr20181120.pdf 60 Soil Survey Staff, Natural Resources Conservation Service, United States Department of Agriculture. Web Soil Survey. Available online at the following link: http://websoilsurvey.sc.egov.usda.gov/. Accessed [09/20/2022]. Speir, S. L., Tank, J. L., Bieroza, M., Mahl, U. H., & Royer, T. v. (2021). Storm size and hydrologic modification influence nitrate mobilization and transport in agricultural watersheds. Biogeochemistry, 156(3), 319–334. https://doi.org/10.1007/s10533-021- 00847-y Stamm, C., Flühler, H., Gächter, R., Leuenberger, J., & Wunderli, H. (1998). Preferential Transport of Phosphorus in Drained Grassland Soils. Journal of Environmental Quality, 27(3), 515–522. https://doi.org/10.2134/jeq1998.00472425002700030006x Tiemeyer, B., Kahle, P., & Lennartz, B. (2009). Phosphorus losses from an artificially drained rural lowland catchment in North-Eastern Germany. Agricultural Water Management, 96(4), 677–690. https://doi.org/10.1016/j.agwat.2008.10.004 Tomer, M. D., Meek, D. W., Jaynes, D. B., & Hatfield, J. L. (2003). Evaluation of nitrate nitrogen fluxes from a tile-drained watershed in central Iowa. Journal of Environmental Quality, 32, 642–653. https://doi.org/10.2134/jeq2003.6420 Ulén, B., & Persson, K. (1999). Field-scale phosphorus losses from a drained clay soil in Sweden. Hydrological Processes, 13(17), 2801–2812. https://doi.org/10.1002/(SICI)1099- 1085(19991215)13:17<2801::AID-HYP900>3.0.CO;2-G Vaughan, M. C. H., Bowden, W. B., Shanley, J. B., Vermilyea, A., Sleeper, R., Gold, A. J., Pradhanang, S. M., Inamdar, S. P., Levia, D. F., Andres, A. S., Birgand, F., & Schroth, A. W. (2017). High‐frequency dissolved organic carbon and nitrate measurements reveal differences in storm hysteresis and loading in relation to land cover and seasonality. Water Resources Research, 53(7), 5345–5363. https://doi.org/10.1002/2017WR020491 Vidon, P., & Cuadra, P. E. (2010). Impact of precipitation characteristics on soil hydrology in tile- drained landscapes. Hydrological Processes, 24(13), 1821–1833. https://doi.org/10.1002/hyp.7627 Vidon, P., & Cuadra, P. E. (2011). Phosphorus dynamics in tile-drain flow during storms in the US Midwest. Agricultural Water Management, 98(4), 532–540. https://doi.org/10.1016/j.agwat.2010.09.010 Welikhe, P., Brouder, S. M., Volenec, J. J., Gitau, M., & Turco, R. F. (2021). Dynamics of dissolved reactive phosphorus loss from phosphorus source and sink soils in tile-drained systems. Journal of Soil and Water Conservation, 00012. https://doi.org/10.2489/jswc.2022.00012 Williams, G. P. (1989). Sediment concentration versus water discharge during single hydrologic events in rivers. Journal of Hydrology, 111(1–4), 89–106. https://doi.org/10.1016/0022- 61 1694(89)90254-0 Williams, M. R., King, K. W., Baker, D. B., Johnson, L. T., Smith, D. R., & Fausey, N. R. (2016). Hydrologic and biogeochemical controls on phosphorus export from Western Lake Erie tributaries. Journal of Great Lakes Research, 42(6), 1403–1411. https://doi.org/10.1016/j.jglr.2016.09.009 Williams, M. R., King, K. W., Duncan, E. W., Pease, L. A., & Penn, C. J. (2018). Fertilizer placement and tillage effects on phosphorus concentration in leachate from fine-textured soils. Soil and Tillage Research, 178, 130–138. https://doi.org/10.1016/j.still.2017.12.010 Williams, M. R., King, K. W., Ford, W., Buda, A. R., & Kennedy, C. D. (2016a). Effect of tillage on macropore flow and phosphorus transport to tile drains. Water Resources Management, 52(4), 2868–2882. https://doi.org/10.1002/2015WR017650.Received Williams, M. R., King, K. W., Ford, W., Buda, A. R., & Kennedy, C. D. (2016b). Effect of tillage on macropore flow and phosphorus transport to tile drains. Water Resources Research, 52(4), 2868–2882. https://doi.org/10.1002/2015WR017650 Williams, M. R., King, K. W., Macrae, M. L., Ford, W., van Esbroeck, C., Brunke, R. I., English, M. C., & Schiff, S. L. (2015a). Uncertainty in nutrient loads from tile-drained landscapes: Effect of sampling frequency, calculation algorithm, and compositing strategy. Journal of Hydrology, 530, 306–316. https://doi.org/10.1016/j.jhydrol.2015.09.060 Williams, M. R., King, K. W., Macrae, M. L., Ford, W., van Esbroeck, C., Brunke, R. I., English, M. C., & Schiff, S. L. (2015b). Uncertainty in nutrient loads from tile-drained landscapes: Effect of sampling frequency, calculation algorithm, and compositing strategy. Journal of Hydrology, 530, 306–316. https://doi.org/10.1016/j.jhydrol.2015.09.060 Williams, M. R., Livingston, S. J., Penn, C. J., Smith, D. R., King, K. W., & Huang, C. (2018). Controls of event-based nutrient transport within nested headwater agricultural watersheds of the western Lake Erie basin. Journal of Hydrology, 559, 749–761. https://doi.org/10.1016/j.jhydrol.2018.02.079 Zhang, Q. (2018). Synthesis of nutrient and sediment export patterns in the Chesapeake Bay watershed: Complex and non-stationary concentration-discharge relationships. Science of The Total Environment, 618, 1268–1283. https://doi.org/10.1016/j.scitotenv.2017.09.221 62 CHAPTER 3 SUBSURFACE DRAINAGE DESIGN AND SUBIRRIGATION AS CLIMATE-SMART STRATEGIES FOR RESILIENT CROP PRODUCTION 63 3.1. Abstract To build a more resilient crop production system in a changing climate, it is crucial to understand how future weather patterns affect subsurface drainage design and whether subirrigation will be needed in the future for crop production. Therefore, the main objectives of this study were to 1) evaluate the change in the optimum drain spacing from using historical (1994 - 2023) and future (2030 - 2059) weather data, and 2) evaluate the efficacy of subirrigation to alleviate yield loss due to drought stress in southeast Michigan, USA. A total of 27 general circulation models with a moderate greenhouse gas emission scenario (shared socioeconomic pathway 2-4.5) were used for climate projections. Simulations were performed using the DRAINMOD model, and the optimum drain spacing was determined based on the maximum average annual return on investment. Results showed that the projected 30-year average annual precipitation is not expected to change significantly while that of the temperature will increase by 2.5°C in the future. Future optimum drain spacings for depths of 75 cm and 125 cm were found to be 300 cm and 600 cm wider than historical spacings, respectively. On average, there was a 23% decrease in 30-year average annual drainage discharge, attributed to an average 17% increase in evapotranspiration. Drought stress is projected to be the primary cause of yield loss in the future, due to increased temperatures and an average 8% deeper water-table depth. Subirrigation shows high potential in reducing year-to-year crop yield variability in the future (decreasing the coefficient of variation for the yield from 0.26 to 0.06, on average) and increasing yield by up to 31%. In the past, subirrigation initiation was feasible in late June with a weir depth of 70 cm. However, in the future, subirrigation is anticipated to be more advantageous when starting sooner in early to mid-June, coupled with a shallower weir depth ranging from 65 cm to 55 cm. In conclusion, a wider drain spacing, providing reduced drainage intensity, along with subirrigation 64 may be needed in the future to mitigate crop yield loss from drought stress. 3.2. Introduction The Midwest USA is recognized as the primary corn and soybeans producer, mainly because of its fertile soil and substantial rainfall. Due to inadequate natural drainage in many Midwestern soils, the implementation of subsurface drainage has significantly enhanced crop yields (Skaggs et al., 1994). Subsurface drainage systems mitigate crop wet stress by lowering the water table and facilitating timely field operations (Arnold, 2004; Willison et al., 2021). However, to maximize the efficacy of a subsurface drainage systems, a proper drainage design (drain depth and spacing) is essential. When designing drainage systems, factors such as the type of crops, soil characteristics, pipe perforation style, and climate should be taken into consideration (Moustafa et al., 1987; Skaggs & Nassehzadeh-Tabrizi, 1986). Choosing the drain depth depends on the soil characteristics and primarily involves installing the pipe in a permeable soil layer away from tillage field operations and away from the restrictive layer. Determining the drain spacing involves the use of the steady-state Hooghoudt equation. Nevertheless, the drain depth and drain spacing are interrelated; shallow drains need narrower spacing to achieve equivalent drainage. Application of the Hooghoudt equation in drainage design has been facilitated by decision-support tools such as the Drain Spacing Tool that uses the site-specific soil and growing season rainfall to estimate the optimum drain spacing (Ghane, 2023). However, these tools mostly use historical weather data to estimate the drain spacing while climate scientists are forecasting changes in precipitation patterns in the future. When using such tools to determine the optimum drain spacing, it is important to consider the impact of climate change on the optimum spacing. Although subsurface drainage plays a vital role in the agricultural water balance in the 65 Midwest, our comprehension of how future climate change may impact the weather parameters (temperature and the quantity and seasonal patterns of precipitation), hydrology (drainage discharge, WTD, and evapotranspiration), drainage design (drain depth and spacing) and crop yield in drained fields remains limited. This understanding is essential to determine future strategies aimed at safeguarding agricultural production while maintaining environmental quality. Southworth et al. (2000) projected that corn yields could decrease by as much as 30% between 2050 and 2059 compared to current yields. Similarly, Petersen (2019) estimated a reduction of around 100 bushels per acre in corn yields for the southeast Michigan and northwest Ohio region from 2090 to 2100 relative to present levels. Deryng et al. (2014) found that the global average corn yield will decrease, ranging from −2.9 ± 2.6% under RCP 2.6 to −12.8 ± 6.7% under RCP 8.5 by the 2080s. Anticipated future climate change also is projected to vary geographically across regions (Pease et al., 2017). According to Intergovernmental Panel on Climate Change (2013), variations in the water balance will be a function of local shifts in temperature and precipitation pattern. Pryor et al. (2014) predicted that by midcentury, Ohio will experience a more pronounced rise in temperature and the frequency of days with heavy precipitation would be substantially different compared to Iowa. Nevertheless, temperature fluctuations and variations in rainfall patterns between Ohio and Iowa might result in different effects on drainage discharge in these respective states. Hence, it is crucial to evaluate the water balance using localized climate projections in order to assess the effect of climate change on the hydrology of a drained field and its drainage design within the Western Lake Erie Basin. Interest in supplemental irrigation, like subirrigation, is increasing in the Midwest due to more frequent summer droughts (Galloway et al., 2014). The region's extensive subsurface 66 drainage systems offer potential for subirrigation adoption, although it remains underutilized (Yu et al., 2020). Subirrigation, by applying water beneath the soil surface, helps maintain optimal water table depth and can utilize existing drainage systems for both drainage and irrigation purposes (Cooper et al., 1999) and (Zucker & Brown, 1998). Research investigating the impact of subirrigation on crop yield consistently demonstrates a significant increase in yield, which tends to stabilize at elevated levels. Investigating the efficacy of subirrigation in a field with Omulga silt loam in southern Ohio, Fisher et al. (1999) reported a 19% increase in corn yield with adoption of subirrigation. Cooper et al. (1999) observed a notable rise in corn production, up to 30%, on Ravenna silt loam and Hoytville silty clay loam soils in Ohio, particularly during dry years. Furthermore, in sandy loam soil in southwestern Ontario, Ng et al. (2002) observed a 64% larger corn yield with subirrigation, while Mejia et al. (2000) found increases ranging from 2.8% to 13.8% in eastern Ontario. Gunn et al. (2018) investigated the effectiveness of subirrigation in increasing corn yield using the historical (from 1984 through 2013) and future (from 2041 through 2070) climate data in fields with different soil series located in northwest Ohio. They reported as high as 26.5% and 36% increase in corn yield in the past and future, respectively, with subirrigation. Conversely, in Woodslee, Ontario, Drury et al. (2009) noted a significant decrease in corn yield under subirrigation practice on Brookston clay loam soil, attributing this outcome to substantial precipitation in August and drainage design. Expanding upon field measurements, modeling studies are necessary to forecast the potential effect of subirrigation under future climate scenarios. The efficacy of subirrigation in improving the crop yield has been investigated using the DRAINMOD model. DRAINMOD is a computer model for field hydrology water balance, capable of simulating various drainage practices, including free subsurface drainage, controlled drainage, and subirrigation, either 67 individually or in combination (Skaggs et al., 2012). In a study on Rains and Portsmouth sandy loam soils in North Carolina, Murugaboopathi et al. (1995) used uncalibrated DRAINMOD simulations over a 37-year period (1950–1986). Their findings revealed a 21% relative yield advantage for corn under subirrigation compared to free subsurface drainage. Gunn et al. (2018) used DRAINMOD to investigate the effectiveness of subirrigation in increasing corn yield using the historical (from 1984 through 2013) and future (from 2041 through 2070) climate data in fields with different soil series located in northwest Ohio. They reported as high as 26.5% and 36% increase in corn yield in the past and future, respectively, with subirrigation. Although climate change is already affecting the hydrology of subsurface drained fields and crop yield, the literature highlights a significant gap in understanding regarding the potential impact of climate change on drainage design and the performance of subirrigation practice. While Ghane et al. (2021) and Ghane and Askar (2021) have investigated the impact of climate change on the drainage design, their focus has been limited to a 30-year historical period from 1990 to 2019, not the future. Also, regarding subirrigation, we only found Gunn et al. (2018) investigating the climate change impact on relative crop yield in subirrigated fields in northwest Ohio, however, they only used three General Circulation Models (GCM) which indicates a need for further investigation in this area. No study also has investigated how subirrigation management, such as weir depth and irrigation start time, might potentially change from the past to the future. Therefore, the main objectives of this study were 1) evaluating the impact of climate change on weather parameters (temperature and precipitation) and hydrology (drainage discharge, water-table depth, and evapotranspiration) of a drained field, 2) assessing the impact of climate change on drainage design (drain depth and drain spacing), and 3) investigating the efficacy of subirrigation to alleviate the adverse impact of climate change in southeast Michigan, USA. The 68 significance of this work lies in its ability to offer insights into how the historical optimum drainage design responds to changing climate conditions, leveraging the latest CMIP6 modeling framework. This study also enhances our understanding of subirrigation’s performance in the future so that informed policies to improve crop yield can be effectively implemented. In general, results from this study help us to enhance the resiliency of subsurface drained fields and improve their economic viability as well as environmental sustainability. 3.3. Materials and Methods 3.3.1. Site description This study was performed on a private farm situated in Blissfield, Lenawee County, southeast Michigan, USA (Figure 3.1) from October 2018 to December 2023. The drainage area was 7.6 ha (18.8 ac) with a field grade of 0.1%. Throughout the entire duration of the study, the drainage system performed under conventional free drainage. The cropping system was corn- soybean rotation and vertical tillage was performed prior to planting corn. Figure 3.1. The Blissfield site is situated within the hydrological boundary of the River Raisin Watershed. This watershed ultimately drains into the western Lake Erie basin at Monroe Harbor. 69 3.3.2. Measurement of drainage discharge and water-table depth The drainage discharge and water-table depth (WTD) midway between the lateral drains were measured to calibrate and validate the DRAINMOD model. The WTD between the laterals was measured hourly in an observation well with a depth of 152cm (5 ft) using a HYDROS-21 water-depth sensor. We utilized two different approaches to measure hourly drainage discharge. We employed a metal-edge sharp-crest 45° V-notch weir, set up within a 25-cm water-level control structure, both from Agri Drain Corp. We used this approach only under specific circumstances: first, when water remained within the V-notch, ensuring it did not surpass the capacity of the weir, and second, when the water level in the downstream of the structure did not exceed the apex height of the V- notch. We installed a HYDROS-21 water-depth sensor from METER Group inside the upstream chamber of the control structure to measure the water level every hour. Subsequently, we used a calibrated for the V-notch weir equation to determine the hourly drainage discharge based on the water-level measurements. We also utilized a TIENET-350 area-velocity sensor from Teledyne ISCO, placed inside a pipe downstream of the control structure. This approach was used when either water overflowed the V-notch weir, surpassing its capacity, or when the water level in the downstream chamber exceeded the height of the V-notch apex. The area-velocity sensor and the V-notch weir were used for high flow and low flow rates, respectively. 70 3.3.3. The DRAINMOD model setup to simulate hydrology, relative crop yield, and subirrigation We utilized DRAINMOD to predict drainage discharge, WTD, crop yield, and evapotranspiration (ET) for the past and future. The DRAINMOD model is a processed-based field-scale one-dimensional model widely employed for assessing the hydrology of fields with subsurface drainage systems (Skaggs, 1978, 1982; Skaggs et al., 2012). DRAINMOD utilizes water-balance equations to simulate hydrological processes with soil, weather, crop, and drainage system parameters as the primary inputs. Besides hydrological simulations, the DRAINMOD model estimates crop yield based on four stress factors as well: drought stress, wet stress, salinity stress, and stress attributed to planting delays due to poor trafficability conditions (Evans et al., 1991). The DRAINMOD’s performance in simulating the hydrology and crop has been tested by researchers across the United States and globally, as documented in studies by (Adhikari et al., 2020; Evans et al., 1991; Gunn et al., 2018; Kanwar et al., 1994; Malakshahi et al., 2020; Skaggs et al., 2012; Sojka et al., 2020; Wang et al., 2006; Youssef et al., 2018). DRAINMOD employs a stress day index method to predict the relative crop yield responses to different soil water conditions (Evans et al., 1991). This feature enables DRAINMOD to assess the efficacy of different water management practices and strategies by quantifying relative crop yield responses. In DRAINMOD, the definition of relative crop yield is the ratio of the yield under specific stress conditions such as wet, drought, or planting delay to the highest potential yield achievable. DRAINMOD derives relative crop yield by assessing predicted stresses, including the quantity and continuation of excess water beyond the water table threshold defined by the user, the magnitude and extent of drought, and the duration of planting delays. Equation (1) is utilized to estimate the overall relative crop yield during a simulated year. 71 𝑌𝑅 = 𝑌𝑅𝑤 × 𝑌𝑅𝑑 × 𝑌𝑅𝑝 (1) where 𝑌𝑅𝑤 is the relative crop yield including excess water stress, 𝑌𝑅𝑑 is the relative crop yield including drought stress, and 𝑌𝑅𝑝 is the relative crop yield including planting delay stress. Predictions of these stresses are derived from linear functions, which consider the severity of both wet and drought stresses throughout the growing season (GS), as well as any planting delays. Previous research has demonstrated that DRAINMOD's predictions of crop relative yield closely align with estimates derived from field measurements, contingent upon the prevailing conditions during growing season, as evidenced by studies conducted by (Ale et al., 2009; Kanwar et al., 1994; Satchithanantham & Ranjan, 2015; Wang et al., 2006). 3.3.3.1. Drainage system description The drainage design inputs used in simulations are shown in Table 3.1. The drainage coefficient and depth to impermeable layer were 2.8 cm/day and 203 cm, respectively. The drain pipes were a 100-mm diameter four-row perforated pipe, with an estimated effective radius of 0.7 cm based on Ghane (2022). The maximum capacity of the subirrigation pump was considered to be 2.8 cm/day which is the same as the drainage coefficient. Table 3.1. Drainage system inputs used in DRAINMOD. Drainage system parameters 82 Drain depth (cm) 1005 Drain Spacing (cm) 2.8 Drainage coefficient (cm/day) Depth to impermeable layer (cm) 203 Effective radius for 4-inch perforated pipe (cm) 0.7 2.8 Subirrigation pump capacity (cm/day) Surface Storage Maximum surface storage (cm) Kirkham depth (cm) 2 1 72 3.3.3.2. DRAINMOD soil input data The dominant soil type in the field was Ziegenfuss clay loam, which is classified as a poorly drained soil (Natural Resources Conservation Service, 2023). Table 3.2 presents soil physical properties derived from both the SSURGO database and previous onsite measurements. To execute DRAINMOD, soil-water characteristics, Green-Ampt infiltration model parameters, the drainage volume–water-table depth and the water upflux–water-table depth relationships are also required. The soil utility package integrated into DRAINMOD, which uses the pedotransfer estimation software, ROSETTA (Schaap et al., 2001), facilitated the estimation of these hydraulic properties. ROSETTA determines the pedotransfer function parameters based on the soil series' texture and water holding capacity. %Sanda %Silta %Claya Table 3.2. Soil physical properties used in the simulation. Top and bottom depth of soil layer (cm) 0 – 23 23 – 76 76 – 102 102 - 170 170 - 203 a: These values were obtained from gSURRGO. b: These values were measured. c: These values are the calibrated values used in DRAINMOD. 37.0 35.5 33.5 32.0 29.0 35.8 37.3 33.5 33.0 27.0 27.2 27.2 25.7 35.0 34.0 Soil texture Bulk density Clay loam Clay loam Clay loam Clay loam Clay loam (g/cm3)b 1.37 1.46 1.50 1.72 1.90 Ks (cm/h)c 2.0 1.0 0.5 0.3 0.2 3.3.3.3. Historical and future climate data used in DRAINMOD To provide DRAINMOD with precipitation data, we utilized the disaggregation utility within DRAINMOD, which evenly disperses daily rainfall across a duration defined by the user (Ale et al., 2009; Singh et al., 2006; Skaggs et al., 2012). We assumed that daily precipitation took place from 01:00 am through 5:00 am, constituting a 4-hour period. Consequently, we disaggregated the quantity of daily precipitation within this specified time window. Daily precipitation and temperature data were measured by ATMOS-41 (METER Group) weather unit in the field. 73 We used 30 years of historical weather parameters (daily precipitation, daily maximum air temperature, and daily minimum air temperature) from 1994 to 2023, sourced from the National Oceanic and Atmospheric Administration (NOAA) and the nearest Enviroweather station. Enviroweather is a weather station network developed and supported by Michigan State University (MSU Enviroweather, 2023). These data were then used to conduct hydrology and crop yield simulations for selecting the optimum drainage design for the past. Subsequently, we utilized 30 years of projected weather data for the future spanning from 2030 to 2059 to predict hydrology and crop yields for determining the optimum drainage design for the future. For the future climate data, the NASA Earth Exchange Global Daily Downscaled Projections (NEX-GDDP) dataset sourced from the NASA Center for Climate Simulation (NCCS) was used in this study. This dataset is comprised of downscaled climate scenarios derived from 35 General Circulation Models (GCMs) (NEX-GDDP-CMIP6, 2023). From these GCMs, we selected 27 based on data availability for the chosen greenhouse gas emission scenario (Table A3.1). The analysis of GCM outputs was conducted within the framework of the Coupled Model Intercomparison Project phase 6 (CMIP6), which is supported by the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (Calvin et al., 2023). All 27 selected GCMs project daily precipitation, daily maximum air temperature, and daily minimum air temperatures at a spatial resolution of 0.25 degrees × 025 degrees (approximately 25 km × 25 km). To downscale the future climate data projections, the NEX-GDDP dataset employs the bias-correction spatial disaggregation which is a statistical downscaling algorithm. The GCM runs were executed for the shared socioeconomic pathway (SSP) greenhouse gas emission scenario, specifically SSP2-4.5, representing a medium pathway for the near future (2030–2059). Under SSP2-4.5, greenhouse gas emissions stabilize at 4.5 W/m2, with global emission reduction policies reaching their peak by 74 2050, and the world population stabilizing at 9 billion (Thomson et al., 2011). After acquiring and processing the future climate data, simulations were conducted using DRAINMOD for each GCM, and we used the equal-weighting multi-model ensemble (MME) approach to analyze the resultant outputs (Pease et al., 2017; Robertson et al., 2004; Shokrana et al., 2023; Tebaldi & Knutti, 2007; Weigel et al., 2010). MME averages the outputs from various GCMs. In this study, we followed meteorological seasons, which defines winter from December to February, spring from March to May, summer from June to August, and fall from September to November. Also, we considered the period from November through May as non-growing season and from June through October as growing season. 3.3.3.4. Crop and trafficability parameters in DRAINMOD Simulation of relative crop yield in DRAINMOD necessitates several parameters, including the root zone’s soil water content at wilting point, the susceptibility factor for drought period, the desired planting date, duration of growing season, the threshold WTD, and the relationship between effective root depth and days after planting (Skaggs, 1980). For the soil water content at wilting point, we initially used the gSSURGO database and subsequently calibrated (Table 3.3). During the calibration and validation periods (October 2018 to December 2023), we used the crop files based on the planting system including cover crops, however, for long-term historical and future simulations, we only considered corn as it is the major crop cultivated in southeast Michigan. Since the average corn planting date in our site was May 15th, we selected this date as the desired planting date for all simulation years. In DRAINMOD, the length of growing season was considered to be 170 days, reflecting the average duration observed in our experimental fields, was selected. Consequently, simulations were performed over a cropping window spanning from April to October, aligning with the selected planting window and growing season length. 75 In DRAINMOD, the relationship between effective root depth and days after planting plays a crucial role in determining the depth at which water is drawn to meet evapotranspiration (ETp) requirements. As the growing season progresses, the effective corn rooting depth gradually increases and reaches its maximum effective root depth of 40 cm after 80 days of planting. Following recommendations by Skaggs (1980) for corn, we considered an effective root depth of 3 cm for a fallow period and a crop yield limiting WTD of 30 cm. Field trafficability parameters exert a significant influence on field operations, including planting activities. DRAINMOD adjusts the planting day based on soil moisture conditions that align with specified trafficability parameters. These parameters encompass factors such as the minimum required air volume for field operations, the threshold rainfall amount to postpone work, and the duration of delay following precipitation events. In this study, the values for these trafficability parameters were selected based on the goodness of fit during the calibration. For the minimum soil air volume required to facilitate field operations, a value of 3 cm was utilized for both planting and harvest periods. Additionally, planting and harvest activities were postponed by 1 day following daily precipitation amounts exceeding 1 cm. 76 30 0.19 40 Table 3.3. Crop and trafficability inputs used in simulations. Crop Limiting water-table depth (cm) Root zone lower limit water table (cm3/cm3) Maximum effective root depth (cm) Relative yield simulation Cropping window Desired planting date Growing season length Yield intercept for crop wet stress Slope for yield vs. crop wet stress Drought stress yield intercept Drought stress yield slope Trafficability Minimum required air volume (cm) First period Second period Minimum rain to delay work (cm) First period Second period Delay after rain to restart work (d) April - October May 15th 170 102 0.75 100 1.22 1 1 1 3 3 3.3.3.5. DRAINMOD’s calibration and validation We calibrated and validated DRAINMOD to predict drainage discharge, WTD, observed from October 2018 to December 2023 at the experimental site. The calibration and validation were done based on goodness of fit for daily drainage discharge, daily WTD, and relative yield. The calibration period for drainage discharge was from October 2018 to December 2022. For WTD, the calibration period started from April 14th, 2021, due to missing data from October 2018 to April 2021. The validation period was from January 2023 to December 2023 for both drainage discharge and WTD. Throughout both the calibration and validation periods, the drainage system operated under free drainage conditions. The primary calibration parameters included the lateral saturated hydraulic conductivity, upward flux, drained volume, maximum surface storage, and Kirkham depth. By fine-tuning these parameters, adjustments were made to refine the relationship between WTD and drained water volume, as well as the upflux and parameters of the Green-Ampt model. These adjustments aimed to enhance the alignment between measured and predicted 77 drainage discharge and WTD. 3.3.4. Evaluating the effect of climate change on drainage design To examine the effect of climate change on the drainage design, we investigated various drainage designs including two drain depths of 75 cm and 125 cm for the past and future with continuous corn cropping system. For each drain depth, the drain spacings varied from 5 m to 50 m with 1-m increment using both historical and future climate data. The drain depths of 75 cm and 125 cm were selected as they represented shallow and deep drains. Ghane et al. (2021), Ghane and Askar (2021), and Skaggs and Nassehzadeh-Tabrizi (1986) expanded the drain spacing range from 5 m to 100 m. However, given the heavy-textured soil in this study, we did not exceed 50 m drain spacing. 3.3.5. Economic analysis of various drainage designs We performed an economic analysis to determine the optimum drain spacing for maximizing the average annual return on investment (AAR) for each drain depth, following the methodology outlined by Ghane et al. (2021). For every depth considered, and across spacings from 5 m to 50 m, AAR in dollars per hectare ($/ha) was computed as follows: 𝐴𝐴𝑅 = 𝐴𝐴𝐶𝑃𝐼 − 𝐴𝐶𝑃𝐶 − 𝐴𝐴𝐷𝑆𝐶 (2) where AACPI is the average annual corn production income ($/ha), ACPC is the annual corn production cost ($/ha), and AADSC is the average annual drainage system cost ($/ha). We assumed the annual corn production cost (ACPC) to be $2,225/ha ($901/ac) according to the 2023 cost (Table 3.4) (LaPorte, 2024). We obtained the average annual corn production income (AACPI) ($/ha) as follows 𝐴𝐴𝐶𝑃𝐼 = 𝑅𝑌 × 𝑃𝑌 × 𝐶𝑃 (3) where RY is the relative yield (decimal), PY is the potential yield (kg/ha), and CP is corn price 78 ($/kg). DRAINMOD applies four crop stresses when estimating the relative corn yield for each year. The corn price was assumed to be $0.25/kg ($6.01/bu) based on the 2023 average corn price (USDA National Agricultural Statistics Service, 2024) with a potential corn yield of 14,000 kg/ha according to LaPorte (2020). The potential yield signifies the average yield achievable without any stresses. The average annual drainage system cost (AADSC) ($/ha) was calculated (Huffman et al., 2013) as follows 𝐴𝐴𝐷𝑆𝐶 = 10000 𝑆0 × 𝐼𝐶𝐷𝑆( 𝐼𝑅𝐴𝐼 2 + 𝐷𝐸𝑃 + 𝑀𝑅) (4) where S0 is the optimum drain spacing (m), ICDS is the initial cost of the drainage system that is $2.33/m (Table 3.4), RIAI is the 5% interest rate on average investment per year (decimal), DEP is the 3.3% depreciation per year based on a 30-year expected lifetime of the system (decimal), and MR is the 0.3% maintenance rate per year (decimal). Table 3.4. Corn price, corn production expenses (LaPorte, 2024), and initial cost of the drainage system used in the economic analysis. Items Income Corn price (CP) Corn production expenses Variable (direct) expenses (seed, fertilizer, crop chemicals, crop insurance, repairs, fuel, labor, etc.) Fixed (indirect) expenses (land rental, depreciation, etc.) Total annual corn production cost (ACPC) (sum of variable and fixed expenses) Initial cost of the drainage system 100-mm-diameter drain pipe materiala Drain pipe installation costb Initial cost of drainage system (ICDS) (sum of material and installation) a The cost of drain pipe material was sourced from a Michigan pipe manufacturer in 2023. b The cost of drain pipe installation was sourced from a Michigan drainage contractor in 2023. $1,446/ha $779/ha $2,225/ha $1.18/m $1.15/m $2.33/m Unit Price $0.25/kg 79 3.3.6. Selecting the optimum drain spacing The drainage design, characterized by drain depth and spacing, yielding the highest AAR (an example depicted in Figure 3.2), was identified as the optimum drainage configuration. Notably, during this process, we observed that altering the annual corn production cost (ACPC) led to vertical shifts in the AAR graph presented in Figure 3.2, while maintaining other economic and input parameters constant. Consequently, variations in ACPC resulted in the same optimum drain spacing, indicating that changes in ACPC did not affect the horizontal positioning of the graph on the x-axis with drain spacing. We also investigated a scenario wherein farmers use the optimum drainage design derived from historical data for future. In essence, this scenario reflects a situation where farmers maintain the historical/current drainage design despite changes in future climate parameters. Therefore, the optimum drainage designs were selected for the past (1994 to 2023) and future (2030 to 2059), and the performance of the historical optimum design in the future (2030 to 2059) representing the scenario in which farmers do not change their design. 80 Figure 3.2. The impact of drain spacing on the 30-year average relative corn yield and average annual return on investment, as estimated by DRAINMOD, for both shallow and deep drains depths. In this example, the average annual return graph exhibits peaks at drain spacings of 500 cm and 1000 cm for drain depths of 75 cm and 125 cm, respectively. These peaks signify the optimum drain spacings corresponding to the maximum average annual return on investment. 3.3.7. Suitability of subirrigation We examined the efficacy of subirrigation as a supplemental irrigation to mitigate yield loss due to drought stress. The subirrigation applied to the scenario, in which farmers do not change the drainage design and design the drainage system based on historical weather data. In addition to the effect of subirrigation on the crop yield, we also investigated how the subirrigation management (start time and weir depth) could be different between past and future. For this purpose, three start times of June 1st (early June), June 15th (mid-June), and June 29th (late June) and four subirrigation weir depths of 70 cm, 65 cm, 60 cm, and 55 cm were investigated. Regardless of the starting time and weir depth, subirrigation was applied during July and August 81 (when ET is the highest) for all scenarios. The maximum subirrigation pump capacity was considered to be equal to the drainage coefficient, which was 2.8 cm/day. 3.3.8. Statistical analysis Some of the datasets did not meet the normality assumption according to the Shapiro-Wilk test (Shapiro & Wilk, 1965). Therefore, the non-parametric two-sample Wilcoxon-rank test (Wilcoxon, 1992) was used to compare and identify the significance of differences between the weather parameters, hydrology and crop yield in past and future. The null hypotheses were no significant change in weather parameters, hydrology and crop yield in past and future. The tests were conducted at a 5% level of significance (α = 0.05). To evaluate the performance of DRAINMOD, we used Nash-Sutcliffe Efficiency (NSE), Kling-Gupta Efficiency (KGE), Percent Bias (PBIAS), and Mean Absolute Error (MAE) as follow 𝑁𝑆𝐸 = 1 − 𝑛 ∑ (𝑂𝑖−𝑃𝑖)2 𝑖=1 𝑛 ∑ (𝑂𝑖−𝑂̅)2 𝑖=1 (5) 𝐾𝐺𝐸 = 1 − √(𝑟 − 1)2 + (𝛼 − 1)2 + (𝛽 − 1)2 (6) 𝑃𝐵𝐼𝐴𝑆 = 𝑛 ∑ (𝑂𝑖−𝑃𝑖) 𝑖=1 𝑛 ∑ (𝑂𝑖) 𝑖=1 × 100 (7) 𝑀𝐴𝐸 = 𝑛 𝑖=1 ∑ |𝑂𝑖−𝑃𝑖| 𝑛 (8) where 𝑂𝑖 is the observed value, 𝑃𝑖 is the predicted value, 𝑟 is Pearson correlation coefficient, 𝛼 is a term representing the variability of prediction errors, and 𝛽 is a bias term. Given that the NSE metric tends to be overly responsive to extreme values while less responsive to lower values (Moriasi et al., 2015), we also incorporated the KGE in our evaluation. Evaluation criteria for assessing the level of agreement between DRAINMOD predictions of drainage discharge and WTD and their corresponding measurements are outlined by Skaggs et al. (2012) and the KGE metric. 82 We also used the Coefficient of Variation (CV) to evaluate the efficiency of subirrigation in stabilizing the crop yield in the future as follow 𝐶𝑉 = 𝜎 𝜇 (9) where 𝜎 is the population standard deviation and 𝜇 is the population mean. 3.4. Results and Discussions 3.4.1. Comparing temperature and precipitation the past and future 3.4.1.1. Air temperature of the past and future The future (2030–2059) and historical (1994–2023) maximum, minimum, and average temperatures were analyzed and compared (Table 3.5). All 27 GCMs indicated statistically significant increase in the monthly maximum, minimum, and average air temperatures in the future in comparison to the past (p-values < 0.001) (Figure 3.3 and Figure A3.1). The 30-year average annual temperature was also significantly increased and different (p-values < 0.001). On a seasonal scale, the most substantial average temperature rise occurred during fall, with a 2.8°C increase, while the smallest rise was observed in winter, with a 2.1°C increase. Additionally, the GCMs forecasted a 2.6°C increase in average seasonal temperatures for spring and a 2.5°C increase for summer. During the historical period, the 30-year average temperatures for non-growing season was 3.3°C and for growing season was 18.6°C. Based on the projections, however, the future non- growing season and growing season temperatures were 5.7°C (2.4°C increase) and 21.2°C (2.6°C increase), respectively. 83 Table 3.5. Comparing the 30-year average monthly temperature and precipitation between historical (1994-2023) and future (2030–2059) climate scenarios. Increases are denoted by positive percentages, while decreases are represented by negative percentages. Historical weather (1994-2023) Future weather (2030-2059) Month Maximum temperature (℃) Minimum temperature (℃) Precipitation (mm) Maximum temperature (℃) Minimum temperature (℃) January February March April May June July August 0.3 2.0 8.1 15.2 21.6 26.9 28.8 27.6 September 24.2 October 16.9 November December 9.4 3.1 Annual 15.3 -8.2 -7.7 -3.0 2.5 8.4 13.9 15.9 15.2 11.0 5.0 -0.7 -4.8 4 53.9 54.5 63.3 87.5 101.0 107.0 88.6 90.2 83.1 76.4 66.1 61.5 933.2 2.6 4.7 10.7 18.0 24.4 29.5 31.7 30.8 27.2 20.3 12.2 5.0 18.1 -6.2 -5.1 -0.6 5.0 10.9 15.9 18.2 17.5 13.7 7.5 1.9 -3.5 6.3 Change in average temperature (℃) Precipitation (mm) 58.3 50.5 75.9 87.0 94.1 90.8 93.3 90.7 75.5 65.2 72.4 74.0 927.7 2.1 2.7 2.5 2.7 2.6 2.3 2.5 2.8 2.9 3.0 2.7 1.6 2.5 Change in precipitation (%) 8.2 -7.4 19.9 -0.5 -6.9 -15.2 5.3 0.5 -9.2 -14.6 9.5 20.5 -0.6 Figure 3.3. Heatmap for maximum and minimum temperatures increase in the future. 84 3.4.1.2. Precipitation of the past and future The future and historical precipitations were compared, revealing that the future 30-year average annual precipitation projected by 27 GCMs is anticipated to decrease insignificantly (p- value = 0.83) from 933 mm that was in past to 928 mm in the future (Table 3.5). Notably, this slight decline in precipitation contrasts with Pease et al. (2017) findings showing a 4.4% increase in precipitation based on their moderate RCP 4.5 scenario. Discrepancies in outcomes could arise due to variations in modeling frameworks (CMIP5 vs. CMIP6) and historical baseline periods (1971 - 2000 vs. 1994 - 2023) utilized in these studies. To confirm our reasoning for the differences, we calculated the long-term average annual precipitation from 1971 to 2000 for Blissfield, which came to 888 mm. This shows that precipitation increased from the older baseline (1971-2000) to the more recent baseline (1994-2023). However, the future prediction shows that precipitation will remain steady compared to the more recent baseline. Therefore, when comparing changes in precipitation with climate change, it is important to compare the future precipitation with the more recent baseline period. Overall, results show that the average long-term precipitation has increased from the period of 1971-2000 to 1994-2023, after which precipitation forecast shows no significant change for 2030-2059. While the 30-year average monthly precipitation amount during winter and spring increased, the precipitation amounts during the summer and fall seasons decreased (Figure 3.4). Projections indicate a 2.1% increase in precipitation during spring and a notable 12.8% increase during winter. Conversely, the 30-year average precipitation is expected to decrease by 3.9% during summer and 7.3% during fall. The percent change in precipitation amount is reported to be the highest in December and the lowest in August (Table 3.5). In the future, non-growing season precipitation is estimated to be 1.4% higher (60.3 cm vs. 59.5 cm) than in the historical period, 85 whereas growing season precipitation is forecasted to decrease by 6.7% (41.6 cm vs. 44.6 cm). Figure 3.4. Changes in average monthly precipitation amounts under historical (1994–2023) and future (2030–2059) periods. Regarding the variation in precipitation projections, it is notable that the distances between the upper and lower whiskers are small for January, February, June, and November. This shows that that the precipitation projections from the GCMs are relatively similar for these months, suggesting more confidence in the predictions of these months (Figure 3.5). However, outliers are apparent in February and November. Conversely, the largest distance between the upper and lower whiskers was observed in May, indicating a substantial difference between the minimum and maximum values projected. 86 Figure 3.5. Box plot illustrating the variation in precipitation predictions from 27 GCMs for the future period (2030–2059). The upper and lower whiskers delineate the range of data, while the upper and lower quartiles represent the 25th to 75th percentile range. 3.4.2. Results of calibration and validation of DRAINMOD DRAINMOD was calibrated and validated based on the goodness of fit between the observed and predicted drainage discharge and WTD (Table 3.6 and Figure 3.6). In the calibration period, NSE and KGE for drainage discharge were 0.71 and 0.73, respectively, indicating a good agreement between the estimated and observed data. The PBIAS was -0.85%, reflecting an excellent (< 5%) level of estimation accuracy. During the validation period, the NSE and KGE for drainage discharge were 0.63 and 0.58, respectively, signifying good and acceptable agreements. The PBIAS value of 0.97 further showed the model's excellent proficiency in estimation. 87 Table 3.6. DRAINMOD simulation statistics based on daily data. Estimated Parameter Period PBIAS Criteria MAE Criteria NSE Criteria KGE Criteria Drainage discharge Water- table depth Calibration -0.85 Excellent 0.05 Excellent 0.71 Good 0.73 Good Validation 0.97 Excellent 0.05 Excellent 0.63 Good 0.58 Acceptable Calibration -12.48 Good 14.79 Acceptable 0.60 Acceptable 0.71 Good Validation -1.18 Excellent 17.44 Acceptable 0.64 Good 0.55 Acceptable Figure 3.6. The observed and estimated daily drainage discharge during the calibration and validation periods. Regarding WTD in the calibration phase, the NSE and KGE were 0.60 and 0.70, respectively, indicating good agreements between the estimated and measured data. The PBIAS was -12.48, indicative of good estimation accuracy. In the validation period for drainage discharge, the NSE and KGE were 0.64 and 0.55, denoting good and acceptable agreements, respectively. The PBIAS of -1.18, being less than 5%, characterized the estimation as excellent (Table 3.6 and Figure 3.7). 88 Figure 3.7. The observed and estimated daily water-table depth during the calibration and validation periods. 3.4.3. Optimum drainage design for past and future 3.4.3.1. Optimum drainage design based on the past and future weathers The drain spacing of 500 cm was selected as the optimum drain spacing for the drain depth of 75 cm (shallow drains) and 1000 cm for the drain depth of 125 cm (deep drains) when using historical weather data (Table 3.7). However, the highest AARs were estimated to be achieved by wider drain spacings for the future. The drain spacing of 800 cm was selected as the optimum drain spacing for the drain depth of 75 cm and 1600 cm spacing for drain depth of 125 cm when using future weather data. Results show that in the future, a wider drain spacing would reduce the drainage intensity, retaining the water in the soil profile for a longer period, and consequently alleviating the crop yield loss due to drought stress. For both past and future, the highest AAR among optimum drainage designs was observed for deep drains. The reasons for this are, first, the lower initial cost of the drainage system due to wider drain spacings, and second, the higher corn yield due to lower yield loss due to drought stress. Although the drainage installation costs were lower due to wider drain spacings in the future, the AARs were lower than the AARs for the historical optimum drain spacings because of lower yield (85% in the past vs 80% in the future). Lower yield in the future could be attributed to the higher temperature and a different precipitation pattern. In general, if a field does not already have a drainage system installed, it may be more 89 advantageous to implement a drainage system designed based on future weather data rather than past data. Table 3.7. Optimum drain spacings based on historical weather (1994-2023), future weather (2030-2059), and future with the historical optimum drain spacings. Scenario Time interval Drain depth = 75 cm shallow drains Spacing (cm) Yield (%) AAR ($/ha) Drain depth = 125 cm deep drains Yield (%) Spacing (cm) AAR ($/ha) Drainage design based on historical weather data Drainage design based on future weather data 1994 – 2023 500 2030 – 2059 800 82 78 328 267 1000 85 581 1600 80 442 3.4.3.2. Past optimum drainage design in the future Another scenario pertains to fields with an existing drainage system designed based on past weather data, where it is nearly impossible to modify the drainage design in such fields. In the case of not changing the drainage design and using the drainage system designed based on the historical weather data, a considerable reduction in AAR from $581/ha to $353/ha was observed for deep drains which was due to lower relative crop yield (85% vs 79%). The AAR for shallow drains also decreased from $328/ha to $206/ha which was due to the decrease in relative yield (82% vs 79%). From this section onward, the analyses and discussions are based on maintaining the existing drainage system designed based on the historical weather data. In other words, we predicted the future hydrology and crop yield under the condition where the drainage system is designed based on the historical data to evaluate its performance in the future. 90 3.4.4. Effect of climate change on hydrology and crop yield 3.4.4.1. Water-table depth for past and future The effect of climate change on WTD for optimum drainage design based on the historical climate was investigated. The WTD was predicted to be significantly deeper in the future for both shallow and deep drains (8% on average; p-value < 0.001), which was mainly due to increased temperature and increased ET (Figure 3.8 and Figure A3.2). In clay loam soils, as the evaporation from the soil surface and crop transpiration increases, the water depletion increases due to increased upward flux via capillary rise. The deep drains had a deeper WTD than the shallow drains because of its deeper drain depth. For shallow drains, the highest difference in the historical and future WTD occurred during fall, while deep drains showed the highest difference in summer. During the growing season, the future WTD for shallow and deep drains were 8.63% and 7.64% deeper than the historical, respectively. Our findings agreed with the findings of Sojka et al. (2020) as they reported that the effect of climate change will be a decrease in the mean groundwater table in the fields equipped with drainage system. Salem et al. (2018) showed that higher increase in temperature and no appreciable change in rainfall in winter will cause further declination of groundwater levels in an agricultural area in Bangladesh. 91 Figure 3.8. The 30-year average monthly water-table depth for historical (1994–2023) and future (2030–2059) periods under free drainage for shallow drains. 3.4.4.2. Actual evapotranspiration for the past and future DRAINMOD simulations showed that climate change affected the 30-year average annual actual ET (Figure 3.9 and Figure A3.3). Our findings showed that the future average annual ET for both drainage designs will significantly increase (shallow drains p-value = 0.0068; deep drains p-value = 0.0029). The 30-year historical average annual ET was 54 cm, while, 30-year future average annual ET was 63 cm (17% increase). On average, ET during the growing season and non- growing season increased 8.1% and 28.8%, respectively. The shallow drains had higher historical and future ETs as well as higher change than deep drains which was mainly due to its shallower water table. The main reason for the increased ET in the future is the increased temperature. As explained in section 3.4.1.1, the future temperature will significantly increase which leads to an 92 increase in ET. Pease et al. (2017) reported increased temperature and ET in the future for the Midwest. Jiang et al. (2020) stated that without the influence of elevated CO2, higher future temperatures would result in increased ET in the future. Van Roosmalen et al. (2009) reported substantially higher ET in the future. Wang et al. (2015) projected higher actual and potential ET in the future in Iowa (from 2045 to 2064). However, some studies reported decreased actual ET for the future. Shokrana et al. (2023) reported a slight decrease in future ET as compared with the historical ET. They attributed the decreased ET to the increased RH and subsequently increased vapor pressure in the future. They utilized RZWQM which uses the Shuttleworth–Wallace equation to estimate the daily potential ET. The RH is one of the inputs of the Shuttleworth– Wallace equation, however, in the Thornthwaite equation that DRAINMOD uses to estimate potential ET, temperature is the main input parameter. Figure 3.9. The 30-year average monthly actual evapotranspiration for historical (1994–2023) and future (2030–2059) periods under free drainage for shallow drains. 93 The average annual number of dry days during growing season for future for both drainage designs was considerably higher than in the past (Table A3.2). A dry day is computed in DRAINMOD whenever actual ET is less than potential ET because of limiting soil water conditions. In the past, dry days during growing season for shallow and deep drains were 24 days and 30 days, respectively. However, in the future, we observed 38 dry days for shallow drains and 43 dry days for deep drains. The reason for having fewer dry days in shallow drains is that, for shallow drains, the water table was shallow, providing the crop with enough water to satisfy the water demand for ET (Ghane & Askar, 2021). However, in the case of deep drains, the crop had limited access to water, causing more days when the crop could not achieve its potential ET rate. The ratio of ET to precipitation will increase in the future. In the past, on average, 59% and 57% of precipitation were lost via ET for shallow and deep drains, respectively, while in the future, the ratios are 70% for shallow drains and 67% for deep drains. 3.4.4.3. Drainage discharge for the past and future The effect of climate change on the drainage discharge for both drainage designs was investigated. For both drain depths, the drainage discharge of the future was significantly (shallow drains p-value = 0.0068; deep drains p-value = 0.0067) lower than the historical period (Figure 3.10 and Figure A3.4). The highest reduction in drainage discharge was observed in fall for shallow drain (72% on average) and in summer for deep drain (77% on average). The growing season also had a higher decrease in drainage discharge than non-growing season, on average 72% for both drainage designs. The decreased drainage discharge is attributed to the decreased precipitation (explained in section 3.4.1.2) and the increased ET (explained in section 3.4.4.2). Pease et al. (2017) and Jiang et al. (2020) also reported a decrease in drainage discharge for the future. They attributed this to the higher ET caused by increased temperatures under climate change, resulting 94 in less water being lost to drainage. However, some studies, such as Shokrana et al. (2023), Wang et al. (2015), and Singh et al. (2009) reported increased drainage discharge in the future. Shokrana et al. (2023) showed more precipitation will turn into the drainage discharge as ET will decrease which is in contrast with our findings provided in section 3.4.4.2. Singh et al. (2009) investigated the effect of climate change on drainage discharge under two future climate scenarios. Given the difference in the nature of projections of each GCM, a wide variety of responses from decreased discharge to increased discharge is expected. In this study, we investigated 27 GCMs and the final results provided are the average of all responses. For both past and future, the highest and lowest drainage discharge was estimated in March and September, respectively, for shallow drain. For deep drain also the highest drainage discharge was predicted to be in March for both past and future. However, the lowest drainage discharge was expected to happen in September in the past and July in the future. The highest variability of estimated drainage discharge based on different GCMs was estimated to be in July from 2.48 cm to 9.62 cm for shallow drains. The ratio of drainage discharge to precipitation and surface runoff to precipitation decreased in the future. In the past, on average, 38% and 42% of precipitation turned into drainage discharge for shallow and deep drains, respectively, while in the future, the ratios are 29% for shallow drains and 33% for deep drains. The runoff was not considerable in both the past and future due to the field’s low slope and high surface storage. On average, only 1.4% of precipitation transformed into runoff in the past, while this proportion is expected to decrease to 0.6% in the future. 95 Figure 3.10. The 30-year average monthly drainage discharge for historical (1994–2023) and future (2030–2059) periods under free drainage for shallow drains. 3.4.4.4. Crop yield for the past and future The historical 30-year average crop yield for shallow drains was 82% and for deep drains was 85%. For the future, the 30-year average yield for shallow and deep drains were 78% and 80%, respectively. However, the difference in yield from past to the future was not significant in both drainage designs (shallow drains p-value = 0.19; deep drains p-value = 0.11). Figure 3.11 depicts the relative yield with excess water stress and drought stress in the past and future (based on four random climate models) under free drainage for shallow drain. As explained in section 3.4.3.1 and can be seen in Figure 3.11, the main reason for the yield loss in the past was excess water, while drought stress was the major factor decreasing the yield in the future. In the past, for 96 shallow drains, the yield losses due to excess water and drought stresses were 2.7% and 13%, respectively. In the future, the yield loss due to excess water is projected to be negligible, while the yield loss due to drought stress is expected to reach 20%. Similarly, in the past for deep drains, the yield losses due to excess water and drought stresses were 1.5% and 13%, respectively. The yield loss due to drought stress is anticipated to be 20%, and excess water is not expected to have an effect on crop yield. Figure 3.3. Relative yield for continuous corn with excess water stress and drought stress in the past and future under free drainage for shallow drain. As an example, we present the future GCM for NorESM2-MM, ACCESS-CM2, CanESM5, and EC-Earth3 climate models. The annual average of working days for the future for both drainage designs was higher than in the past. In DRAINMOD, a day is counted as a working day if the drained or water free pore space is greater than a threshold value, if there is less than a given amount of rain that day, and if rainfall sufficient to delay field work has not occurred within a given number of days. In the past, the working days for shallow and deep drains were 85 days and 99 days, respectively. 97 However, in the future, we observed 92 working days for shallow drains and 106 working days for deep drains. The deep drains had more working days for both past and future. The reason is that for deep drains, the water table was the deepest providing a better condition for machinery traffic, in comparison to shallow drains which had the shallowest water table. One way to alleviate the drought stress is to provide supplemental irrigation such as subirrigation which is discussed in section 3.4.5. 3.4.5. Efficiency of subirrigation practice in the past and future Subirrigation was applied to the past and future under different management scenarios to investigate how subirrigation management would change from the past to the future (Table 3.8). For both past and future, subirrigation was able to significantly (p-value < 0.001 for both drainage designs) increase the crop yield, however, this increase was more pronounced in the future than the past (Figure 3.12). Regardless of the weir depth and irrigation start time, in the past, subirrigation increased the relative crop yield by 11% under shallow drains and by 7% under deep drains. However, in the future, subirrigation is expected to increase the relative crop yield on average by 19% (ranging from 8% to 36%) and 17% (ranging from 3% to 31%) under shallow drains and deep drains, respectively (Table 3.8). Table 3.8. Relative yield, percent increase, optimum subirrigation management in the past and future. Period Relative yield under free drainage (%) Relative yield under subirrigation (%) Percent increase Optimum irrigation start time Optimum weir depth (cm) Drain depth = 75 cm, Drain spacing = 500 cm Historical 82 93 11 Future 78 (61 - 90) 97 (93 - 99) 19 (8 - 36) 29-Jun 29-Jun Drain depth = 125 cm, Drain spacing = 1000 cm Historical 85 92 7 Future 80 (61 – 94) 97 (89 - 99) 17 (3 - 31) 29-Jun 15-Jun Values inside the parenthesis are lower and upper limits. 70 60 70 60 98 The higher efficiency of subirrigation in the future compared to the past can be attributed to the fact that, historically, both excess water and drought stresses contributed to yield loss. While subirrigation effectively mitigated drought stress during dry years in the past, it introduced unwanted wet stress to the crop in years with sufficient precipitation. Consequently, the cumulative impact of subirrigation resulted in a slight increase in relative crop yield compared to free drainage in the past. However, in the future, increased temperatures and shifts in precipitation patterns indicate that drought stress will be the predominant factor influencing crop yield loss. Wet stress is expected to play a negligible role in this future scenario. As a result, the application of subirrigation in the future has the potential to significantly enhance relative crop yield. Gunn et al. (2018) and Allred et al. (2014) found a similar response to subirrigation practice in Ohio and confirmed its efficacy in increasing the crop yield. Murugaboopathi et al. (1995) also found that corn subirrigation practice could increase the crop yield in North Carolina soils. 99 Figure 3.12. Relative yield (%) for past and future under free drainage and subirrigation. 100 Our results showed that subirrigation not only increased the yield for both drainage designs, but also it stabilized the yield (Figure 3.12). On average, in the future, the coefficient of variation for relative yield decreased from 0.27 (ranging from 0.14 to 0.44) under free drainage to 0.06 (ranging from 0.02 to 0.20) with subirrigation in shallow drain. For deep drains, subirrigation managed to reduce the coefficient of variation for the relative yield from 0.24 (ranging from 0.10 to 0.43) to 0.05 (ranging from 0.01 to 0.10), on average. This occurred because subirrigation can decrease crop yield loss due to drought stress, which is projected to be the main reason for crop yield loss in the future. Gunn et al. (2018) reported less variation when adopting subirrigation practice. Cooper et al. (1999) also reported subirrigation can stabilize the corn yield at high levels, resulting in higher long term average yields. The average annual amount of water pumped into the subirrigation system also was considerably different in the past and future. In the past, the average annual pumped water was 6.9 cm and 6.5 cm for shallow and deep drains. In the future, the average annual water pumped into the subirrigation system was 9.9 cm (ranging from 6 cm to 15.7 cm) for shallow drains and 10.5 cm (ranging from 6.7 cm to 15.6 cm) for deep drains. 3.4.5.1. Subirrigation Management As explained in sections 3.4.1.1 and 3.4.1.2, not only did projections show increased temperature, but also, they showed considerably lower precipitation in June in the future. As a result, there has been a significant shift in subirrigation management practices from the past to the future (Table 3.8). Based on historical data, for shallow drains, the highest crop yield was achieved with a weir depth of 70 cm and irrigation starting in late June. However, due to increased drought stress in the future, optimum results are anticipated with a shallower weir depth (60 cm). A similar trend is evident for the drainage design of deep drains. In the past, the optimum weir depth and 101 irrigation start time were determined to be 70 cm and late June, respectively. Contrastingly, in the future, the ideal conditions involve a weir depth of 60 cm and a mid-June start for irrigation. 3.4.5.2. Subirrigation limitation Subirrigation introduces the capability to replenish water into the drainage system, mitigating yield losses attributed to drought stress. However, it necessitates a separate water source that meets requirements for both quantity and quality. Drainage water recycling (DWR) is an emerging practice wherein subsurface and surface drainage water are captured and stored for subsequent use as supplementary irrigation during dry spells (Frankenberger et al., 2017). Despite the myriad benefits that DWR offers, including enhanced crop yield, improved water quality, better downstream flow management, and support for wildlife, its implementation requires considerable investment (Hay et al., 2021). Constructing reservoirs large enough to provide adequate storage incurs considerable expenses. However, the cost-benefit evaluation of drainage water recycling has not been investigated. It is imperative to develop research-based irrigation management strategies tailored for poorly drained soils, particularly for farmers well-versed in drainage but lacking experience in irrigation. Strategies should focus on pinpointing sites where costs are minimized, and economic benefits are maximized. 3.5. Conclusions The effect of climate change on the drainage design and the efficacy of subirrigation practice in alleviating yield loss due to drought stress in southeast Michigan were investigated. The optimum drainage designs were compared for the past and future. The key findings are summarized as follows: 102 • Projections indicated a temperature increase of 2.4 °C during the non-growing season (from November through May) and a 2.7 °C rise during the growing season (from June through October) compared to the historical period. • According to the projections, the average annual precipitation is expected to remain relatively stable in the future. Anticipated seasonal changes include a 1.4% increase in precipitation during the non-growing season compared to the historical period and a 6.7% decrease in precipitation during the growing season. Projections indicate a 2.1% increase in precipitation during spring and a notable 12.8% increase during winter. • In the past, the highest annual return on investments for drain depths of 75 cm and 125 cm occurred when the drain spacings were 500 cm and 1000 cm, respectively. However, for the future, the optimum drain spacing is projected to be wider at 800 cm for a drain depth of 75 cm and 1600 cm for a drain depth of 125 cm. • The 30-year average annual drainage discharge and water-table depth for both shallow and deep drains are anticipated to decrease significantly, attributed to the projected increase in evapotranspiration in the future. • Subirrigation is expected to reduce drought stress and increase the yield up to 31% in the future. • In the past, subirrigation initiation was feasible in late June with a weir depth of 70 cm. However, in the future, subirrigation is anticipated to be more advantageous when starting in mid-June, coupled with a shallower weir depth. In conclusion, wider drain spacings would be more efficient for the future. Due to shift in precipitation and increased temperature, a wider drain spacing can hold the water in the field by reducing the drainage intensity. In the case of maintaining the past or the current drainage design 103 for the future, the adoption of subirrigation may be a valuable practice to alleviate potential yield loss resulting from the drought stress. The value of this study lies in its contribution to understanding the future performance of existing drainage designs and the efficacy of subirrigation as a supplementary irrigation method, utilizing the latest CMIP6 modeling framework. 104 REFERENCES Adhikari, N., Davidson, P. C., Cooke, R. A., & Book, R. S. (2020). Drainmod-Linked Interface for Evaluating Drainage System Response to Climate Scenarios. Applied Engineering in Agriculture, 36(3), 303–319. https://doi.org/10.13031/aea.13383 Ale, S., Bowling, L. C., Brouder, S. M., Frankenberger, J. R., & Youssef, M. A. (2009). Simulated effect of drainage water management operational strategy on hydrology and crop yield for Drummer soil in the Midwestern United States. Agricultural Water Management, 96(4), 653–665. https://doi.org/10.1016/j.agwat.2008.10.005 Allred, B. J., Gamble, D. L., Clevenger, W. B., LaBarge, G. A., Prill, G. L., Czartoski, B. J., Fausey, N. R., & Brown, L. C. (2014). Crop Yield Summary for Three Wetland Reservoir Subirrigation Systems in Northwest Ohio. Applied Engineering in Agriculture, 889–903. https://doi.org/10.13031/aea.30.10501 Arnold, L. A. (2004). Effects of Drain Depth on Nitrogen Losses in Drainage in Shallow Water Table Soils. North Carolina State University. Calvin, K., Dasgupta, D., Krinner, G., Mukherji, A., Thorne, P. W., Trisos, C., Romero, J., Aldunce, P., Barrett, K., Blanco, G., Cheung, W. W. L., Connors, S., Denton, F., Diongue-Niang, A., Dodman, D., Garschagen, M., Geden, O., Hayward, B., Jones, C., … Ha, M. (2023). IPCC, 2023: Climate Change 2023: Synthesis Report. Contribution of Working Groups I, II and III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change [Core Writing Team, H. Lee and J. Romero (eds.)]. IPCC, Geneva, Switzerland. (P. Arias, M. Bustamante, I. Elgizouli, G. Flato, M. Howden, C. Méndez-Vallejo, J. J. Pereira, R. Pichs-Madruga, S. K. Rose, Y. Saheb, R. Sánchez Rodríguez, D. Ürge-Vorsatz, C. Xiao, N. Yassaa, J. Romero, J. Kim, E. F. Haites, Y. Jung, R. Stavins, … C. Péan, Eds.). https://doi.org/10.59327/IPCC/AR6-9789291691647 Cooper, R. L., Fausey, N. R., & Johnson, J. W. (1999). Yield Response of Corn to a Subirrigation/Drainage Management System in Northern Ohio. Journal of Production Agriculture, 12(1), 74–77. https://doi.org/10.2134/jpa1999.0074 Deryng, D., Conway, D., Ramankutty, N., Price, J., & Warren, R. (2014). Global crop yield response to extreme heat stress under multiple climate change futures. Environmental Research Letters, 9(3), 034011. https://doi.org/10.1088/1748-9326/9/3/034011 Drury, C. F., Tan, C. S., Reynolds, W. D., Welacky, T. W., Oloya, T. O., & Gaynor, J. D. (2009). Managing Tile Drainage, Subirrigation, and Nitrogen Fertilization to Enhance Crop Yields and Reduce Nitrate Loss. Journal of Environmental Quality, 38(3), 1193– 1204. https://doi.org/10.2134/jeq2008.0036 Ehsan Ghane. (2023). Drainage Design Tools. Michigan State University Extension. 105 Evans, R. O., Skaggs, R. W., & Sneed, R. E. (1991). STRESS DAY INDEX MODELS TO PREDICT CORN AND SOYBEAN RELATIVE YIELD UNDER HIGH WATER TABLE CONDITIONS. Transactions of the ASAE, 34(5), 1997–2005. https://doi.org/10.13031/2013.31829 Fisher, M. J., Fausey, N. R., Subler, S. E., Brown, L. C., & Bierman, P. M. (1999). Water Table Management, Nitrogen Dynamics, and Yields of Corn and Soybean. Soil Science Society of America Journal, 63(6), 1786–1795. https://doi.org/10.2136/sssaj1999.6361786x Frankenberger, J., Reinhart, B., Nelson, K., Bowling, L., Hay, C., Youssef, M., & Allred, B. (2017). Questions and answers about drainage water recycling for the Midwest.: Vol. ABE-156. Purdue University Extension. Galloway, J. N., Schlesinger, W. H., Clark, C. M., Grimm, N. B., Jackson, R. B., Law, B. E. ;, Thornton, P. E., Townsend, A. R., & Martin, R. (2014). Ch. 15: Biogeochemical Cycles. Climate Change Impacts in the United States: The Third National Climate Assessment. https://doi.org/10.7930/J0X63JT0 Geneva, S. (2013). Intergovernmental panel on climate change. Ghane, E., & Askar, M. H. (2021). Predicting the effect of drain depth on profitability and hydrology of subsurface drainage systems across the eastern USA. Agricultural Water Management, 258, 107072. https://doi.org/10.1016/j.agwat.2021.107072 Ghane, E., Askar, M. H., & Skaggs, R. W. (2021). Design drainage rates to optimize crop production for subsurface-drained fields. Agricultural Water Management, 257, 107045. https://doi.org/10.1016/j.agwat.2021.107045 Gunn, K. M., Baule, W. J., Frankenberger, J. R., Gamble, D. L., Allred, B. J., Andresen, J. A., & Brown, L. C. (2018). Modeled climate change impacts on subirrigated maize relative yield in northwest Ohio. Agricultural Water Management, 206, 56–66. https://doi.org/10.1016/j.agwat.2018.04.034 Hay, C. H., Reinhart, B. D., Frankenberger, J. R., Helmers, M. J., Jia, X., Nelson, K. A., & Youssef, M. A. (2021). Frontier: Drainage Water Recycling in the Humid Regions of the U.S.: Challenges and Opportunities. Transactions of the ASABE, 64(3), 1095–1102. https://doi.org/10.13031/trans.14207 Huffman, R. L., Fangmeier, D. D., Elliot, W. J., Workman, S. R., & Schwab, G. O. (2013). Soil and water conservation engineering. American Society of Agricultural and Biological Engineers. Jiang, Q., Qi, Z., Xue, L., Bukovsky, M., Madramootoo, C. A., & Smith, W. (2020). Assessing climate change impacts on greenhouse gas emissions, N losses in drainage and crop production in a subsurface drained field. Science of The Total Environment, 705, 135969. https://doi.org/10.1016/j.scitotenv.2019.135969 106 Kanwar, R. S., Melvin, S. W., & Sanoja, J. (1994). Simulating Corn Yields for Two Iowa Soils. Applied Engineering in Agriculture, 10(3), 373–379. https://doi.org/10.13031/2013.25866 LaPorte, J. P. (2020). Crop Budget Estimator Tool for Grains [WWW Document]. Malakshahi, A.-A., Darzi- Naftchali, A., & Mohseni, B. (2020). Analyzing water table depth fluctuation response to evapotranspiration involving DRAINMOD model. Agricultural Water Management, 234, 106125. https://doi.org/10.1016/j.agwat.2020.106125 Mejia, M. N., Madramootoo, C. A., & Broughton, R. S. (2000). Influence of water table management on corn and soybean yields. Agricultural Water Management, 46(1), 73–89. https://doi.org/10.1016/S0378-3774(99)00109-2 Moriasi, D. N., Gitau, M. W., Pai, N., & Daggupati, P. (2015). Hydrologic and Water Quality Models: Performance Measures and Evaluation Criteria. Transactions of the ASABE, 58(6), 1763–1785. https://doi.org/10.13031/trans.58.10715 Moustafa, A. T. A., Ahmed, W. H., Chang, A. C., & Hoffman, G. J. (1987). Influence of Depth and Spacing of Tile Drains on Crop Productivity in Nile Delta. American Society of Agricultural Engineers, 30(5), 1374–1377. MSU Enviroweather. (2023). Https://Legacy.Enviroweather.Msu.Edu/Weather.Php?Commodity=&stn=drf. Murugaboopathi, C., Skaggs, R. W., & Evans, R. O. (1995). Subirrigation for corn and soybean production for North Carolina soils (H. W. Belcher & F. M. D’Itri, Eds.). Lewis Publishers. Natural Resources Conservation Service, U. S. D. of Agriculture. (2023). Soil Survey Geographic (SSURGO) Database for Blissfield, Lenawee County, MI. Ng, H. Y. F., Tan, C. S., Drury, C. F., & Gaynor, J. D. (2002). Controlled drainage and subirrigation influences tile nitrate loss and corn yields in a sandy loam soil in Southwestern Ontario. Agriculture, Ecosystems & Environment, 90(1), 81–88. https://doi.org/10.1016/S0167-8809(01)00172-4 Pease, L. A., Fausey, N. R., Martin, J. F., & Brown, L. C. (2017). Projected climate change effects on subsurface drainage and the performance of controlled drainage in the Western Lake Erie Basin. Journal of Soil and Water Conservation, 72(3), 240–250. https://doi.org/10.2489/jswc.72.3.240 Petersen, L. (2019). Impact of Climate Change on Twenty-First Century Crop Yields in the U.S. Climate, 7(3), 40. https://doi.org/10.3390/cli7030040 107 Pryor, S. C., Scavia, D., Downer, C., Gaden, M., Iverson, L., Nordstrom, R., Patz, J., & Roberson, G. P. (2014). Chapter 18: Midwest. In Climate Change Impacts in the United States: The Third National Climate Assessment. Robertson, A. W., Lall, U., Zebiak, S. E., & Goddard, L. (2004). Improved Combination of Multiple Atmospheric GCM Ensembles for Seasonal Prediction. Monthly Weather Review, 132(12), 2732–2744. https://doi.org/10.1175/MWR2818.1 Salem, G. S. A., Kazama, S., Shahid, S., & Dey, N. C. (2018). Impacts of climate change on groundwater level and irrigation cost in a groundwater dependent irrigated region. Agricultural Water Management, 208, 33–42. https://doi.org/10.1016/j.agwat.2018.06.011 Satchithanantham, S., & Ranjan, R. S. (2015). Evaluation of DRAINMOD for Potato Crop under Cold Conditions in the Canadian Prairies. Transactions of the ASABE, 307–317. https://doi.org/10.13031/trans.58.10300 Schaap, M. G., Leij, F. J., & van Genuchten, M. Th. (2001). rosetta : a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. Journal of Hydrology, 251(3–4), 163–176. https://doi.org/10.1016/S0022-1694(01)00466-8 Shapiro, S., & Wilk, M. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3–4), 591–611. Shokrana, M. S. Bin, Ghane, E., Abdalaal, Y., & Nejadhashemi, A. P. (2023). Predicting the effect of weir management on the discharge of a controlled drainage system in a changing climate. Agricultural Water Management, 289, 108534. https://doi.org/10.1016/j.agwat.2023.108534 Singh, R., Helmers, M. J., Kaleita, A. L., & Takle, E. S. (2009). Potential Impact of Climate Change on Subsurface Drainage in Iowa’s Subsurface Drained Landscapes. Journal of Irrigation and Drainage Engineering, 135(4), 459–466. https://doi.org/10.1061/(ASCE)IR.1943-4774.0000009 Singh, R., Helmers, M. J., & Qi, Z. (2006). Calibration and validation of DRAINMOD to design subsurface drainage systems for Iowa’s tile landscapes. Agricultural Water Management, 85(3), 221–232. https://doi.org/10.1016/j.agwat.2006.05.013 Skaggs, R. W. (1978). A WATER MANAGEMENT MODEL FOR SHALLOW WATER TABLE SOILS. Skaggs, R. W. (1980). Drainmod: reference report; methods for design and evaluation of drainage-water management systems for soils with high water tables. . In Soil Conserv. Service. 108 Skaggs, R. W. (1982). Field Evaluation of a Water Management Simulation Model. Transactions of the ASAE, 25(3), 0666–0674. https://doi.org/10.13031/2013.33592 Skaggs, R. W., Brevé, M. A., & Gilliam, J. W. (1994). Hydrologic and water quality impacts of agricultural drainage∗. Critical Reviews in Environmental Science and Technology, 24(1), 1–32. https://doi.org/10.1080/10643389409388459 Skaggs, R. W., & Nassehzadeh-Tabrizi, A. (1986). Design Drainage Rates for Estimating Drain Spacings in North Carolina. American Society of Agricultural Engineers, 29(6), 1631– 1640. Skaggs, R. W., Youssef, M. a., & Chescheir, G. M. (2012). DRAINMOD: Model Use, Calibration, and Validation. Transactions of the ASABE, 55(4), 1509–1522. https://doi.org/10.13031/2013.42259 Sojka, M., Kozłowski, M., Kęsicka, B., Wróżyński, R., Stasik, R., Napierała, M., Jaskuła, J., & Liberacki, D. (2020). The Effect of Climate Change on Controlled Drainage Effectiveness in the Context of Groundwater Dynamics, Surface, and Drainage Outflows. Central-Western Poland Case Study. Agronomy, 10(5), 625. https://doi.org/10.3390/agronomy10050625 Southworth, J., Randolph, J. C., Habeck, M., Doering, O. C., Pfeifer, R. A., Rao, D. G., & Johnston, J. J. (2000). Consequences of future climate change and changing climate variability on maize yields in the midwestern United States. Agriculture, Ecosystems & Environment, 82(1–3), 139–158. https://doi.org/10.1016/S0167-8809(00)00223-1 Tebaldi, C., & Knutti, R. (2007). The use of the multi-model ensemble in probabilistic climate projections. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 365(1857), 2053–2075. https://doi.org/10.1098/rsta.2007.2076 Thomson, A. M., Calvin, K. V., Smith, S. J., Kyle, G. P., Volke, A., Patel, P., Delgado-Arias, S., Bond-Lamberty, B., Wise, M. A., Clarke, L. E., & Edmonds, J. A. (2011). RCP4.5: a pathway for stabilization of radiative forcing by 2100. Climatic Change, 109(1–2), 77– 94. https://doi.org/10.1007/s10584-011-0151-4 USDA National Agricultural Statistics Service. (2024, April). Iowa Cash Corn and Soybean Prices (USDA NASS). Https://Www.Extension.Iastate.Edu/Agdm/. Van Roosmalen, L., Sonnenborg, T. O., & Jensen, K. H. (2009). Impact of climate and land use change on the hydrology of a large‐scale agricultural catchment. Water Resources Research, 45(7). https://doi.org/10.1029/2007WR006760 Wang, X., Mosley, C. T., Frankenberger, J. R., & Kladivko, E. J. (2006). Subsurface drain flow and crop yield predictions for different drain spacings using DRAINMOD. Agricultural Water Management, 79(2), 113–136. https://doi.org/10.1016/j.agwat.2005.02.002 109 Wang, Z., Qi, Z., Xue, L., Bukovsky, M., & Helmers, M. J. (2015). Modeling the impacts of climate change on nitrogen losses and crop yield in a subsurface drained field. Climatic Change, 129(1–2), 323–335. https://doi.org/10.1007/s10584-015-1342-1 Weigel, A. P., Knutti, R., Liniger, M. A., & Appenzeller, C. (2010). Risks of Model Weighting in Multimodel Climate Projections. Journal of Climate, 23(15), 4175–4191. https://doi.org/10.1175/2010JCLI3594.1 Wilcoxon, F. (1992). Individual comparisons by ranking methods. In Breakthroughs in statistics: Methodology and distribution (pp. 196–202). Springer New York. Willison, R. S., Nelson, K. A., Abendroth, L. J., Chighladze, G., Hay, C. H., Jia, X., Kjaersgaard, J., Reinhart, B. D., Strock, J. S., & Wikle, C. K. (2021). Corn yield response to subsurface drainage water recycling in the midwestern United States. Agronomy Journal, 113(2), 1865–1881. https://doi.org/10.1002/agj2.20579 Youssef, M. A., Abdelbaki, A. M., Negm, L. M., Skaggs, R. W., Thorp, K. R., & Jaynes, D. B. (2018). DRAINMOD-simulated performance of controlled drainage across the U.S. Midwest. Agricultural Water Management, 197, 54–66. https://doi.org/10.1016/j.agwat.2017.11.012 Yu, F., Frankenberger, J., Ackerson, J., & Reinhart, B. (2020). Potential Suitability of Subirrigation for Field Crops in the U.S. Midwest. Transactions of the ASABE, 63(5), 1559–1570. https://doi.org/10.13031/trans.13783 Zucker, L. A. , & Brown, L. C. (1998). Agricultural drainage: Water quality impacts and subsurface drainage studies in the Midwest . In Ohio State University Extension (Vol. 871). 110 APPENDIX A: CHAPTER 1 Table A1.1. Sampling strategies in some previous studies investigating P transport through subsurface drainage systems. Study P form Sampling interval Nash et al. (2020) DRP 1 day, event- based1 Compositing scenario 4 aliquots per day (every 6 hours); daily composite Purpose of research Controlled drainage evaluation Carstensen et al. (2019) DRP and TP 7 days – Saadat et al. (2018) DRP and TP 7 to 14 days Pease et al. (2018) DRP and TP 1 day, event- based1 Hourly aliquot; weekly composite 4 aliquots per day (every 6 hours); daily composite Daly et al. (2017) Clement & Steinman (2017) DRP, DP, TRP, and TP 30 days – DRP and TP 30 days and 1 day – King et al. (2016) DRP and TP 1 day King et al. (2015) DP and TP 1 h and 1 day 4 aliquots per day (every 6 hours); daily composite 4 aliquots per day (every 6 hours); daily composite Daigh et al. (2015) DRP and TRP Twice weekly – Hoffmann et al. (2020) TP 7 days 8 aliquots per day (every 3 hours); weekly composite Stamm et al. (1998) DRP 30 min, event- based 2 aliquots in 30 min (every 15 min) Controlled drainage evaluation Controlled drainage evaluation Quantifying the impact of agricultural crop production on surface and subsurface water quality P dynamics in surface and subsurface soils Investigating the role of the subsurface drainage system in P transport Quantifying the impact of FGD gypsum on P concentration and load in tile discharge Investigation of the relationship between P transport and macropore flow Evaluating subsurface drainage TRP concentration and yield in a drained field Evaluating re-established riparian wetlands to mitigate nutrient loss Investigating the effect of preferential flows on P transport TP= total phosphorus; TRP= total reactive phosphorus; DRP= dissolved reactive phosphorus; DP= total dissolved phosphorus. 1- event-based, samples are not collected on a regular schedule but during or after precipitation events. 111 Table A1.2. The relative error in P load estimation and approximate final cost for each sampling strategies over the study period (305 days). Sampling strategy 1 h 3 h 6 h 12 h 24 h 48 h 72 h 7 d 14 d 1 d 2 d 3 d 7 d 1 mm 2 mm 3 mm 5 mm 1 mm 2 mm 3 mm 5 mm Time- proportional discrete Time- proportional composite Flow- proportional discrete Flow- proportional composite Relative error (%) 0 0.2 1.4 5.3 12.2 19.2 27.2 42.7 51 12.4 19.7 27.7 43 0.2 3.1 2.9 5.4 0.5 7.2 4.8 1.9 Number of samples 7320 2440 1220 610 305 154 100 44 22 305 154 100 44 389 190 127 85 66 34 24 16 Cost* Sample analysis Automated sampler Flow sensor and data logger Final $73,200 $24,400 $12,200 $6,100 $3,050 $1,540 $1,000 $440 $220 $3,050 $1,540 $1,000 $440 $3,890 $1,900 $1,270 $850 $660 $340 $240 $160 $5,000 $5,000 $5,000 $5,000 $0 $0 $0 $0 $0 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $5,000 $83,200 $34,400 $22,200 $16,100 $8,050 $6,540 $6,000 $5,440 $5,220 $13,050 $11,540 $11,000 $10,440 $13,890 $11,900 $11,270 $10,850 $10,660 $10,340 $10,240 $10,160 *All sampling strategies included a cost for an automated sampler, except for the 1 to 14-day discrete time-proportional sampling. The cost of $10 per sample was included in the analysis. The cost of a flow sensor was included for all strategies 112 Figure A1.1. Two photos of HydroCycle-PO4, with and without shield (Right Image: courtesy of seabird.com). 113 Figure A1.2. Hourly TRP concentration and drainage discharge in a subset of reference dataset from May 25, 2019, to June 19, 2019. 114 Figure A1.3. Hourly TRP concentration and drainage discharge in a subset of reference dataset from April 7, 2020, to April 21, 2020. 115 Figure A1.4. TRP load (kg/ha) estimation using different sampling intervals (dashed line is the reference TRP load). 116 Figure A1.5. Time-proportional discrete sampling strategy: water samples collected using different sampling intervals. To demonstrate that long sampling intervals often miss sharp increases in P concentration, see storm event on January 24th, 2019. The 24-h discrete sample was taken farthest away from the peak TRP concentration, whereas the 6-h discrete sample was taken closest to the peak. 117 Figure A1.6. Flow-proportional discrete sampling strategy: to demonstrate that long sampling intervals often miss sharp increases in P concentration, see storm event on January 24th, 2019. The 1-mm discrete sample was taken near the peak TRP concentration, thereby better representing the variation in TRP concentration compared to the 5-mm flow interval that missed the peak TRP concentration. 118 REFERENCES Carstensen, M. v., Børgesen, C. D., Ovesen, N. B., Poulsen, J. R., Hvid, S. K., & Kronvang, B. (2019). Controlled drainage as a targeted mitigation measure for nitrogen and phosphorus. Journal of Environmental Quality, 48(3), 677–685. https://doi.org/10.2134/jeq2018.11.0393 Clement, D. R., & Steinman, A. D. (2017). Phosphorus loading and ecological impacts from agricultural tile drains in a west Michigan watershed. Journal of Great Lakes Research, 43(1), 50–58. https://doi.org/10.1016/j.jglr.2016.10.016 Daigh, A. L. M., Zhou, X., Helmers, M. J., Pederson, C. H., Horton, R., Jarchow, M., & Liebman, M. (2015). Subsurface Drainage Nitrate and Total Reactive Phosphorus Losses in Bioenergy-Based Prairies and Corn Systems. Journal of Environmental Quality, 44(5), 1638–1646. https://doi.org/10.2134/jeq2015.02.0080 Daly, K., Tuohy, P., Peyton, D., Wall, D. P., & Fenton, O. (2017). Field soil and ditch sediment phosphorus dynamics from two artificially drained fields on poorly drained soils. Agricultural Water Management, 192, 115–125. https://doi.org/10.1016/j.agwat.2017.07.005 Hoffmann, C. C., Zak, D., Kronvang, B., Kjaergaard, C., Carstensen, M. V., & Audet, J. (2020). An overview of nutrient transport mitigation measures for improvement of water quality in Denmark. Ecological Engineering, 155(April), 105863. https://doi.org/10.1016/j.ecoleng.2020.105863 King, K. W., Williams, M. R., Dick, W. A., & LaBarge, G. A. (2016). Decreasing phosphorus loss in tile-drained landscapes using flue gas desulfurization gypsum. Journal of Environmental Quality, 45(5), 1722–1730. https://doi.org/10.2134/jeq2016.04.0132 King, K. W., Williams, M. R., Macrae, M. L., Fausey, N. R., Frankenberger, J., Smith, D. R., Kleinman, P. J. A., & Brown, L. C. (2015). Phosphorus transport in agricultural subsurface drainage: A review. Journal of Environmental Quality, 44(2), 467–485. https://doi.org/10.2134/jeq2014.04.0163 Nash, P. R., Singh, G., & Nelson, K. A. (2020). Nutrient loss from floodplain soil with controlled subsurface drainage under forage production. Journal of Environmental Quality, 49(4), 1000–1010. https://doi.org/10.1002/jeq2.20072 Pease, L. A., King, K. W., Williams, M. R., LaBarge, G. A., Duncan, E. W., & Fausey, N. R. (2018). Phosphorus export from artificially drained fields across the Eastern Corn Belt. Journal of Great Lakes Research, 44(1), 43–53. https://doi.org/10.1016/j.jglr.2017.11.009 Saadat, S., Bowling, L., Frankenberger, J., & Kladivko, E. (2018). Nitrate and phosphorus transport through subsurface drains under free and controlled drainage. Water Research, 142, 196–207. https://doi.org/10.1016/j.watres.2018.05.040 119 Stamm, C., Flühler, H., Gächter, R., Leuenberger, J., & Wunderli, H. (1998). Preferential Transport of Phosphorus in Drained Grassland Soils. Journal of Environmental Quality, 27(3), 515–522. https://doi.org/10.2134/jeq1998.00472425002700030006x 120 APPENDIX B: CHAPTER 2 Figure A2.1. Two photos of HydroCycle-PO4, with and without shield (Right Image credit: seabird.com). 121 Figure A2.2. Due to the strong relationship between flow and P concentration, we considered the entire flow as the flow event. Hence, if we had subtracted the baseflow from the event flow, the resulting baseflow would not have represented the typical baseflow prior to the event because of the high P concentration. 122 Figure A2.3. The TRP concentration in drainage discharge during summer was more scattered (R- square = 0.03) compared to other seasons. The reason could be a combination of differences in frequency, intensity, and duration of storm events such that during intense and short-duration storm events most likely only P molecules in preferential flow pathways were transported to the monitoring point while in low-intensity and long-duration storm events, water had the opportunity of moving through both preferential flow pathways and soil matrix, thereby washing P molecules in both pathways. Therefore, for the same quantity of drainage discharge, different amounts of P concentrations were measured. 123 Figure A2.4. An example of different time intervals for time proportional discrete sampling, which shows that the 6-h interval is the most suitable sampling interval for capturing the rapid variation in P concentration. The 24-h discrete sample was taken farthest away from the peak TRP concentration, whereas the 6-h discrete sample was taken closest to the peak. 124 Table A2.1. The flow events sorting from the largest to the smallest contribution to P load. The highlighted rows show the events that were not included for rising limb, falling limb, and HI and FI analyses. Date Event Number Event TRP load (kg/ha) Contribution to P loss (%) Water Yield (mm) 11.24 7.39 13.91 9.45 6.52 3.66 7.44 2.94 1.14 6.28 3.70 8.48 17.26 4.31 3.36 3.22 6.96 8.09 4.24 4.10 2.38 4.54 13.40 2.43 16.60 1.45 4.98 12.76 22.29 7.07 23-Jan-19 3-Feb-19 4-Feb-19 6-Feb-19 7-Feb-19 12-Feb-19 14-Feb-19 21-Feb-19 22-Feb-19 23-Feb-19 19-Apr-19 20-Apr-19 30-Apr-19 26-May-19 27-May-19 29-May-19 30-May-19 1-Jun-19 5-Jun-19 10-Jun-19 16-Jun-19 20-Jun-19 26-Oct-19 30-Oct-19 31-Oct-19 18-Nov-19 29-Dec-19 29-Dec-19 11-Jan-20 24-Jan-20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 0.0712 0.0748 3.8609 4.0567 0.090 0.049 0.020 0.025 0.069 0.023 0.010 0.045 0.008 0.028 0.053 0.017 0.013 0.007 0.026 0.029 0.015 0.016 0.011 0.029 0.068 0.012 0.072 0.003 0.027 0.080 0.198 0.035 4.86 2.66 1.11 1.33 3.74 1.26 0.54 2.42 0.44 1.52 2.90 0.94 0.73 0.40 1.39 1.58 0.81 0.86 0.62 1.59 3.71 0.68 3.91 0.19 1.46 4.34 10.74 1.90 125 Table A2.1 (cont’d) 27-Jan-20 19-Mar-20 27-Mar-20 28-Mar-20 7-Apr-20 18-Apr-20 25-Apr-20 15-May-20 17-May-20 10-Jun-20 16-Jul-20 31 32 33 34 35 36 37 38 39 40 41 0.008 0.040 0.031 0.019 0.019 0.002 0.002 0.067 0.212 0.019 0.001 0.44 2.19 1.67 1.05 1.01 0.14 0.10 3.61 11.50 1.04 0.07 2.11 5.95 5.08 4.24 4.02 1.61 1.34 8.46 27.46 2.91 0.81 126 APPENDIX C: CHAPTER 3 Modeling Center Centre National de Recherches Meteorologiques Australian Community Climate and Earth System Simulator Fondazione Centro Euro-Mediterraneo sui Cambiamenti Climatici Beijing Climate Center Climate System Model The Canadian Earth System Model Table A3.1. Summary of CMIP6 projections obtained from GCMs used in NASA NEX-GDDP database for SSP245 emission scenario. Model ID ACCESS-CM2 ACCESS-ESM1-5 BCC-CSM2-MR CanESM5 CMCC-CM2-SR5 CMCC-ESM2 CNRM-CM6-1 CNRM-ESM2-1 EC-Earth3 EC-Earth3-Veg-LR FGOALS-g3 GFDL-CM4 GFDL-ESM4 GISS-E2-1-G INM-CM4-8 INM-CM5-0 MIROC6 MIROC-ES2L Atmosphere and Ocean Research Institute (The University of Tokyo), National Institute for Environmental Studies, and Japan Agency for Marine-Earth Science and Technology Institute Pierre Simon Laplace Climate Modeling Center Korea Institute of Ocean Science and Technology Earth System Model Flexible Global Ocean-Atmosphere-Land System Model European Community Earth System Model Geophysical Fluid Dynamics Laboratory Institute for Numerical Mathematics Goddard Institute of Space Studies IPSL-CM6A-LR KIOST-ESM MPI-ESM1-2-HR MPI-ESM1-2-LR MRI-ESM2-0 NESM3 NorESM2-LM NorESM2-MM TaiESM1 Max Planck Institute for Meteorology The Meteorological Research Institute Earth System Model The NUIST Earth System Model (NESM) version 3 Norwegian Earth System Model Taiwan Earth System Model 127 Table A3.2. The average annual number of dry days and working days during the growing season predicted by DRAINMOD for each GCM. Shallow drain (75 cm) Deep drain (125 cm) GCM ACCESS-CM2 ACCESS-ESM1-5 BCC-CSM2-MR CanESM5 CMCC-CM2-SR5 CMCC-ESM2 CNRM-CM6-1 CNRM-ESM2-1 EC-Earth3 EC-Earth3-Veg-LR FGOALS-g3 GFDL-CM4 GFDL-ESM4 GISS-E2-1-G INM-CM4-8 INM-CM5-0 IPSL-CM6A-LR KIOST-ESM MIROC6 MIROC-ES2L MPI-ESM1-2-HR MPI-ESM1-2-LR MRI-ESM2-0 NESM3 NorESM2-LM NorESM2-MM TaiESM1 Average Working Days 106 108 106 110 106 104 105 105 103 107 107 107 105 108 105 106 109 108 107 108 105 107 105 107 109 103 102 106 Dry Days Working Days Dry Days 61 50 51 51 67 38 30 50 28 45 33 46 33 39 33 37 54 51 43 41 37 33 35 39 52 50 30 43 56 43 45 44 63 30 43 44 41 39 25 39 24 32 25 28 49 44 37 34 31 27 28 34 46 43 22 38 96 94 93 94 102 89 90 86 88 93 93 95 89 95 85 88 99 97 94 88 89 94 93 93 95 89 86 92 128 Figure A3.1. The box plots illustrate the variation in GCM predictions for maximum future temperature (left) and minimum future temperature (right) during 2030–2059. It is evident that the monthly average minimum and maximum temperatures in the future are projected to be higher than in the past. Additionally, the shorter quartiles for both minimum and maximum temperatures in the future indicate increased confidence in the projected temperature increase. 129 Figure A3.2. The 30-year average monthly water-table depth for historical (1994–2023) and future (2030–2059) periods for deep drains under free drainage. 130 Figure A3.3. The 30-year average monthly actual evapotranspiration for historical (1994–2023) and future (2030–2059) periods under free drainage for deep drains. 131 Figure A3.4. The 30-year average monthly drainage discharge for historical (1994–2023) and future (2030–2059) periods under free drainage for deep drains. 132