FILLING IN THE GAPS: MODELING THE ROLE OF GROUNDWATER IN LAKE ERIE’S NUTRIENT BUDGET By Alexis Ann Lanier A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Geological Sciences – Master of Science 2022 ABSTRACT Lake Erie is a hotspot for large harmful algal blooms, which damage human health, degrade natural habitats, and impair industries reliant on the lake. The Maumee River watershed, the largest in the Great Lakes, often acts as a major driver for these blooms, as it is the largest contributor of nutrients to the lake, mainly attributed to intense agricultural activity. Consequently, surficial transport of phosphorus and nitrogen within the Maumee River watershed has been extensively studied. However, there has been very little research into the role of groundwater here, especially groundwater modeling studies. Here, I evaluate the literature that has explored nutrient transport to Lake Erie, with a focus on the Maumee River watershed, and examine groundwater nutrient transport. This knowledge will inform nutrient management decisions, especially those regarding future and legacy nutrient loads. In Chapter 1, I review the current state of literature on hydrologic nutrient modeling in the Lake Erie Basin. I highlight common themes in the literature and detail prominent gaps. Specifically, I focus on the role of groundwater in nutrient modeling studies within the Maumee River watershed and recommend future directions for research. In Chapter 2, I create a spatially explicit, process-based groundwater model of the Maumee River watershed. This model allows me to quantify the contributions of groundwater in the context of total basin loading. I then quantify the role of legacy nutrient accumulation by reducing input loads in a projected future scenario. This research completes the nutrient budget by highlighting ‘hidden’ groundwater nutrient loads and informs the timescale of subsurface nutrient management. This thesis is dedicated to my cats, Nilla and Pizzelle, and to Shay Steele. Thank you for being there for me through it all. iii ACKNOWLEDGMENTS This thesis is the fruit of two years of hard work, done not just by me, but by many supporters along the way as well. It was a very challenging and rewarding experience, and I am grateful for the things I learned as well as those who made these experiences possible. I want to acknowledge and thank the mentors, friends, and family that have accompanied me along this journey. I will begin by thanking my committee, Drs. Dave Hyndman, Anthony Kendall, Sherry Martin, and Bruno Basso. Thank you all for your guidance and patience, your advice and insight, and your mentorship and friendship. I would like to thank my primary advisor, Dr. Dave Hyndman, for encouraging and motivating me. His optimism both reassured and challenged me to reach for goals I had previously thought were too difficult, and his big-picture view helped to widen my perspective on the project. Although we have yet to meet in-person due to restrictions from the COVID-19 pandemic, I hope to one day shake his hand to express my gratitude. I would like to thank Dr. Anthony Kendall for his unwavering support in technical details, conceptual development, and personal ideals. His feedback is always honest and sincere, which is a necessary part of bettering yourself and your work. His knowledge and willingness to learn made him my go-to person to troubleshoot difficult issues and develop new ideas, and I have learned so much from him. His inclusivity and personal ideals are an inspiration, and made the lab a highly welcoming place, even if it was virtual. I am very grateful for his time, guidance, and friendship. I would like to thank Dr. Sherry Martin for her conceptual support, especially with the nutrient side of the model, and for her unfailing kindness. She made a point to check in on me to ensure that I felt welcome and heard, especially early on, which went a long way to making me feel included, respected, and a part of the team. She is a role model, who helped to iv show me the importance of science communication to any audience, and who taught me a great deal. I would also like to thank Dr. Bruno Basso for his support and conceptual contributions to the nutrient model. Thank you to other members of the Department of Earth and Environmental Sciences, especially faculty members who have enriched my experience, and staff members, including Elizabeth McElroy, Ami McMurphy, Pam Robinson, Dallas Coryell, and Judi Smelser, who make this work possible. Thank you to my students, who have encouraged me to continue my work in science education. The members of the Hydrolab are well-deserving of my praise and gratitude as well. They opened their arms to me during a global pandemic, and welcomed me to an environment which helped me grow both academically and personally. Dr. Anthony Kendall, the current leader of the lab, has cultivated a supportive and inclusive space for the students. Previous mentors in the lab, Drs. Dave Hyndman, Sherry Martin, and Leanne Hancock, helped to make this possible. I want to thank my friends and labmates Ally Brady, Chanse Ford, Bailey Hannah, Quercus Hamlin, Brent Heerspink, Alex Kuhl, Ben McCarthy, Behnaz Mirzendehdel, Jeremy Rapp, Jake Stid, and Luwen Wan. Their continuing positivity, support, and encouragement was an essential piece of my journey. Alex, thank you for your guidance and advice. You are a role model and an inspiration to us all. Brent and Quercus, thank you for making a special point to reach out and become friends early on, as well as for your friendship, encouragement, and advice through the end. Ally, thank you for your kindness and steadfastness in standing up for what you believe in. Behnaz and Luwen, thank you for your optimism and constant willingness to bounce ideas around and collaborate. Jake, thank you for your creativity, positivity, and constant sense of camaraderie. Jeremy, thank you for reminding us to take care of ourselves and to have a little v fun along the way. Finally, I would like to thank my family. Thanks to my mom and dad, Kay and Troy Lanier, for their never-ending love and support. You made me who I am today, and I couldn’t have done any of it without you. Mom, thank you riding with me for the stressed or worried or angry or excited phone calls, for visiting me, and for your encouragement to tread my own path. Shay Steele, thank you for everything. You were there with me through it all and supported me all the way. You make me a better person every day, and I couldn’t have done it without you. I cannot describe how grateful I am to have you in my life. Last, but not least, I would like to thank my cats, Nilla and Pizzelle. They were there for me, by my side or in my lap, throughout the entire process. They are a source of comfort and unconditional love. Again, I thank you all for your support. vi TABLE OF CONTENTS LIST OF TABLES .................................................................................................................. viii LIST OF FIGURES .................................................................................................................. ix LIST OF ABBREVIATIONS......................................................................................................x CHAPTER 1: REVIEW OF HYDROLOGIC NUTRIENT MODELING IN THE LAKE ERIE BASIN ........................................................................................................................................1 1. Abstract...............................................................................................................................1 2. Introduction.........................................................................................................................2 3. Background .........................................................................................................................4 4. Methods ............................................................................................................................ 13 5. Results .............................................................................................................................. 14 6. Discussion ......................................................................................................................... 25 7. Conclusions ....................................................................................................................... 30 CHAPTER 2: QUANTIFYING THE HIDDEN COMPONENT OF NUTRIENT BUDGETS: MODELING GROUNDWATER N AND P LEGACIES IN A MAJOR LAKE ERIE WATERSHED .......................................................................................................................... 32 1. Abstract............................................................................................................................. 32 2. Introduction....................................................................................................................... 33 3. Methods ............................................................................................................................ 36 4. Results .............................................................................................................................. 61 5. Discussion ......................................................................................................................... 69 6. Conclusions ....................................................................................................................... 74 7. Acknowledgments ............................................................................................................. 75 REFERENCES ......................................................................................................................... 76 APPENDIX .............................................................................................................................. 91 vii LIST OF TABLES Table 1.1. Study locations in literature review. .......................................................................... 15 Table 1.2. Nutrient forms in literature review. ........................................................................... 15 Table 2.1. Summary of model hydraulic conductivities. ............................................................ 48 Table 2.2. Drain spacing parameters. ......................................................................................... 54 Table 2.3. Summary of model hydraulic conductivities. ............................................................ 54 Table A1. Lake Erie articles in Chapter 1. ................................................................................. 91 viii LIST OF FIGURES Figure 1.1. Groundwater relevance of Maumee nutrient modeling studies. ................................ 21 Figure 2.1. Map of the Maumee River Watershed...................................................................... 38 Figure 2.2. Hydrogeologic characteristics of the Maumee River Watershed. .............................. 39 Figure 2.3. Methodology flow diagram. .................................................................................... 42 Figure 2.4. MODPATH groundwatershed delineation. .............................................................. 44 Figure 2.5. Model error by geologic type. .................................................................................. 48 Figure 2.6. Drain types used in model. ...................................................................................... 51 Figure 2.7. Tile drain conductance............................................................................................. 52 Figure 2.8. Surrogate conceptual model..................................................................................... 53 Figure 2.9. Predicting ditch and tile drain conductance. ............................................................. 57 Figure 2.10. Depth to water. ...................................................................................................... 62 Figure 2.11. Flow modeling validation. ..................................................................................... 63 Figure 2.12. Phosphorus concentrations..................................................................................... 64 Figure 2.13. Nitrogen concentrations. ........................................................................................ 66 Figure 2.14. Nutrient validation. ................................................................................................ 68 ix LIST OF ABBREVIATIONS BMP Best management practice DEM Digital elevation model DRP Dissolved reactive phosphorus EPA Environmental Protection Agency gSSURGO Gridded Soil Survey Geographic database LHM Landscape Hydrology Model L-THIA-LID Long-Term Hydrologic Impact Assessment-Low Impact Development MDEQ Michigan Department of Environmental Quality MT Metric ton N Nitrogen NASS National Agricultural Statistics Service NOAA National Oceanic and Atmospheric Administration P Phosphorus SENSE Spatially Explicit Nutrient Source Estimate SPARROW SPAtially Referenced Regression on Watershed Attributes SRP Soluble reactive phosphorus SWAT Soil and Water Assessment Tool TKN Total Kjeldahl nitrogen TN Total nitrogen TP Total phosphorus USDA United States Department of Agriculture USGS United States Geological Survey x WLEEM Western Lake Erie Ecosystem Model xi CHAPTER 1: REVIEW OF HYDROLOGIC NUTRIENT MODELING IN THE LAKE ERIE BASIN 1. Abstract Eutrophication, caused by excess nutrient inputs, is a concern for water quality in freshwater and marine settings across the world, leading to hypoxic dead zones, toxic drinking water, and economic decline for local industry. Lake Erie has undergone an intense eutrophication over the past several decades, making it a prime target for nutrient reduction efforts. Modeling is an essential tool to quantify, predict, and manage loads to the lake, so a diverse variety of models exist for nutrient management in this basin. In this chapter, we compile and assess results from hydrologic nutrient transport models in the Lake Erie Basin, with a focus on the Maumee River watershed, the highest nutrient contributor to the lake. The role of groundwater in Maumee River watershed models is also evaluated. Models were found through an extensive literature search, and we searched articles for groundwater-relevant key words to quantify the extent to which groundwater was considered in each identified Maumee watershed model. Modeling results consistently showed that agriculture was the primary source of nutrients to Lake Erie, and suggested that tile drainage and no-till farming exacerbates nutrient loads. Cover crops and filter strips were recommended, but do not impact subsurface nutrient accumulation and transport. Subsurface nutrient accumulation was estimated to be substantial in areas in the Lake Erie Basin, especially for phosphorus. However, no models were found that explicitly simulate groundwater flow and transport in a process-oriented manner for the Maumee River watershed, nor did most papers consider groundwater. Future research should seek to estimate groundwater nutrient loading to the lake and quantify lag times associated with groundwater transport and legacy nutrient accumulation. 1 2. Introduction Lake Erie is the shallowest of the Great Lakes, containing the smallest volume of water despite having the second smallest area. However, owing to its warm shallow waters, it is naturally the most biologically productive of the Great Lakes, meaning it can support the greatest amount of living material, primarily concentrated in the form of algae (Botts and Krushelnicki, 1995). Early and intense development of the area surrounding Lake Erie further increased the lake’s productivity, leading to nutrient loading and eutrophication. Lake Erie receives the highest annual phosphorus and nitrogen load of all the Great Lakes, driving an aggressive eutrophication of the lake (Robertson et al., 2019a; Hamlin et al., 2020). Eutrophication can cause many problems for life reliant on the lake. Of particular concern in Lake Erie are harmful algal blooms. Algal blooms are common under natural conditions but are generally limited in size by the nutrient (i.e., phosphorus and nitrogen) loads supplied in the water. With eutrophication, extra nitrogen and phosphorus are injected into the system, allowing algal blooms to vastly increase in size. These blooms cause deoxygenation of the underlying water, leading to hypoxic lake conditions, and release the toxin microcystin, which can be harmful to both human health and aquatic life. The Environmental Protection Agency estimates that in the Lake Erie basin, harmful algal blooms pose a threat to the quality of drinking water for about 10 million people relying on Lake Erie for their water supply (US EPA, 2018). At times, the microcystin concentration in drinking water supplies has been so high that warnings were put out to avoid contact with the water and public water supplies were shut off. One such case in Toledo, Ohio during August of 2014 caused about 400,000 people to be without water for nearly three days (Jetoo et al., 2015). Another in September 2013 in Carroll Township caused the 2 water treatment plant to be shut down for 2 weeks, affecting about 2,000 people (Smith et al., 2015). In addition to concerns related to eutrophication, the nutrients themselves, namely nitrate and nitrite, can act as drinking water contaminants. In areas with high agricultural production and high groundwater extraction for use as drinking water, especially where private household wells are common, nitrate and nitrite concentrations can often exceed guideline values (World Health Organization, 2011). Excessive nitrate or nitrite intake can lead to methaemoglobinaemia, which blocks oxygen transport in the blood (World Health Organization, 2011). Nitrate/nitrite contamination is of particular concern for bottle-fed infants, who can develop cyanosis or blue- baby syndrome as a consequence of methaemoglobinaemia (World Health Organization, 2011). To combat declining water quality in the Great Lakes, in 1972 the U.S. and Canada signed the Great Lakes Water Quality Agreement, a commitment to restoring and protecting the Great Lakes. One of the aims of the agreement is to designate and remediate Areas of Concern, which are defined as “geographic area[s]... where significant impairment of beneficial uses has occurred as a result of human activities at the local level” (Gov'ts of Canada and the U.S.A., 2012, p. 21). The Maumee River was designated an Area of Concern in 1987 due to large amounts of agricultural runoff entering the river and moving into Lake Erie, causing excessive phosphorus, nitrogen, and sediment loading to the lake. It is the largest watershed draining into the Great Lakes, with about 90% of this land being used for agriculture (Brown et al., 2020; NASS, 2021). The Maumee River, on average, is the highest single tributary contributor of total phosphorus, soluble reactive phosphorus (SRP), and total nitrogen to Lake Erie (Robertson and Saad, 2011; Maccoux et al., 2016). Because of this, nutrient loading from the Maumee River, 3 especially in the spring, is generally considered to be a driver of algal bloom severity in the western Lake Erie basin (Baker et al., 2014; Kane et al., 2014; Cousino et al., 2015). Due to its importance to Lake Erie water quality, the Maumee River basin has been extensively studied. However, most of these studies focus on the surficial transport of nutrients to and within the Maumee River rather than on subsurface transport (e.g., Bosch et al., 2011; Baker et al., 2014; Kane et al., 2014; Stow et al., 2015; Cousino et al., 2015; Culbertson et al., 2016; Verhamme et al., 2016; Kujawa et al., 2020; Apostel et al., 2021). Most studies that examine the groundwater transport of nutrients tend to use direct measurements and mass balance approaches to quantify groundwater transport of nutrients rather than modeling (e.g., Knights et al., 2017; Oldfield et al., 2020; Rakhimbekova et al., 2021). Modeling allows more flexibility to test how various hydrogeological parameters facilitate the transport of bioavailable nutrients and provide predictions for how nutrient transport might change under future conditions. This chapter aims to assess the current state of literature on hydrologic nutrient modeling in the Lake Erie Basin, placing particular emphasis on the Maumee River watershed and the role of groundwater in those studies. 3. Background 3.1 Nutrient Inputs While both nitrogen and phosphorus are important to eutrophication, the limiting nutrient in a system depends on the typical nitrogen to phosphorus ratio. If this ratio is high, the limiting nutrient tends to be phosphorus, while if this ratio is low, nitrogen tends to be the limiting nutrient. In Lake Erie, the balance of nitrogen and phosphorus has shifted due to an increase in nitrogen loading concurrent with a decrease in the loading of bioavailable phosphorus (Hartig and Wallen, 1984). This shift has caused the primary limiting nutrient in the lake to change from 4 nitrogen to phosphorus. Bioavailable phosphorus loading from the Maumee River Basin has been shown in several studies to be closely interconnected with phytoplankton and cyanobacterial biomass in Lake Erie (Logan, 1987; Baker et al., 2014; Kane et al., 2014; Newell et al., 2019). However, the nitrogen-phosphorus balance in Lake Erie has also been shown to be delicate and sensitive to change. Should this balance shift again, the limiting nutrient could revert to nitrogen. Chaffin et al. (2014) suggest that nitrogen limitation already plays a role in later parts of the algal growth season in Lake Erie. In addition, a prolonged shift in nitrogen to phosphorus ratios could influence which taxa live in the lake (Hartig and Wallen, 1984; Jankowiak et al., 2019; Mooney et al., 2020). Newell et al. (2019) has shown that reduced forms of nitrogen are necessary for toxin production in non-nitrogen fixing cyanobacteria such as those prevalent in western Lake Erie. The proportion of reduced nitrogen forms within total nitrogen loading from the Maumee has increased concurrent with the re-eutrophication of Lake Erie, and this proportion is significantly correlated with cyanobacterial bloom biomass in western Lake Erie (Newell et al., 2019). Even if a shift in taxa or the nitrogen-phosphorus balance does not occur, nitrogen control is important to maintain the biological integrity of stream ecosystems in the Midwest and to meet drinking water quality standards (Logan, 1987; Robertson and Saad, 2011). Nitrogen has been shown to drive invasive plant growth in the Great Lakes Basin (Hannah et al., 2020). There is ample evidence to support the importance of controlling both nitrogen and phosphorus inputs to protect Great Lakes water quality. 3.1.1 Nutrient forms Not all nutrient inputs to Lake Erie can be used by phytoplankton. Nutrient forms that are readily available for use by phytoplankton are referred to as bioavailable. It is the relative concentrations of these forms that primarily drive changes in algal growth rather than changes in 5 total nutrient loads (Logan, 1987; Baker et al., 2014; Maccoux et al., 2016; Newell et al., 2019). For example, re-eutrophication of Lake Erie beginning in the 1990s occurred despite a decline in the amount of total phosphorus (Baker et al., 2014). Bioavailable forms of phosphorus are primarily considered to be concentrated in forms of dissolved phosphorus, especially dissolved reactive phosphorus (DRP), also called soluble reactive phosphorus. Some bioavailable forms (i.e., extractable NaOH) make up approximately 20-35% of particulate phosphorus (Sonzogni et al., 1982; Baker et al., 2014). However, not all bioavailable phosphorus, particularly in particulate form, is able to be used by algae immediately, as it must first be dissolved, desorbed, or mineralized, often converting to DRP in the process. The likelihood of these processes occurring depends on several factors, including the DRP concentration of the receiving waters and the location of the particle within the water (Sonzogni et al., 1982). Nitrogen input to Lake Erie is present in many forms, several of which are bioavailable. Despite potential bioavailability in both chemically reduced and oxidized forms of nitrogen, however, nitrate and total nitrogen seem to have little to no correlation with algal bloom biomass in Lake Erie (Kane et al., 2014), while organic nitrogen or total Kjeldahl nitrogen (TKN) is significantly correlated with western Lake Erie algal bloom biomass (Newell et al., 2019). TKN includes organic nitrogen as well as ammonium, which generally includes all nitrogen forms except nitrate, nitrite, and atmospheric nitrogen (N2). Forms of TKN are chemically reduced, while non-TKN forms tend to be oxidized (excepting atmospheric nitrogen). Nitrogen is a necessary component of toxin production in many varieties of non-nitrogen fixing cyanobacteria, but cyanobacteria prefer reduced forms of nitrogen inputs (Newell et al., 2019). Ammonium is particularly important, as low levels inhibit toxin production (Newell et al., 2019). However, all forms should be monitored due to their high reactivity. 6 3.2 Nutrient Sources Sources of nitrogen and phosphorus can roughly be divided into point and nonpoint source categories. Point sources include singular locations which lead to nutrient plumes (e.g., wastewater treatment plants), while nonpoint sources are distributed across the landscape (e.g., fertilizer). Early efforts to reduce nutrient loads to Lake Erie addressed both nonpoint and point source loading, with specified goals for point sources in the Great Lakes Water Quality Agreement (Baker et al., 2014). However, due to the relative ease with which point sources are identified and controlled in comparison to nonpoint sources, point source loading has decreased over time concurrent with increases in nonpoint source loading (Baker et al., 2014). Consequently, nonpoint sources comprise the largest proportion of nutrient loading to Lake Erie (Maccoux et al., 2016; Hamlin et al., 2020). Loads from groundwater are typically lumped with nonpoint or tributary sources rather than being explicitly quantified. Without addressing nonpoint source pollution, nutrient reduction efforts will not meet Lake Erie’s water quality goals (Baker et al., 2014; Maccoux et al., 2016). 3.2.1 Agricultural sources In Lake Erie, nonpoint source agricultural runoff and discharge is a major source of both phosphorus and nitrogen (Watson et al., 2016). Chemical fertilizer production, referring to synthetically produced fertilizer, typically contains one or more of the elements nitrogen, phosphorus, and potassium to replenish nutrients lost in the soil from the previous growing season and increase production. Organic fertilizers, including animal manure, treated sewage sludge, plant waste, and animal waste, are also sometimes spread across fields. Often, however, it is difficult for farmers to know how much fertilizer should be applied. To maximize production, farmers err on the side of applying more fertilizer than is needed. As a result, surface 7 runoff from fertilized fields is enriched in nutrients from the overapplication of fertilizer and manure, which can then percolate into the underlying aquifer (Luscz et al., 2015; Robinson, 2015). Extensive tile drainage can reduce the amount of nutrients in groundwater by increasing loading directly to surface waters (Robinson, 2015). In addition to manure applied to agricultural fields, nutrient-rich runoff is common from concentrated manure production at locations with intensive livestock husbandry, such as at concentrated animal feeding operations (Luscz et al., 2015; Robinson, 2015). Manure not applied to cropland is often stored in lagoons designed to treat the waste (Robinson, 2015). These storage locations can become significant point sources of nutrients. Lagoon linings can leak, leading to nutrient infiltration to soil water and groundwater. Large rainfall events can also overwhelm the system, leading to overflow and large amounts of nutrient-rich surface runoff. 3.2.2 Urban sources Cities and other concentrated urban areas can act as a significant source of nutrients, largely due to poor infrastructure maintenance and substantial amounts of impervious surfaces with inadequate drainage. Wastewater treatment systems, generally considered point source pollution, are a major source of urban pollution (Maccoux et al., 2016). Effluent sludge, containing the solid sewage components, is often spread on agricultural fields as fertilizer or deposited in a landfill, where it becomes agricultural or urban runoff respectively (Robertson and Saad, 2011). Treated liquid wastewater effluent can be reused for agricultural irrigation or urban activities such as construction or dust control. However, oftentimes treated wastewater from household septic systems is directly discharged into the subsurface as a distributed grid of point sources, which can act like nonpoint sources due to their density in the landscape (Robinson, 2015; Oldfield et al., 2020; Rakhimbekova et al., 2021). Consequently, in rural areas, heavy 8 reliance on domestic septic systems can lead to significant nutrient input to groundwater, especially if septic systems are leaky or poorly maintained (Robinson, 2015). Other urban sources include landfills, industrial emissions, and non-agricultural fertilization. Overland runoff from landfills is rich in chemicals including metals and nutrients, especially nitrogen. Landfill leaching to groundwater is also quite common in landfills built before the enactment of the 1971 Protection of the Environment Act, which required the use of liners and other mitigation techniques (Robinson, 2015). Likewise, improperly stored industrial waste can act as a long-term source of nutrients both to groundwater and to surface waters (Robinson, 2015). Non-agricultural fertilization can be a large portion of total urban nutrient loading, including both lawn and golf course fertilization, and contributing more often to surface runoff than groundwater due to the surrounding impervious cover (Luscz et al., 2015). 3.2.3 Environmental sources Encompassing a large range of primarily nonpoint sources, environmental sources are those sourced in, or delivered via, natural processes. Tributary discharge is by far the largest of these sources, exacerbated greatly by both agricultural and urban runoff (Robertson and Saad, 2011; Maccoux et al., 2016). Tributary discharge, consequently, is difficult to separate from other sources. Additionally, due to the nonpoint source nature of overland runoff and large drainage areas of some Lake Erie tributaries, sourcing tributary loads is difficult to do (Robertson and Saad, 2011). Instead, many studies simply opt to account for tributary loads by watershed. Atmospheric deposition and nitrogen fixation can also act as significant sources of nutrients. Both nitrogen and phosphorus can be delivered to landscapes via wet atmospheric deposition, whereby nutrients are added to the landscape by precipitation, or via dry deposition, 9 whereby nutrients are transported on the wind by dust particles that are deposited on the land surface (Luscz et al., 2015). Atmospheric deposition delivers a much higher proportion of total nitrogen loads to the landscape than of total phosphorus loads (Luscz et al., 2015). Another pathway for nitrogen to enter the landscape is through fixation, which occurs when lightning strikes or when bacteria or legumes transform atmospheric nitrogen into organic nitrogen (Robertson and Saad, 2011; Lewandowski et al., 2015). Legacy nutrients, or nutrients present in the subsurface, is a source which is often overlooked. Even with reduced nutrient input, there is a large reservoir of phosphorus and nitrogen in the groundwater, plus additional mass can desorb from soil grains and enter the aqueous phase to be transported by groundwater or surface water to lakes, especially in wet conditions (Sharpley et al., 2013; King et al., 2017; Van Meter et al., 2021). Groundwater carrying these nutrients can either be directly discharged into Lake Erie, or indirectly discharged through baseflow to rivers and streams that discharge into Lake Erie (Robinson, 2015). Local biota in the lake have complex relationships with algal blooms and nutrient dynamics. For Lake Erie in particular, dreissenid mussels graze algae, reducing populations present in blooms. However, they also excrete significant amounts of nitrogen and phosphorus into the water column, fueling algal growth (Zhang et al., 2008). In addition, increased competition with other algal species may decrease the size or toxicity of blooms (Zhang et al., 2008). Competition is in large part determined by the relative amounts of nitrogen and phosphorus in the lake, with some ratios favoring more toxic species (Newell et al., 2019). 3.2.4 Subsurface nutrient mobility While nutrients may be present in groundwater, they do not always make it to Lake Erie due to reduced subsurface mobility and attenuation. Phosphorus is often thought of as largely 10 immobile in groundwater, especially in the vadose zone, due to its ability to readily sorb to soils in oxic conditions (Lewandowski et al., 2015). However, all available sorption sites can be filled with consistently high phosphorus loads, especially in silicon-dominant sediments, enhancing phosphorus mobility (Lewandowski et al., 2015). In the saturated zone, reducing redox conditions lower the sorption capacity for phosphorus, thus increasing the amount of dissolved phosphorus (Lewandowski et al., 2015). Phosphorus sorption is also diminished with increasing pH and decreasing oxygen content, enhancing phosphorus mobility in groundwater (Lewandowski et al., 2015; Robertson et al., 2019b). While sorption is the primary means by which phosphorus is removed from groundwater, mineralization to apatite in the presence of calcium can lead to more permanent phosphorus removal (Robertson et al., 2019b). The mobility of nitrogen is much more variable due to its variety of forms in differing redox states. Nitrate (NO3-), one of the most common contaminants in groundwater, is highly mobile in oxic aquifers and easily flushed from rivers and streams with precipitation (Lewandowski et al., 2015; King et al., 2017). Under oxic conditions, nitrate can travel through groundwater with minimal sorption to soil or other aquifer materials, but is degraded in anaerobic aquifers through denitrification, primarily in upper layers where dissolved organic carbon is in highest abundance to act as an electron donor (Lewandowski et al., 2015). Nitrite (NO2-) acts similarly to nitrate, and the two are often considered together, although nitrite is usually present at lower concentrations. Ammonium (NH4+) is significantly less mobile in groundwater due to its ability to sorb to clays and other cation exchangers in aquifers; it is also much less prevalent in groundwater than nitrate, due in part to its oxidation to nitrate in aerobic conditions via nitrification and its conversion to atmospheric dinitrogen (N2) gas in anoxic conditions via microbial oxidation (Lewandowski et al., 2015). Dissolved organic nitrogen, 11 understudied and not well understood, may come from some combination of natural and anthropogenic sources in relatively high concentrations (Kroeger et al., 2006; Lewandowski et al., 2015). Some studies suggest that dissolved organic carbon, of which dissolved organic nitrogen is a type, experiences substantial losses in the vadose zone due to mineralization, sorption, and microbial transformation (Cronan and Aiken, 1985; Pabich et al., 2001). Although it is suspected that dissolved organic nitrogen plays a larger role than previously thought in lacustrine nutrient budgets, nitrate present in oxic aquifers is generally of the greatest concern due to its high mobility. The temporary removal of phosphorus from groundwater by sorption leaves significant legacy phosphorus mass in the system, while the high mobility of nitrate is typically thought to leave little nitrogen remaining in subsurface storage (Sharpley et al., 2013; Robinson, 2015; King et al., 2017; Van Meter et al., 2021). Evidence of phosphorus accumulation in the subsurface can be found in many studies, such as Muenich et al. (2016), King et al. (2017), Robertson et al. (2019b), and Van Meter et al. (2021). However, nitrogen accumulation in the subsurface is more contested in literature. Some sources disregard nitrogen accumulation, often using the mobility and dominance of nitrate to contrast with phosphorus. For example, King et al. (2017) shows evidence of nitrate flushing in the Maumee and Portage Rivers to contrast it with dissolved reactive phosphorus concentrations, which have little to no trend. Van Meter et al. (2016), however, argues that accumulation of organic nitrogen in the soil constitutes a major portion of the ‘missing’ nitrogen component observed in nitrogen mass balance analyses. Wang et al. (2013, 2016) also predict delays in nitrogen delivery due to subsurface storage. 12 4. Methods The purpose of this review is to analyze previous efforts to model nutrient transport in the Lake Erie Basin with a focus on the Maumee River Watershed. I will assess the modeled efforts to explain the cause of lake re-eutrophication, especially models created for the Maumee River Watershed, and detail knowledge gaps in the literature to provide a foundation for future research and nutrient management. Emphasis is placed on the inclusion of groundwater in hydrologic models, shown in preliminary literature review to be highly understudied. To guide these objectives, I pose several questions: Which nutrient reduction strategies are likely to be most effective in the future? What gaps exist in the current modeling literature surrounding nutrient loads to Lake Erie? Has groundwater been well-characterized and simulated in Lake Erie nutrient modeling? To answer these questions, literature was found through Web of Knowledge. Relevant articles were found by searching key words including “Lake Erie”, “Maumee”, “model”, “nutrient”, “phosphorus”, “nitrogen”, “nitrate”, “groundwater” and “surface water” in various combinations. Literature was analyzed in two phases. First, the literature was examined for any quantitative hydrologic nutrient model within the Lake Erie Basin. Each article was individually assessed for its relevance, and articles that did not use modeling or that did not quantitatively assess nutrient loading to the lake were excluded from study. The remaining papers were then individually examined to determine study area, model type, nutrient studied, and conclusions. I compiled the articles, placing particular emphasis on quantitative conclusions of the paper, for example estimated loading rates and relative loading proportions or changes. For the second phase, I quantitatively analyzed the number of existing groundwater models of the Maumee with searches for “Maumee” and “model” with “nutrient”, “phosphorus”, or “nitrogen” in Web of 13 Knowledge by topic. Resulting papers were not filtered for relevance; all resulting articles were considered. I then searched the resulting papers for the words “groundwater” (as well as its variations, including “ground water” and “ground-water”), “baseflow”, and “aquifer” to evaluate their relevance to groundwater. I individually assessed each paper containing these keywords to determine if a hydrologic model was built of the area. The models were also assessed to determine if they successfully accounted for groundwater. For the sake of brevity and the avoidance of repetition, models of the same type (e.g., Soil and Water Assessment Tool or SWAT) were considered together. 5. Results 5.1 Lake Erie Basin Of the 83 articles identified from the Web of Knowledge search on nutrient models in the Lake Erie Basin, 24 were deemed relevant to the topic at hand and analyzed for the study (Appendix). A large number of these studies examined part or all of the Maumee River watershed (Table 1.1). The most common nutrients modeled in the literature, in descending order, were total phosphorus, soluble reactive phosphorus, total nitrogen, and nitrate/nitrite (Table 1.2). Most articles focused on surface water or surface runoff of nutrients, while only 5 articles primarily examined subsurface processes, such as legacy nutrient accumulation in soil and septic systems as a nutrient source. Commonly studied topics included agricultural management and best management practices (BMPs) (10 articles), responses of nutrient loading to climate change (5 articles), and nutrient source attribution (4 articles). 14 Table 1.1. Study locations in literature Table 1.2. Nutrient forms in literature review. Number and proportion of studies review. Number and proportion of examined in each location. Some articles studies examining each nutrient form. studied more than one watershed, so the total Some articles studied more than one number of studies examined here is not nutrient form, so the total number of equivalent to the sum of the total studies shown studies examined here is not equivalent below. Most studies examined all or part of the to the sum of the total studies shown Maumee River watershed. below. Most studies modeled total phosphorus. # of % of Watershed/location studies studies # of % of Nutrient form All of Lake Erie 1 4.2 studies studies All of western Lake Erie 1 4.2 Total P (TP) 17 70.8 basin Total soluble P 3 12.5 All or large part of 3 12.5 (TDP, SP) Canadian Lake Erie basin Soluble reactive P 14 58.3 All of Michigan lower 1 4.2 (SRP, DRP, PO4) peninsula Particulate P 1 4.2 Maumee watershed 9 37.5 Total N (TN) 10 41.7 Grand watershed 3 12.5 Total soluble N 2 8.3 St. Clair-Detroit 4 16.7 (TDN, SN) watershed Dissolved organic N 2 8.3 Sandusky watershed 4 16.7 Total organic N 1 4.2 Raisin watershed 1 4.2 Dissolved 1 4.2 inorganic N Ammonia (NH3) 2 8.3 Nitrate/nitrite 6 25 (NO3/NO2) Articles that examined BMPs and agricultural management revealed a few common themes. On both a field-scale and watershed-scale, tile drains consistently increased nutrient loading in both nitrogen and phosphorus, and in many cases had higher nutrient loads than surface runoff (Nelson et al., 2018; Guo et al., 2020; Apostel et al., 2021). The only exception to this was Apostel et al. (2021), which predicted that field-scale total phosphorus and dissolved reactive phosphorus would have lower loads in tile drains than in surface runoff. However, Guo et al. (2020) found that tile drains transported higher dissolved reactive phosphorus loads than surface runoff. No-till farming techniques, despite their reputed advantages, were also generally 15 found to increase nutrient loads (Forster et al., 2000; Daloğlu et al., 2012; Wang and Boegman, 2021). These increases were independent of nutrient type and spatial scale. However, the effects were dependent on the pathway examined. Forster et al. (2000) found increased nitrate in runoff under no-till farming scenarios, but decreased organic nitrogen and decreased total phosphorus in soil. Increases in nitrate for this study were minimal, 5-14%, compared to decreases in soil organic nitrogen and total phosphorus (44-47% and 45-46% respectively). Daloğlu et al. (2012) found increases in dissolved reactive phosphorus loads from the Sandusky watershed to Lake Erie under no-till scenarios ranging from 20-25%, while conventional tillage led to decreases of 5-20% compared to baseline. Wang and Boegman (2021) found more conservative changes, with an increase of 2.5% in total phosphorus loads and a decrease of 1.54% in phosphate loads. The most effective BMPs in the mitigation of nutrient loading, outside of the avoidance of tile drains and no-till agriculture, tended to include cover crops and filter strips (Merriman et al., 2018; Dagnew et al., 2019; Igras and Creed, 2020; Wang and Boegman, 2021). Improvements in nutrient loads necessarily depended on percent implementation, with higher adoption rates achieving lower loads, and implementation location, with implementation in high source areas being most effective (Merriman et al., 2018; Dagnew et al., 2019; Igras and Creed, 2020; Wang and Boegman, 2021). Modeled responses of nutrient loads to future climate change were highly variable, both within and between studies (Culbertson et al., 2016; Wallace et al., 2017; Mehan et al., 2019; Kujawa et al., 2020; Wang et al., 2021). All nutrient forms examined by these studies were predicted to both increase in some models and decrease in some models. In the Maumee River watershed, soluble phosphorus loads are predicted by some to decrease by as much as 25% by the end of the century (Wallace et al., 2017), but by others to increase by up to 42% (Culbertson 16 et al., 2016; Kujawa et al., 2020). Mehan et al. (2019) provides a range of estimates for soluble phosphorus loads in the Maumee River watershed by the end of the century, varying between a decrease of up to 60% to an increase of up to 75%. Phosphorus delivered via tile drains seems to have some consistency, with models from Mehan et al. (2019) and Wang et al. (2021) both predicting a decrease in soluble phosphorus tile loads (35-60% decrease and 22% decrease predicted respectively). Nitrogen predictions are inconclusive, with most predictions near a net zero change or with extremely large changes. Models from Wallace et al. (2017) predict no significant change in soluble nitrogen, while models from Mehan et al. (2019) predict a decrease in nitrate tile loads ranging from 25-75% and an increase in surface runoff nitrate loads ranging from 50-100% by the end of the century. In terms of total nitrogen loads, Wallace et al. (2017) predicts a range from no significant change to an increase of 20%, while similarly Kujawa et al. (2020) predicts a decrease of 1%, with estimates ranging from a decrease of 31% to an increase of 22%. Nutrient source attribution studies examined a range of areas, sources, and topics, with results being largely incomparable (Luscz et al., 2015; Jarvie et al., 2017; Hu et al., 2019; Liu et al., 2021). However, some overlaps in the results can be found despite differences in study area. Findings tended to show that agricultural land uses generally contributed the highest nutrient loads to the areas studied (Luscz et al., 2015; Liu et al., 2021). Liu et al. (2021) found that agricultural sources, namely manure and inorganic fertilizer, contributed the highest amount of nitrogen to the Grand River watershed (about 23% and 25% respectively), followed by biological fixation (about 35%), atmospheric deposition (about 11%), and domestic waste (about 6%). Similarly, Luscz et al. (2015) found both the highest loading rates and highest nutrient contributions from agriculture in the lower peninsula of Michigan. However, Luscz et al. (2015) 17 found higher contributed proportions for agriculture than Liu et al. (2021), totaling about 60% for nitrogen (43% and 17% for chemical fertilizer and manure respectively) and about 87% for phosphorus (59% and 28% for chemical fertilizer and manure respectively. Outside of atmospheric deposition of nitrogen (about 33%), all other contributions were about 5% or lower. Loading rates from agricultural land uses for both total nitrogen (118 kg per hectare per year) and total phosphorus (17 kg per hectare per year) were over three times higher than the next highest land use loading rates, found in pastureland (36 kg per hectare per year for total nitrogen and 4 kg per hectare per year for total phosphorus) (Luscz et al., 2015). Urban loading rates were next highest (31 and 2 kg per hectare per year for total nitrogen and total phosphorus respectively), followed by the lowest rates in unmanaged land (17 and 0.2 kg per hectare per year for total nitrogen and total phosphorus respectively; Luscz et al., 2015). Hu et al. (2019) studied the urban St. Clair-Detroit River watershed, examining urban nutrient load proportions and load proportions to the Detroit River. In urban land uses, the Detroit wastewater treatment plant was by far the largest source of total phosphorus loads, making up 56% of total loads (Hu et al., 2019). This was followed by other point sources, contributing 25%; runoff, contributing 10%; and combined sewer overflows, contributing 9% (Hu et al., 2019). More generally, Jarvie et al. (2017) investigated the cause of the rapid increase in soluble reactive phosphorus loads in 2002 from the Sandusky, Raisin, and Maumee River watersheds. Jarvie et al. (2017) determined that the primary driver of this stepwise increase was increased soluble reactive phosphorus delivery from sources, with about 65% of the load increase sourced in delivery and only about 35% from increased runoff volumes. Some other articles focused on subsurface nutrient sources (Oldfield et al., 2020; Van Meter et al., 2021; Rakhimbekova et al., 2021), while one model was designed to predict nutrient 18 loads (Thomas et al., 2018) and one quantified the amount of nutrients excreted by invasive mussels (Zhang et al., 2008). Oldfield et al. (2020) and Rakhimbekova et al. (2021) quantified the impact of septic systems on the Canadian side of the Lake Erie basin. Oldfield et al. (2020) determined that septic systems contribute about 1.7-2.3% of total phosphorus loads to the lake from Canada and predict future increases, but these estimates have a relatively large degree of error. Rakhimbekova et al. (2021) expand on this research to look at only septic systems near the shoreline, estimating that septic systems within 100 meters of the Lake Erie shoreline contribute up to 17 metric tons per year of phosphorus and 110 metric tons per year of nitrogen. Instead of examining a particular source, Van Meter et al. (2021) focused on the fate of phosphorus in the subsurface as accumulated legacy phosphorus. Since 1900, Van Meter et al. (2021) estimates that about 95.5% of net phosphorus inputs to the Grand River watershed have been retained in the subsurface, leaving about 4.5% to be exported to downstream waters. Of the phosphorus accumulated in the landscape, about 40% is in the form of soil organic phosphorus, 30% is strongly sorbed in soil, 19% is weakly sorbed, 6% is in reservoirs and riparian zones, 5% is in landfills, and less than 1% has accumulated in groundwater (Van Meter et al., 2021). Outside of the subsurface, Thomas et al. (2018) modeled land use descriptors as predictor variables for various nutrient forms. Thomas et al. (2018) found that wastewater treatment plant population density was the most common best predictor (for dissolved organic nitrogen, dissolved Kjeldahl nitrogen, soluble reactive phosphorus, and total phosphorus), followed by the percent of pervious surfaces in the headwaters (for total dissolved nitrogen and total nitrogen). Top predictors for other nutrient forms were percent of impervious surfaces within the watershed (for ammonia), percent of forage cropland in the entire stream buffer (for nitrate/nitrite), and percent of forage cropland in the stream buffer within 600 meters upstream of the water quality 19 sampling sites (for dissolved phosphorus). Zhang et al. (2008) modeled dreissenid mussel excretion, a source within the lake itself which was not included by any other article examined in this study. Zhang et al. (2008) found that significant portions of the total phosphorus and nitrogen measured in the water column were excreted daily by dreissenid mussels in the western Lake Erie basin, ranging from 19.4% to 55.5% for total phosphorus and from 11.8% to 14.6% for total nitrogen. Proportions remained under 5% for the central and eastern basins (Zhang et al., 2008). 5.2 Maumee River Watershed In this literature review of nutrient models in the Maumee River Watershed, 59 were found using the keywords “Maumee” and “model” with “nutrient”, “phosphorus”, or “nitrogen”. Of the 59 resulting articles, only 20 mentioned any of the groundwater-relevant keywords. Of those 20, four had mentions only within the citations, 14 had one to four mentions, and two had 5 or more mentions of the groundwater-relevant keywords (Figure 1.1). Of the 20 articles which mention key words, eight used one model type, the Soil and Water Assessment Tool (SWAT) (Bosch et al., 2011; Verma et al., 2015; Cousino et al., 2015; Culbertson et al., 2016; Merriman et al., 2018; Arhonditsis et al., 2019; Liu et al., 2019; Kujawa et al., 2020); six used statistical models (including one article with 5 or more mentions) (Sinsabaugh et al., 1997; Stow et al., 2015; Zhang et al., 2016; Sharma et al., 2018; Choquette et al., 2019; Fang et al., 2019); three used other hydrologic models distinct from SWAT, either surface water models (Yuan et al., 2011; Liu et al., 2017) or crop models (Gunn et al., 2018); and three did not utilize a model, including one article with five or more mentions (Calhoun et al., 2002; Sinsabaugh and Shah, 2010; River and Richardson, 2019). 20 Figure 1.1. Groundwater relevance of Maumee nutrient modeling studies. Sankey diagram showing how many times Maumee nutrient modeling articles mention groundwater- relevant keywords and what type of models were used. Articles were found through Web of Knowledge. Created using sankeymatic.com. The relatively small number of articles with significant mentions of groundwater-relevant keywords indicate that the majority of previous studies focus on surface water nutrient transport, mentioning groundwater only in passing if at all (e.g., Bosch et al., 2011; Baker et al., 2014; Kane et al., 2014; Stow et al., 2015; Cousino et al., 2015; Culbertson et al., 2016; Maccoux et al., 2016; Verhamme et al., 2016; King et al., 2017; Newell et al., 2019; Kujawa et al., 2020). This suggests that in studies based in the Maumee, groundwater is generally disregarded as a nutrient 21 transport pathway. The large majority of known studies that use modeling to examine nutrient transport in and from the Maumee River Basin, except Verhamme et al. (2016), Robertson and Saad (2011), Robertson et al. (2019a), and Liu et al. (2017), utilize the Soil and Water Assessment Tool (SWAT) for modeling (e.g., Bosch et al., 2011; Verma et al., 2015; Cousino et al., 2015; Culbertson et al., 2016; Merriman et al., 2018; Arhonditsis et al., 2019; Liu et al., 2019; Kujawa et al., 2020). Arhonditsis et al. (2019) claim that over 70% of Lake Erie watershed modelling studies use SWAT. While it is a powerful tool for simulating surface processes, SWAT’s treatment of groundwater is simplistic and highly limited (Shao et al., 2019; Molina-Navarro et al., 2019). The SWAT model is divided into a series of individual Hydrologic Response Units, each of which has lumped parameters and effective no-flow groundwater boundaries. The use of a lumped parameter model prevents the characterization of intermixing and assumes constant groundwater travel times within each Hydrologic Response Unit. SWAT uses a linear reservoir model to characterize groundwater flow, which means that the relationship between change in storage and discharge is assumed to be constant (Neitsch et al., 2011). This simplistic assumption is violated anytime the storage-discharge relationship constant changes in the system, causing the groundwater portion of the model to fail. This constant must also be calibrated for each Hydrologic Response Unit based on the surface runoff response rather than true groundwater storage, and thus the model is unable to output groundwater head. The combination of lumped parameter modeling and linear reservoir modeling leads to a conflict in model resolution and accuracy. Models with higher ‘resolution’, or smaller Hydrologic Response Units, can more finely characterize spatial variability, but also tend to have less accuracy and reliability in groundwater results due to there being less points per Hydrologic Response Unit used to 22 calibrate the storage-discharge relationship constant. Additionally, SWAT treats groundwater as two layers with lumped parameters: an unconfined shallow aquifer layer and a confined deep aquifer layer (Neitsch et al., 2011; Shao et al., 2019; Molina-Navarro et al., 2019). Water is permitted to flow from the unconfined layer into the confined layer, but is not permitted to flow back out and is assumed to discharge somewhere outside of the modeled area (Neitsch et al., 2011; Shao et al., 2019). Consequently, SWAT is unable to accurately characterize relationships between surface water and groundwater in areas with high groundwater contributions, such as the Maumee River watershed (Shao et al., 2019; Molina-Navarro et al., 2019). This leads to underestimation of groundwater flow and, accordingly, groundwater nutrient contributions (Shao et al., 2019). Other models used to characterize the Maumee River watershed are the Spatially Referenced Regression on Watershed Attributes (SPARROW) model (Robertson and Saad, 2011; Robertson et al., 2019a), the Western Lake Erie Ecosystem Model (WLEEM) (Verhamme et al., 2016), and the Long-Term Hydrologic Impact Assessment-Low Impact Development (L- THIA-LID) model (Liu et al., 2017). Again, although powerful, these models are unable to characterize the complexities of groundwater flow. SPARROW is a surface water model, focusing on the transport of nutrients in streams, rivers, and lakes, which by nature does not include a groundwater modeling component (Robertson and Saad, 2011; Robertson et al., 2019a). The WLEEM model developed by Verhamme et al. (2016) couples four individual sub- models to create a surface sediment transport model, able to link algal biomass to nutrient influxes. However, like SPARROW, the model does not incorporate a groundwater flow component. The L-THIA-LID model used in Liu et al. (2017), like the others, is a surface water model that does not characterize groundwater interaction. The models that have been utilized in 23 the Maumee River watershed to date are insufficient to understand groundwater’s contribution to Lake Erie nutrient inflows. Those studies on nutrient transport in the Maumee River watershed that do not utilize hydrologic modeling, primarily using direct measurements/observations, mass balance approaches, or geochemical analyses, also largely overlook the groundwater component (e.g., Baker et al., 2014; Kane et al., 2014; Stow et al., 2015; Maccoux et al., 2016; King et al., 2017; River and Richardson, 2019; Newell et al., 2019). Broadly, these studies focus on nutrient bioavailability and surficial input of nutrients to Lake Erie, aiming to characterize nutrient transport trends and describe how chemical form affects bioavailability. However, neglected groundwater-surface water interactions have a significant effect on trends and complicate bioavailability, adding additional layers of complexity to an already intricate system. Without this groundwater component, a major piece of the system is missing, leaving emergent properties to remain unseen. The only study found which examines groundwater transport of nutrients in the Maumee River watershed considers the whole Great Lakes Basin and does not estimate nutrient loads from the watershed (Knights et al., 2017). Knights et al. (2017) uses a water budget approach to estimate groundwater discharge, validated by field measurements 12 km away from the point at which the Maumee River discharges into the lake, and land use data to assess vulnerability to nutrient loading via groundwater. Modeling can be used to estimate missing groundwater nutrient loads, in addition to assessing lag times associated with legacy nutrient build-up in the subsurface. Studies in other areas have found groundwater to be a significant component of nutrient transport to lakes (Brock et al., 1982; Nakayama and Watanabe, 2008; Lewandowski et al., 2015; Meinikmann et al., 2015; Schilling et al., 2016), so it is reasonable to assume that significant 24 groundwater nutrient transport to Lake Erie is possible. Schilling et al. (2016) used land cover, groundwater quality measurements, and mass balance equations to estimate groundwater loading of nitrate and phosphorus to West Lake Okoboji in Iowa, concluding that groundwater could be responsible for as much as 80% of total nitrate in the lake and 10% of phosphorus in the lake. Meinikmann et al. (2015) used direct phosphorus measurements and previously acquired volume fluxes to assess groundwater phosphorus fluxes to Lake Arendsee, Germany, concluding that lacustrine groundwater discharge could provide more than 50% of overall external phosphorus loads. Lewandowski et al. (2015) review groundwater nutrient transport to lakes, detailing why groundwater is traditionally disregarded in nutrient budgets, factors affecting nutrient transport in groundwater, and how groundwater can be an important contributor in lake nutrient budgets. Nakayama and Watanabe (2008) build a process-based model, incorporating interactions between surface water and groundwater, which showed that groundwater is an important source of nutrients in the eutrophication of Lake Kasumigaura in Japan. Early studies also recognize the importance of groundwater, such as Brock et al. (1982). This study quantified groundwater nutrient flux into Lake Mendota, Wisconsin using direct measurements and concluded that groundwater is a significant source of nutrient loading to the lake, especially for phosphorus. These studies refute notions that groundwater cannot be a significant source of both phosphorus and nitrogen to freshwater lakes. 6. Discussion Lake Erie nutrient inputs are extensively studied and modeled. However, most of these studies focus on individual watersheds within the Lake Erie basin rather than the basin as a whole. The Maumee River is generally reported as delivering the highest nutrient loads to Lake Erie, and as such is often considered to be a primary driver of eutrophication; this explains its 25 prevalence in Lake Erie nutrient models (e.g., Robertson and Saad, 2011; Baker et al., 2014; Kane et al., 2014; Maccoux et al., 2016; Robertson et al., 2019a). In addition, most models tended to focus on surface nutrient loading, despite the demonstrated importance of groundwater nutrient delivery in other lacustrine settings (e.g., Brock et al., 1982; Nakayama and Watanabe, 2008; Lewandowski et al., 2015; Meinikmann et al., 2015; Schilling et al., 2016). The most common incorporation of subsurface nutrient transport was in the SWAT model. In the Maumee River Watershed, SWAT modeling is insufficient to accurately characterize the contributions of groundwater, both direct and indirect (Shao et al., 2019; Molina-Navarro et al., 2019). Of those articles that did model subsurface nutrient transport, none quantified the indirect contribution of nutrients from groundwater via baseflow or tile drains. Although many studies stressed the importance of tile drains, this importance was almost never placed in context of groundwater. In addition, legacy nutrient storage remains highly understudied, despite nearly 96% of phosphorus inputs and 30% of nitrogen inputs ultimately ending up stored in the subsurface (Van Meter et al., 2021; Liu et al., 2021). The accumulation of legacy nutrients in the subsurface presents a significant setback for efforts to reduce nutrient loading to surface water. The presence of legacy nutrients can lead to a lag time between land use changes (e.g., changes in farming practices) and observable effects in nutrient loads, compounded by already increased travel times in groundwater (Meals et al., 2010; Wang et al., 2013, 2016; Muenich et al., 2016; Van Meter et al., 2016, 2018, 2021; Vero et al., 2018; Martin et al., 2021; Liu et al., 2021). Muenich et al. (2016) used the SWAT model to estimate that even after total cessation of fertilizer use, DRP loads in the Maumee River watershed would take three to five years to reach target loads for spring in average weather conditions, and total phosphorus would take even longer. In the Grand River watershed, Liu et 26 al. (2021) estimated mean travel times for legacy nitrogen in the subsurface ranging from 5 to 34 years. In the Mississippi River Basin, Van Meter et al. (2016) also estimate long residence times and biogeochemical lag times for nitrogen in the subsurface, especially in the form of soil organic nitrogen. If discounted in management plans, these lag times can lead to substantial demotivation in nutrient load reduction efforts, which severely harms chances that harmful algal bloom biomass will decrease. Agricultural studies examining the most effective BMPs for reducing nutrient loads were very common for the agricultural Lake Erie basin. Studies generally found that tile drains increase nutrient loads to the lake, likely due to the immediate export pathway opened up in tile drains (Nelson et al., 2018; Guo et al., 2020; Apostel et al., 2021). Instead of allowing nutrient- laden soil water to infiltrate downwards, encountering attenuation pathways via sorption, mineralization, microbial metabolism, denitrification, etc., the water and nutrients enter the tile drains and are carried directly to surface water bodies, where nutrient export to the lake is much faster. No-till farming also tended to increase nutrient loads, but overall conclusions were mixed (Forster et al., 2000; Daloğlu et al., 2012; Wang and Boegman, 2021). No-till farming, compared to conventional tillage, neglects to incorporate fertilizer into the soil and leaves a crop residue on the soil surface, which reduces soil erosion. Soil stratification leads to a buildup of nutrients at the surface, resulting in increased nutrient loads in runoff (Forster et al., 2000; Daloğlu et al., 2012; Wang and Boegman, 2021). The most effective BMPs were cover crops and filter strips (Merriman et al., 2018; Dagnew et al., 2019; Igras and Creed, 2020; Wang and Boegman, 2021). Filter strips are a buffer of grass or other vegetation surrounding either the source fields or the destination waterbody, and act to reduce runoff velocity and promote infiltration and nutrient uptake. Cover crops are planted on fields in the off-season to reduce soil erosion and increase 27 plant uptake of residual nutrients in the soil. However, these proposed solutions are present only for surface sources and transport pathways, ignoring groundwater and legacy nutrients. Overall, the prediction of changes in nutrient loads due to climate change proves to be difficult using most available models, which provide highly variable and inconsistent results. Future nutrient loading under changing climatic conditions is highly dependent on outside variables, which are in turn reliant on the future decisions of stakeholders. The rate of fertilizer application and use of BMPs was highly influential in the Culbertson et al. (2016) models, while other models did not attempt to predict future changes in nutrient sourcing (Wallace et al., 2017; Mehan et al., 2019; Kujawa et al., 2020; Wang et al., 2021). Kujawa et al. (2020) determined that the primary source of uncertainty in nitrogen modeling came from the climate model, while the primary source of uncertainty in phosphorus modeling was the hydrologic model. Hydrologically, phosphorus pathway partitioning and subsurface phosphorus sorption/desorption are difficult to model and highly dependent on initial conditions, such as soil phosphorus concentrations, which are challenging to obtain sufficient measurements for (Kujawa et al., 2020). Nitrogen transport, on the other hand, is dependent largely on its form, which is in turn determined by redox reactions and climatic conditions, such as O2 and CO2 availability. Despite these uncertainties, authors gave several possible explanations for changing nutrient loads with climate change. Increased runoff and river discharge can be expected to increase nutrient loads delivered to Lake Erie (Culbertson et al., 2016; Wallace et al., 2017; Mehan et al., 2019; Kujawa et al., 2020; Wang et al., 2021), as has been seen in previous decades (Jarvie et al., 2017). However, Culbertson et al. (2016) and Wang et al. (2021) predict lower nutrient loads due to increased plant uptake of nutrients from a longer growing season, decreased plant temperature stress, and increased radiation-use efficiency from increased CO2 concentrations. 28 Determining the source of nutrients in the landscape is essential to mitigate future water quality degradation of the Great Lakes. While some increases in nutrient loading are related to direct effects of climate change (Culbertson et al., 2016; Wallace et al., 2017; Mehan et al., 2019; Kujawa et al., 2020; Wang et al., 2021), the primary driver of increasing nutrient loads is the increase of nutrient delivery at the sources (Jarvie et al., 2017). By concentrating efforts to reduce nutrient loading in high-source areas, we can greatly increase the likelihood of reaching nutrient loading targets (Dagnew et al., 2019; Wang and Boegman, 2021). Modeling studies conclusively find that agricultural sources, especially the application of fertilizer, are the primary source of nutrients in the Lake Erie basin (Luscz et al., 2015; Liu et al., 2021). Further research is necessary to determine the optimal amount of fertilizer necessary for crop growth. Ideally, a tool should be developed and distributed to farmers to advise them on how much fertilizer is necessary and to discourage overapplication, accounting for subsurface legacy nutrients. In urban areas, nutrient loading can be decreased by focusing efforts on wastewater effluent. The largest proportion of nutrients in urban areas comes from point sources, especially wastewater treatment plants (Hu et al., 2019). With better management of effluent before it is discharged into surface waters, the urban nutrient footprint can be greatly reduced. In addition, the use of green infrastructure and the reduction of impervious surface area in urban areas can reduce the amount of nutrient-laden runoff entering surface waters (Thomas et al., 2018). However, nutrient sources within the lake itself should not be neglected. Invasive dreissenid mussels, introduced in the late 1980s to early 1990s, excrete significant amounts of ammonia and phosphate into Lake Erie (Zhang et al., 2008). While mussel grazing of algae slightly decreases algal biomass, the excretion of nutrients more than compensates, leading to a net increase of algal biomass due to the presence of the mussels (Zhang et al., 2008). 29 7. Conclusions This chapter compared the results of a variety of nutrient models within the Lake Erie basin, especially within the Maumee River Watershed. Nutrient reduction efforts within the basin are necessary to mitigate eutrophication of the lake, which causes harmful algal blooms and leads to lake water hypoxia and the release of toxic microcystin to the water. Models generally agree that agricultural areas are the primary source of nutrients, and that the use of no-till farming and tile drains exacerbates nutrient loading. Cover crops and filter strips are recommended as best management practices to reduce loads, focusing efforts on high-source areas and encouraging high adoption rates. However, these recommendations from the literature are only surface solutions, which would have little to no impact on legacy nutrients already present in the subsurface. Legacy nutrient accumulation in the subsurface and dreissenid mussel excretion should be considered when setting future nutrient reduction targets, as both will slow and diminish the effects of nutrient load reduction. Predicting future changes in nutrient loading due to climate change is highly uncertain, with hydrologic modeling driving uncertainty in phosphorus loading and climate modeling driving uncertainty in nitrogen modeling. Efforts should be made to improve both hydrologic and climate models to take the inevitable effects of climate change into account when planning nutrient reduction strategies and targets. In addition, modeling of the entire Lake Erie basin is lacking; future research should study the basin as a whole to examine emergent behavior. However, it is also clear that the Maumee River Watershed, as the highest tributary contributor of nutrients to the lake, should be the priority watershed studied in the basin. Subsurface processes are also understudied, with groundwater seldom modeled. No models were found that explicitly simulate groundwater flow and nutrient 30 transport in the Lake Erie basin. The role of groundwater, especially in its interactions with tile drains, should be another target for future studies. 31 CHAPTER 2: QUANTIFYING THE HIDDEN COMPONENT OF NUTRIENT BUDGETS: MODELING GROUNDWATER N AND P LEGACIES IN A MAJOR LAKE ERIE WATERSHED 1. Abstract The Maumee River watershed is the largest source of nutrients to Lake Erie and is a primary driver of harmful algal blooms in the lake. Surficial transport of nutrients in the watershed has been extensively studied, but little is known about the role of groundwater in this system. More widely, the interactions of groundwater nutrient transport with tile drainage are not well characterized or modeled in the literature. To assess the importance of nutrients in groundwater, we built a spatially explicit, process-based groundwater and nutrient transport model using MODFLOW and MT3D. We also incorporate a surrogate model to simulate the farm-scale behavior of tile drainage within the larger model. Results show that 1.5–4.3% of total phosphorus and nitrogen loads from the Maumee River watershed can be attributed to groundwater, with most of these loads being delivered as baseflow to streams (50–66%) or via tile drains (34–48%). If we assume all phosphorus loading from groundwater is in the form of soluble reactive phosphorus (SRP), as much as 18% of SRP loads from the watershed may be attributable to groundwater. Simulations indicate that 264 km2 (1.4%) of the study area exceeds nitrogen concentration limits set by the World Health Organization, putting groundwater-reliant communities at risk. If nutrient leaching to groundwater were to stop today, except from atmospheric deposition, 82–99% of today’s groundwater loads would likely continue loading 100 years into the future due to the compounding effects of long groundwater travel times and biogeochemical legacies. By quantifying the contributions of groundwater, often unintentionally hidden in nutrient source studies as nonpoint or tributary sources, we can inform nutrient management decisions with realistic expectations for future loads. 32 2. Introduction Eutrophication is a major threat to fresh and coastal waters and can lead to serious consequences to human health, natural ecosystems, and industry. For decades, nations across the world have recognized the effects of eutrophication and harmful algal blooms, ranging from deoxygenation and dead zones to toxin release and drinking water contamination, and made efforts to reduce anthropogenic nutrient loading (Nolan and Hitt, 2006; Paerl and Otten, 2013; Watson et al., 2016; Wurtsbaugh et al., 2019). Despite these attempts, many locations still experience severe eutrophication (Nakayama and Watanabe, 2008; Dobiesz et al., 2010; Kane et al., 2014; Watson et al., 2016; Van Meter et al., 2018; Murray et al., 2019), and nutrient loads in some areas still trend upwards (Oelsner and Stets, 2019). Part of the reason for this continued degradation of global waters may be related to the overlooked role of nutrient loading from groundwater (Holman et al., 2008; Lewandowski et al., 2015; McDowell et al., 2015; Rosenberry et al., 2015; Robinson, 2015; Brookfield et al., 2021). Solutions to quantifying the importance of groundwater in the nutrient cycle remain incomplete, with considerable gaps in knowledge surrounding modeling and legacy nutrient stores (Lewandowski et al., 2015; Robinson, 2015; Brookfield et al., 2021). Several excellent reviews and case studies indicate groundwater can act as a significant source of nutrient loading to lakes, major tributaries, and coastal settings (e.g., Brock et al., 1982; Belanger et al., 1985; Oberdorfer et al., 1990; Valiela et al., 1990; Hagerthey and Kerfoot, 1998; Kang et al., 2005; Holman et al., 2008; Nakayama and Watanabe, 2008; Meinikmann et al., 2013, 2015; Kidmose et al., 2015; Lewandowski et al., 2015; McDowell et al., 2015; Rosenberry et al., 2015; Robinson, 2015; Schilling et al., 2016; Martin et al., 2017; Kazmierczak et al., 2020; Förster et al., 2021; Brookfield et al., 2021). However, there is limited work on direct 33 quantification of the importance of groundwater in nutrient transport (Holman et al., 2008; Robinson, 2015). The body of literature on groundwater nutrient transport and retention primarily uses direct measurement and mass balance approaches to approximate groundwater nutrient contribution (e.g., Belanger et al., 1985; Oberdorfer et al., 1990; Hagerthey and Kerfoot, 1998; Kang et al., 2005; Meinikmann et al., 2013, 2015; Schilling et al., 2016; Kazmierczak et al., 2020; Förster et al., 2021). While beneficial and valid, approaches based on direct measurements and mass balance rely on extensive data collection, are highly limited in scale, and are largely unable to predict nutrient loading under alternative scenarios. Modeling allows the flexibility to both quantify groundwater nutrient contributions over large areas and test how changes in inputs and land management practices affect nutrient fluxes. Groundwater nutrient modeling in lacustrine settings is commonly overlooked in favor of surficial modeling. The most common tool to simulate nutrient transport is the Soil and Water Assessment Tool (SWAT) (Wei and Bailey, 2021; e.g, Bosch et al., 2011; Verma et al., 2015; Cousino et al., 2015; Culbertson et al., 2016; Merriman et al., 2018; Arhonditsis et al., 2019; Liu et al., 2019; Kujawa et al., 2020). While SWAT includes groundwater flow in its simulations, this portion of the model is limited and underestimates groundwater nutrient contributions (Molina-Navarro et al., 2019; Shao et al., 2019). Briefly, SWAT uses a lumped parameter model in combination with linear reservoir equations to simulate groundwater flow (Neitsch et al., 2011). This combination prevents intermixing, assumes constant groundwater travel times within lumped units, ignores local variability, and calibrates groundwater storage based on the surface response (Molina-Navarro et al., 2019; Shao et al., 2019). In addition, the model cannot output groundwater head and assumes that the confined aquifer discharges outside of the model domain. Paradoxically, groundwater results for SWAT are less accurate with smaller lumped units due to 34 less calibration points per unit, forcing a compromise between groundwater accuracy and resolution. Another model is needed to accurately assess and quantify the role of groundwater in nutrient transport and retention. Modeling the relationship between groundwater and nutrient transport is complicated by the presence of tile drainage and legacy nutrient accumulation in the subsurface. Tile drainage reduces groundwater travel times, thereby reducing the amount of time nutrients are exposed to subsurface attenuation mechanisms, such as sorption, mineralization, microbial degradation, and denitrification (Schilling et al., 2015). Traditional methods for modeling tile drainage in groundwater simulations involve setting a high conductance value on the drains, allowing large volumes of water to pass quickly out of them. However, this strategy does not reflect farm-scale drain behavior, especially local water table mounding between drains. Consequently, the high- conductance technique can lead to an underestimation of average head in drained cells. In addition, models that do not incorporate time-variant nutrient input from previous decades of nutrient application may neglect to simulate subsurface nutrient accumulation in the subsurface, termed legacy nutrients. Legacy nutrients can act as both a sink and source of nutrients depending on sorption gradients. Accurate understanding of modeled tile drain-groundwater relationships and the magnitude and timescale of legacy nutrient retention is needed in the literature. The western Lake Erie basin is an ideal location to simulate groundwater nutrient loading due to its history of harmful algal blooms and rich data record (Watson et al., 2016). Lake Erie, as part of one of the largest collections of freshwater in the world, is an important resource. However, Lake Erie has a long history of biological productivity, exacerbated by large amounts of nutrient-rich runoff from intense agricultural activity in its surrounding watershed (Botts and 35 Krushelnicki, 1995). We still fail to consistently meet nutrient reduction targets, such as the phosphorus loading target of 11,000 metric tons annually set in 1972 by the Great Lakes Water Quality Agreement, exceeded in 2007 and 2011 (Gov'ts of Canada and the U.S.A., 2012; Maccoux et al., 2016). Moreover, phosphorus loads do not show a consistent decrease (Stow et al., 2015; Maccoux et al., 2016; Oelsner and Stets, 2019), with increases evident in bioavailable forms of phosphorus and nitrogen (Kane et al., 2014; Stow et al., 2015; Watson et al., 2016; Newell et al., 2019). Given these trends, it comes as no surprise that harmful algal bloom frequency and severity show increasing trends (Kane et al., 2014; Watson et al., 2016; Newell et al., 2019; Sayers et al., 2019). Discharge from the Maumee River, the highest tributary contributor of nutrients to Lake Erie, often drives these eutrophic events (Robertson and Saad, 2011; Baker et al., 2014; Kane et al., 2014; Cousino et al., 2015; Maccoux et al., 2016). Here, we develop a minimally calibrated, spatially explicit, process-based nutrient transport and groundwater model using MODFLOW and MT3D to quantify the impact of nutrient transport and retention in groundwater on Lake Erie management efforts. Particular attention is paid to tile drainage and legacy nutrient accumulation. Our primary objective is to quantify groundwater nutrient contributions from the Maumee River watershed to Lake Erie to more accurately quantify the magnitude and timescale of these contributions related to nutrient load reduction goals. We also assess the magnitude of legacy nutrient accumulation in the Maumee River watershed and the role of tile drainage in groundwater nutrient transport. 3. Methods 3.1 Site Description The Maumee River watershed, the largest in the Laurentian Great Lakes, is located primarily in northwest Ohio, USA, stretching into neighboring states including northeast Indiana 36 and southeast Michigan (Figure 2.1). The Maumee River is the top contributor of nitrogen and phosphorus to Lake Erie (Robertson and Saad, 2011; Maccoux et al., 2016; Verhamme et al., 2016; US EPA, 2018; Newell et al., 2019), providing nearly 50% of total annual phosphorus loads to the western basin, about 94% of which is considered to be nonpoint source (Maccoux et al., 2016). Phosphorus reduction goals are set at 40% of 2008 loads to the lake, with reductions evenly proportioned to each watershed (US EPA, 2018). Algal bloom season in the western Lake Erie Basin generally starts in June and ends in October, with blooms usually reaching their highest extents in late August and early to mid-September (Wynne and Stumpf, 2015). Cyanobacteria in the blooms, primarily the genera Microcystis, Anabaena, and Planktothrix, release a toxin called microcystin as well as other cyanotoxins (World Health Organization, 2011). This nitrogen-rich compound (Chaffin et al., 2014) can cause liver damage and is a possible carcinogen (World Health Organization, 2011). Effects of microcystin can propagate through the food chain, but human contact with microcystin typically occurs through ingestion of affected drinking water, or less commonly through recreational activities in affected waters (World Health Organization, 2011). 37 Figure 2.1. Map of the Maumee River Watershed. The surface watershed is shown in light green and the modeled groundwatershed is in darker green. Major cities with over 30,000 people are located at the red dots. The Maumee River is shown as the dark blue line, while the light blue lines are other major rivers in the surface watershed. The watershed is covered primarily by agricultural farmland with a few large urban areas (Figure 2.2a), which drive its large nutrient loads. Cropland covers approximately 90% of the surface watershed, with the major crop being a corn-soybean rotation (Brown et al., 2020; NASS, 2021). Quaternary geology of the Maumee River watershed consists of glacial deposits from the Last Glacial Maximum deposited during the Wisconsinan, with fine-grained lacustrine deposits dominating the area (Michigan Natural Features Inventory, 1999; Pavey et al., 1999; Gray and Walls, 2002). Coarse lacustrine deposits of sand and gravel, glacial outwash, dune 38 sand, glacial till deposits, and peat and muck also cover parts of the watershed (Figure 2.2b). Bedrock in the area consists largely of Silurian and Devonian limestone and dolomite deposits, with some Silurian through Mississippian shale, sand, and evaporite deposits in the northwest and south of the study area (Figure 2.2c). Limestone and dolomite deposits tend to be karstic and/or fractured in northwest Ohio and western Indiana, resulting in a highly productive aquifer unit (Ohio EPA, 2014). Recharge is estimated to range from 0 to 376 cm per year, with highest rates in the Oak Openings region and in upland headwaters. A B Figure 2.2. Hydrogeologic characteristics of the Maumee River Watershed. A) Land cover and land use. Most of the watershed is covered by cropland (orange), with some large urban areas (red). B) Quaternary geology, deposited primarily during the Wisconsinan glaciation. Lacustrine deposits (blue) dominate the watershed, with till in moraine deposits. C) Bedrock geology. The watershed largely has a carbonate limestone bedrock, which is highly karstic in nature and highly productive as an aquifer. Interbedded deposits of shale cover the northwest region. D) Recharge to groundwater averaged over the 2000-2014 time period, in meters per year, obtained from the Landscape Hydrology Model (LHM) (Hyndman et al., 2007). 39 Figure 2.2 (cont’d) C D Before the late 1800s, poor drainage in the watershed combined with post-glacial lake formation resulted in a vast wetland known as the Great Black Swamp along the southern side of the Maumee River (Kaatz, 1955; Mitsch, 2017). The Great Black Swamp is estimated to have had an area of about 4,000 km2, consisting intermittently of forested swamps and marshland (Mitsch, 2017). Standing water was present in the swamp nearly constantly, greatly hindering travel through the area year-round (Kaatz, 1955). First efforts to drain portions of the swamp began in 1823, but did not gain much traction until the late 1800s, when public drainage laws were passed in 1859 and utilization of natural clay deposits fueled tile production through the 1870s (Kaatz, 1955). By 1930, the swamp was nearly totally drained and agricultural production rapidly increased to become “the most completely cropped area in the state (Kaatz, 1955, p. 33)”. 40 3.2 Modeling Methods We used the USGS MODFLOW, MODPATH, and MT3D-USGS software to simulate saturated groundwater flow and transport in a multi-step process illustrated in Figure 2.3. First, a preliminary groundwater model was created to define more robust boundary conditions for the regional Maumee River groundwatershed. Then, we created a series of high-resolution flow models at the farm-scale to directly simulate agricultural tile drainage. From these, we created a statistical surrogate model for tile and county drain conductances that we used in our regional model. With this surrogate model embedded, we then calibrated the groundwater recharge and hydraulic conductivity values in the model by manually applying a multiplier to model-wide recharge, using baseflows as calibration targets, and to hydraulic conductivity zones with a substantial head bias. This final calibrated model was used to provide flows for MT3D-USGS. Finally, we ran transport simulations for total phosphorus and total nitrogen concentrations and loads in groundwater, accounting for phosphorus sorption and legacy inputs. 41 Figure 2.3. Methodology flow diagram. First, an initial groundwater flow model was created using MODFLOW. This was used, in conjunction with MODPATH particle tracking, to delineate a groundwatershed of the study area. Next, we created surrogate farm-scale models using MODFLOW to simulate county and tile drain conductance, correlating conductances to soil hydraulic conductivity. The groundwatershed and conductance equations were then input to an updated MODFLOW model of the groundwatershed, which was calibrated using observed head and baseflows. Finally, the flow model was used to create an MT3D nutrient model to simulate nitrogen and phosphorus concentrations in the domain. 42 3.2.1 Initial groundwater flow model Our saturated groundwater model encompassed the groundwatershed of the Maumee River, conceptually divided into two geologic layers: a shallow unconsolidated glacial drift aquifer, and a deep bedrock aquifer. To delineate the spatial extent of the groundwatershed, we created a preliminary model domain by extending the surface watershed out to the next nearest major river where possible and connecting them with smaller flow lines in the National Hydrography Dataset (Figure 2.4A) (USGS and National Geospatial Program, 2020). A preliminary groundwater flow model was created with MODFLOW-2005, a saturated three- dimensional groundwater flow model developed and maintained by USGS (Harbaugh et al., 2017). The model was developed using the coding language Python and the package FloPy (Bakker et al., 2021), maintained by USGS. Cell size, layer elevations, boundary conditions, hydraulic conductivity, and groundwater recharge were set equal to initial model values (see Table 2.1 and Recharge estimation and calibration). Agricultural drain tiles and adjacent county drainage ditches were simulated using the MODFLOW DRN package with an assumed width and sediment thickness, using hydraulic conductivity values from the uppermost cell to estimate conductance (see equation 1). 3.2.2 Defining the groundwatershed A groundwatershed is the land area within which aquifers will contribute water to the surface watershed, in this case the Maumee River watershed. We used MODPATH (Pollock, 2017), a post-processing particle-tracking software used in conjunction with MODFLOW, to track the ending locations and travel times of 100,000 particles randomly distributed across the model domain. The starting locations for each particle whose tracking ended within the surface watershed were compiled and included in the groundwatershed. We created and dissolved a 500- 43 meter buffer surrounding each starting point to be included, and combined areas covered by these buffers and the surface watershed. We then manually smoothed a few irregularities in the area, such as thin concavities, to avoid isolated sets of cells. This final estimated groundwatershed was used as the final model domain (Figure 2.4A). The travel times for each of those particles to reach their ultimate destination is shown in Figure 2.4B and range from less than a year to over 500 years. A B Figure 2.4. MODPATH groundwatershed delineation. A) Comparison of the model domains used, with the Maumee River surface watershed shown for reference (light green). Preliminary model domain used for delineation of the groundwatershed is shown in dark green. Groundwatershed, delineated using MODPATH (Pollock, 2017), was used as the final model domain (medium green). B) Groundwater travel times, as modeled by MODPATH (Pollock, 2017). Travel times range from less than a year in locations close to an outlet (e.g., tile drain, stream) to over 500 years. 44 3.2.3 Maumee groundwatershed flow model Using the boundary defined from the MODPATH simulations, we then defined the groundwatershed flow model. This model was created using MODFLOW-NWT (Niswonger et al., 2011), a saturated three-dimensional groundwater flow model using a Newton-Raphson formulation of MODFLOW-2005, intended to improve solutions in models involving drying and rewetting of cells in unconfined aquifers. As with the initial groundwater flow model, this model was developed using the coding language Python and the FloPy package (Bakker et al., 2021). All model boundaries except those adjacent to Lake Erie were considered to be no-flow cells. Cells bordering Lake Erie were modeled as constant-head cells, with long-term average lake levels estimated at 174.5 meters (NAVD88) from the NOAA Marblehead, Ohio lake monitoring station (NOAA). The domain was divided into 240-meter by 240-meter cells. Surface elevation was defined using a 1-arc second Digital Elevation Model (DEM) from the USGS National Elevation Dataset (USGS, 2018). Bedrock elevations were obtained from a USGS bedrock topography map (Soller and Garrity, 2018). After resampling to the model grid, the minimum layer thickness was set by lowering the bedrock elevation to 2 meters below the DEM in areas where the Quaternary sediment thickness was less than 2 meters. Because the complete bedrock portion of the groundwater flow system is not well-characterized, we assigned the bedrock layer a constant thickness of 200 meters. Providing for some substantial thickness allowed model convergence in areas with very thin Quaternary sediments overlying the bedrock. While this simplification is not representative of true bedrock geometries, bedrock flow is expected to contribute only a very small portion of Lake Erie nutrient loads. The inclusion of the bedrock aquifer allows for improved simulation in areas of minimal drift thickness. 45 3.2.4 Hydraulic conductivity estimation and calibration In both layers, hydraulic conductivity was estimated in groups based on local geologic classification (Table 2.1) (Geologic Survey Division of the MDEQ, 1987; Michigan Natural Features Inventory, 1999; Pavey et al., 1999; Gray and Walls, 2002; Slucher et al., 2006; Gray et al., 2010). Initial hydraulic conductivity values were only available for wells in Michigan, so this information was used to generalize Quaternary shallow aquifer layer estimates for each geologic type in the Michigan domain (Wellogic, 2005). Patterns in measured hydraulic conductivity values matched patterns in Quaternary geology type (Figure 2.2B, Michigan Natural Features Inventory, 1999; Pavey et al., 1999; Gray and Walls, 2002), so we assumed that these values were representative of the upper shallow aquifer layer and that each geologic type could be characterized by a single value. The data were first transformed by a natural log to approximate a normal distribution. Outliers from this distribution at unrealistic values were eliminated to reduce the impact of bedrock wells included in the dataset. For lacustrine deposits, measurements were highly variable; likely because a default hydraulic conductivity value was assigned by drillers when the wells encountered rock, vastly inflating the number of measurements at unrealistically low values. After these outliers were truncated, the median of log-transformed hydraulic conductivity values for each geologic type available was used as the model’s initial estimate. For geologic types not present in the Michigan portion of the model domain, namely sand, peat and muck, and coarse-to-medium till, no measured values were available to estimate from within the study area. For coarse-to-medium grained till, a value was chosen between the estimated values for medium-grained till and coarse-grained till. For sand and peat/muck, an initial central estimate was chosen as the mid-point between minimum and maximum hydraulic conductivity values within the Michigan lower peninsula model (Martin et al., 2021). 46 In the bedrock aquifer layer, initial hydraulic conductivity values for each geologic type were selected based on published ranges (Freeze and Cherry, 1979; Domenico and Schwartz, 1990). A roughly central value was chosen from the karst and reef limestone category for limestone and dolomite, given its karstic nature (Ohio EPA, 2014). A high value was chosen for shale and sandstone to avoid convergence issues due to dramatic changes in hydraulic conductivity across adjacent cells. For sandstone, the highest literature value within published ranges (Freeze and Cherry, 1979; Domenico and Schwartz, 1990) was used. For shale, we used a value about two orders of magnitude below the next lowest initial estimate. We then manually calibrated hydraulic conductivity values to better match simulated and observed head values, obtained from state well databases and compiled in Hamlin et al. (2020), by calibrating a zone-specific multiplier. A multiplier was applied to geologic types that had a substantial bias in head (Figure 2.5). All values of till in addition to sandstone had low initial simulated heads, so a hydraulic conductivity multiplier of less than one was used. For sandstone and till (except coarse till), a multiplier of 0.008 was used, determined through calibration to observed heads. We optimized a separate multiplier of 0.0001 for coarse till, as heads were still abnormally low at higher multipliers. We estimated vertical anisotropy, or the ratio of horizontal to vertical hydraulic conductivity, using a priori knowledge in conjunction with published values (Chapuis and Gill, 1989). These values were also calibrated using a multiplier of 1.02, calibrated to observed head. 47 Table 2.1. Summary of model hydraulic conductivities. Initial estimates of horizontal hydraulic conductivity, in meters per day, and vertical anisotropy for each geologic type. Estimates were later calibrated based on observed head values (compiled by Hamlin et al., 2020) using hydraulic conductivity multiplier values. Model Category Horizontal K (m/day) Vertical Anisotropy Quaternary Lacustrine (fine) 0.3048 5 * 1.02 Lacustrine (coarse) 3.687114 4 * 1.02 Till (fine) 8.150889 * 0.008 5.5 * 1.02 Till (fine to medium) 8.16617 * 0.008 5 * 1.02 Till (medium) 10.321953 * 0.008 4.7 * 1.02 Till (medium to coarse) 15.068 * 0.008 4 * 1.02 Till (coarse) 15.24 * 0.0001 3.5 * 1.02 Peat and muck 10.5325 1.5 * 1.02 Outwash 20.784177 1.5 * 1.02 Sand 25.4 1 * 1.02 Bedrock Limestone/dolomite 17.28 3 * 1.02 Sandstone 4.32 * 0.008 2 * 1.02 Shale/evaporites 0.001 5 * 1.02 A B Figure 2.5. Model error by geologic type. A) Model error (m) in Quaternary geology before calibration. Simulated heads for all types of till are low, with coarse till heads especially low. B) Model error in Quaternary geology after calibration. A multiplier of 0.008 was applied to the hydraulic conductivity of all till types except coarse, for which a multiplier of 0.0001 was used. C) Model error in bedrock geology before calibration. Simulated heads in sandstone are low. D) Model error in bedrock geology after calibration. A multiplier of 0.008 was applied to the hydraulic conductivity of sandstone. 48 Figure 2.5 (cont’d) C D 3.2.5 Recharge estimation and calibration Variable average annual recharge rates across the study area were obtained from the Landscape Hydrology Model (LHM), averaged over the 2000-2014 time period (Hyndman et al., 2007). LHM separates recharge estimates into wetland and upland deep percolation; however, wetland deep percolation values were initially up to 220 meters of recharge per year. Unreasonable values in the wetland deep percolation dataset, at or above 3 meters per year per cell, were truncated and set to 0 meters per year. We chose to set these cells to 0 as these anomalously high values are likely due to error with initial model estimates. The wetland and upland deep percolation estimates were then summed to get total recharge estimates. Unfortunately, LHM outputs do not cover the entire study area. To extrapolate recharge estimates to the rest of the model domain, land cover data from the National Land Cover Database, using classes matching those used in LHM, were used to extend the dataset (Dewitz, 2019). Recharge estimates from LHM for each land cover class were averaged; these averaged values were applied to cells in each land cover class outside of the LHM domain. Calibration was performed on average annual recharge rates by applying a spatially constant multiplier of 0.18, selected using residuals between simulated and observed baseflows 49 (Read et al., 2017). This was done to calibrate the model to observed flows and nutrient concentrations taken during baseflow conditions, as the steady-state flow model is unable to output seasonal flows. Simulated flows were spatially aggregated for each observation point by delineating the watershed of each point and summing cell-by-cell outflows from drain cells within the watershed, following the methods of Mooney et al. (2020). Observation points with total simulated drain flows of 0, primarily watersheds in dry cells, were excluded from flow comparison. 3.2.6 Tile drain modeling We assumed that all streams, subsurface tile drains, and surface ditches act as groundwater sinks and modeled them as drain cells using the MODFLOW DRN package (Figure 2.6). Drain cells are a boundary condition in which water can leave the system when head exceeds the drain elevation at a rate controlled by the conductance. The assumption that surface drainage in the region is largely gaining (and thus can be well-represented by the DRN package) was supported by initial model runs, in which man-made subsurface drains were not included. In these initial runs, much of the model area was covered by flooded cells, especially surrounding the stream drain cells. Flooding in these cells indicates that groundwater levels are higher than land surface elevations and, consequently, above stream levels, providing evidence that the hydraulic gradient is from groundwater to stream channels. Streams were mapped using the ArcGIS Pro 2.8.1 Flow Accumulation tool, which creates a raster of the number of cells draining into each model cell. Cell size of the flow accumulation raster was equivalent to the DEM (USGS, 2018) cell size (1 arc-second, or approximately 27 meters). Cells with 4,000 or more other cells draining into them were considered as stream cells. This threshold was qualitatively chosen to resemble National Hydrography Dataset stream flowline lengths (USGS and National 50 Geospatial Program, 2020). Stream channel elevation was set equal to DEM elevation (USGS, 2018) at those cells, and streams were given a constant width of 2 meters. Stream conductance was calculated by multiplying initial Quaternary hydraulic conductivity by stream area, then divided by an assumed streambed sediment thickness of 1 meter (Equation 1). Surface ditches were mapped using the National Hydrography Dataset’s canal features (USGS and National Geospatial Program, 2020) and given a constant width of 2 meters. Surface ditch elevation was set at 1.5 meters below the DEM. Subsurface tile drains were mapped using USGS areal coverage data from Nakagaki and Wieczorek (2016). Tile drain elevation was set at 1 meter below the DEM to account for the typical burial depth of tile drains in the region (Smith et al., 2015; Gunn et al., 2018). Figure 2.6. Drain types used in model. Tile drains, which typically drain agricultural fields using subsurface pipes, flow into surface ditches. These ditches are usually dug or cut at the edges of agricultural plots, and drain to nearby streams. To estimate conductance of subsurface tile drains and surface ditches, we built a surrogate model using MODFLOW to simulate farm-scale drain behavior at the regional model scale. Typical groundwater models use a constant high conductance for tile drains, which causes head in drained locations to lie at tile drain elevation. However, this underestimates head, as it 51 does not account for mounding of the water table between tile drains (Figure 2.7). This can have a significant impact on water fluxes to the drains, especially in an area as heavily tile drained as the Maumee River watershed. Figure 2.7. Tile drain conductance. Typical models use an assumed high conductance in tile drains, resulting in average head in drained cells to lie at or very near the drain elevation. However, in reality areally-averaged head in a tile drained area would lie above tile drain elevations due to water table mounding between tile drains. The surrogate model covered a 720-meter by 720-meter area with a constant ground surface elevation to represent an individual agricultural field. Within the field, we placed three ditches at 2 meters below the land surface: one at the top of the field, making up the upper model boundary; one at the bottom of the field, making up the lower model boundary; and one approximately in the center of the field, connected to the top and bottom ditches. We then filled the field with tile drains, each connecting to its closest drainage ditch so that a diamond-shaped pattern was formed (Figure 2.8). Tile drains were given a variable elevation based on distance from the connected drainage ditch, starting at an elevation of 1.5 meters below the land surface at the ditches and increasing elevation by an assumed factor of 1:1000 with increasing distance 52 from the ditch. We placed specified head boundaries on the eastern and western edges of the model. We chose a sample area within the larger Maumee model to place the surrogate model in, from which we obtained recharge estimates. We used the resultant heads from the larger Maumee model (run without drains or ditches) as the specified head boundary inputs. We chose a 1-meter cell size to allow adequate space between each tile drain. Figure 2.8. Surrogate conceptual model. The dark red lines represent tile drains, each of which lead to their nearest surface ditch. Surface ditches are shown in bright red, and lie in the center of the field as well as at the upper and lower boundaries. Left and right boundaries were modeled using specified head cells, with head taken from cells in the initial Maumee model. The surrogate model area covers the area of nine Maumee model grid cells. The centermost cell, highlighted in blue, was used to calculate conductances to avoid impacts from the boundary conditions. All cells in the tile drained area, colored light blue, were used to calculate conductance for the tile drains. The two columns of cells surrounding the surface ditches, colored dark blue, were used to calculate conductance for the surface ditches. We tested several tile drain spacings, corresponding to different values of soil hydraulic conductivity within the larger study area. With all other parameters held constant (Table 2.2), we 53 used the Michigan State University Drain Spacing Calculator to estimate appropriate drain spacings based on input values of hydraulic conductivity (www.egr.msu.edu/bae/water/drainage/ node/192). We estimated soil hydraulic conductivity within the larger study area using the gridded Soil Survey Geographic (gSSURGO) database (Soil Survey Staff, 2019). We then transformed the data using a natural log and chose values within the roughly normal distribution of log values to test (Table 2.3). We used a constant ratio of horizontal to vertical hydraulic conductivity of 5:1 for all simulations. Table 2.2. Drain spacing parameters. Table 2.3. Summary of model hydraulic Parameters used to estimate tile drain spacing conductivities. Saturated soil hydraulic using the Michigan State University Drain conductivity values and corresponding drain Spacing Calculator (www.egr.msu.edu/bae/ spacings used in the surrogate tile drain water/drainage/node/192). Saturated soil model. Drain spacings were determined using hydraulic conductivity was the only parameter the Michigan State University Drain Spacing varied in the Calculator. Calculator (www.egr.msu.edu/bae/water/ drainage/node/192). Parameter Estimated value Saturated hydraulic Drain spacing Drainage coefficient 0.0127 m/day conductivity (m/day) (m) Minimum water table 0.3048 m 0.053 3.05 depth 0.059 3.05 Design drain depth 0.9144 m 0.072 3.66 Depth of restrictive 30.48 m 0.088 4.27 layer 0.107 4.57 Effective radius 0.018034 m 0.118 5.18 0.131 5.49 0.160 6.40 0.195 7.32 0.239 8.53 0.291 9.75 0.322 10.7 0.356 11.3 0.494 14.6 0.531 15.5 0.875 23.5 1.44 35.1 2.38 53.0 3.92 79.9 6.47 119 10.7 176 54 We calculated farm-scale drain and ditch conductance within the surrogate model using equation 1, where C is conductance, K is soil hydraulic conductivity, W is width of the tile drain or ditch, L is the intersecting cell length of the tile drain or ditch, and T is the thickness of the sediments below the drain or ditch. We used a tile drain width of 0.1016 meters, a ditch width of 2 meters, a tile drain sediment thickness of 1 meter, and a ditch sediment thickness of 2 meters. 𝐾∗𝑊∗𝐿 𝐶𝑓𝑎𝑟𝑚 = 1 𝑇 A simulation was run for each tile drain spacing and soil hydraulic conductivity value, resulting in 21 runs. To calculate the tile drain and ditch conductance of each simulation at the scale of the regional Maumee model (240 meter cell size), we used equation 2, where F is flow out of the drain, H̅ is the average head within the cell, and E is the drain elevation. We used the Zone Budget package to separate tile drains and ditches so conductances for each zone could be calculated separately. We defined a 240-meter by 240-meter area, representing a cell in the regional model, and divided it into two zones. One zone contained the single row of cells representing the ditch in the middle of the field, while the other zone contained all other cells representing the tile drained area. Flow from the drains and ditches for each zone was input to equation 2, which used average head from the entire sample area for both ditches and drains. We calculated drain conductance using drain elevation averaged across the whole zone, while ditch conductance simply used the constant input ditch elevation. ∑ 𝐹𝑓𝑎𝑟𝑚 𝐶𝑟𝑒𝑔𝑖𝑜𝑛 = ̅ 2 𝐻𝑓𝑎𝑟𝑚 − 𝐸𝑟𝑒𝑔𝑖𝑜𝑛 Once we established a relationship between conductance and hydraulic conductivity for tile drains and ditches at the regional model scale (Figure 2.9), we used the representative equations for these relationships (equations 3 and 4) to determine conductance based on 55 gSSURGO hydraulic conductivity values (Soil Survey Staff, 2019). However, since the surrogate model assumed total field tile drain coverage, we reduced the calculated conductance value for each cell by the proportion of total drain coverage for that cell (equation 5). For ditches, however, the surrogate model only assumed a 240 m2 area per 240-meter by 240-meter cell, as ditches had a width of 1 meter along the length of the cell. To scale ditch conductance values, then, we multiplied the calculated conductance per cell by the ditch area in that cell divided by the 240 m2 assumed in the surrogate model (equation 6). 0.8184 𝐶𝑡𝑖𝑙𝑒 = 503.5 ∗ 𝐾𝑠𝑜𝑖𝑙 + 122.9 3 2 𝐶𝑑𝑖𝑡𝑐ℎ = −3.757 ∗ 𝐾𝑠𝑜𝑖𝑙 + 114.5 ∗ 𝐾𝑠𝑜𝑖𝑙 + 1.137 4 𝐴𝑡𝑖𝑙𝑒 5 𝐶𝑐𝑒𝑙𝑙 = 𝐶𝑒𝑞𝑛 ∗ 𝐴𝑐𝑒𝑙𝑙 𝐴𝑑𝑖𝑡𝑐ℎ 6 𝐶𝑐𝑒𝑙𝑙 = 𝐶𝑒𝑞𝑛 ∗ 𝐴𝑧𝑜𝑛𝑒 56 Figure 2.9. Predicting ditch and tile drain conductance. Relationship between saturated soil hydraulic conductivity and conductance for tile drains (blue) and surface ditches (red) as modeled by MODFLOW surrogate model. The surrogate model took a field-scale portion of the initial regional Maumee model and varied soil hydraulic conductivity with drain spacing. Conductances were then calculated for tile drains and surface ditches using equation 2. These values were then graphed against hydraulic conductivity to develop a relationship between the two values, used to estimate farm-scale conductance in the regional Maumee model. 3.2.7 Nutrient modeling Using results from the MODFLOW model, we created a nutrient transport and attenuation model with MT3D-USGS (Bedekar et al., 2016). We modeled total phosphorus (TP) and total nitrogen (TN) in groundwater, accounting for advection, dispersion, source-sink mixing, and phosphorus sorption. Initial background concentrations of 0.03 and 0.5 mg/L were used for phosphorus and nitrogen respectively (Dubrovsky and Hamilton, 2010; Welch et al., 2010). Input nutrient concentrations to groundwater were extracted from the Spatially Explicit Nutrient Source Estimate maps (SENSEmap) and flux models (SENSEflux) (Luscz et al., 2015, 2017; Hamlin et al., 2020; Martin et al., 2021). Briefly, the SENSE framework quantifies 57 nitrogen and phosphorus inputs, both from point sources and nonpoint sources (i.e., chemical agricultural fertilizer, manure, chemical non-agricultural fertilizer, septic tanks, nitrogen fixation via legumes, and nitrogen deposition), to the Great Lakes Basin circa 2010 and uses them to statistically model the transport and fate of nutrients in surface water and groundwater. More detailed methods on the SENSE framework are described in Luscz et al., 2015; Luscz et al., 2017; Hamlin et al., 2020; and Martin et al., 2021. Nutrient loads to groundwater modeled by SENSEflux were divided by spatially explicit recharge estimates (scaled using an optimized multiplier) and used as baseline nutrient concentrations in recharge to the MT3D model. To account for the effects of legacy nutrient inputs, we scaled baseline nutrient concentrations by source for each stress period modeled from 1900-2012. Legacy land uses and land covers are a significant predictor of legacy nutrient concentrations (Martin et al., 2017), so most sources were scaled by changes in land use and land cover. For legacy fertilizer and manure inputs, we scaled inputs to groundwater (equation 7) by TP and TN estimates in fertilizer and manure respectively (Falcone, 2021). Stress periods were set based on available data in Falcone (2021), with legacy inputs constant for each stress period and the year given in Falcone (2021) situated in the middle of the stress period (e.g., stress periods of 1948-1951, 1952-1956, 1957- 1963). Inputs for the stress period from 1900-1948 were constant and set equal to earliest available data in 1950. TP inputs from atmospheric deposition were scaled in the same manner by total fertilizer usage, as we assumed that atmospheric TP is primarily sourced in agriculture, and TN was scaled on historic atmospheric deposition data (Hegglin et al., in prep.). Nitrogen fixation was scaled by historical soybean acreage (Haines et al., 2018). We scaled input nutrients from septic tanks by population, linearly interpolated between census years (Manson et al., 2021). 58 𝑠𝑐𝑎𝑙𝑒𝑟𝑙𝑒𝑔𝑎𝑐𝑦 𝑦𝑒𝑎𝑟 𝑙𝑜𝑎𝑑𝑙𝑒𝑔𝑎𝑐𝑦 = 𝑙𝑜𝑎𝑑𝑆𝐸𝑁𝑆𝐸 ∗ 7 𝑠𝑐𝑎𝑙𝑒𝑟2010 To model dispersion, we used MODPATH to calculate the flow lengths of each cell to its destination surface water, used as the longitudinal dispersivity (Pollock, 2017). We used 0.1 as the ratio of horizontal transverse to longitudinal dispersivity, and 0.05 as the ratio of vertical to longitudinal dispersivity. We used published estimates for the ideal molecular diffusion coefficients for both phosphorus and nitrogen (Cussler, 1997) and scaled those estimates using porosity to obtain the effective diffusion coefficients (Ullman and Aller, 1982; Cheng et al., 2014). Aquifer porosity was estimated using gSSURGO soil porosity maps in the lowest horizon for the Quaternary aquifer layer (Soil Survey Staff, 2019), and using published porosity estimates for each geologic type in the bedrock aquifer layer (Freeze and Cherry, 1979). A simple advection model was applied as well, using a Courant number of 1. We also implemented a simple source-sink mixing package, using measured concentrations in Lake Erie (0.158 mg/L for TP and 0.871 mg/L for TN) as constant concentrations in the constant head boundary cells. For TP, this measured concentration came from Rowland et al. (2020). For TN, we used the observation site for baseflow concentration nearest the mouth of the Maumee River (Read et al., 2017). To model phosphorus sorption, we used a Langmuir sorption isotherm. Isotherm parameters were determined through a combination of optimization and published estimates for soils in the Maumee River basin (McCallister and Logan, 1978). Simulated concentrations and loads were compared against samples taken during baseflow conditions, and a linear regression was performed to evaluate goodness-of-fit (Read et al., 2017). Baseflow conditions were defined by the samples taken during baseflow season, July 15 to September 15, and with flows below the 5-year (2008-2012) 30th percentile. Observed loads were aggregated for each observation point 59 by multiplying cell-by-cell drain fluxes by simulated concentrations, then summing loads out of each drain cell contained in the point’s watershed (Mooney et al., 2020). To compare simulated loads to observed loads in baseflow, we also summed upstream seasonally constant point loads, obtained from SENSEmap (Hamlin et al., 2020), in each simulated point’s watershed with simulated aggregate loads. Loads were then divided by aggregated drain fluxes to get the simulated concentrations at each observation point. Observation points with total simulated drain fluxes of 0, primarily those watersheds lying in dry cells, were excluded from comparison. We also removed clear outliers from the observation data, having simulated concentrations or loads (with point loads) 0.5-2 orders of magnitude larger than other observation points. Values ranging from 1000 to 5000 m3/kg were tested as the Langmuir adsorption constant, but the model was found to be relatively insensitive to this parameter, with R2 values for both concentrations and loads varying by only a fraction of a percent regardless of the value of maximum adsorption. A final value of 2000 m3/kg was chosen based on the average of published estimates for soils (McCallister and Logan, 1978). A value slightly below that of soil was chosen, as soils tend to be more adsorptive than aquifer materials. We found modeled results to be much more sensitive to the adsorption maximum, so a range of values from 10-8 to 10-4 mg/mg were used in the final modeled results. Nitrogen was modeled conservatively, assuming no subsurface attenuation, as the dominant form (nitrate, NO3-) is highly mobile in groundwater (Lewandowski et al., 2015). To examine lag times in nutrient delivery, we modeled future concentrations and loads from groundwater with highly reduced nutrient inputs. Our scenario eliminated all future nutrient leaching to groundwater, except from atmospheric deposition, which was assumed to be constant at estimated SENSEflux values. This scenario builds on previous models, which quantify legacy nutrient accumulation and estimate ‘current’ 2012 nutrient concentrations in groundwater. 60 4. Results 4.1 Flow Modeling Groundwater flow model results are consistent with observed static groundwater levels. Root mean squared error of head was 10.8 meters. Groundwater head generally follows patterns in surface elevation, with flow towards the Maumee River and Lake Erie from the north and south. Heads ranged from 173.5 to 426.1 meters. Flooded cells are primarily in stream networks, but a swath of flooded cells also covers an area of low elevation just northwest of the former Great Black Swamp and a patch in the northernmost point of the domain. Depth to water ranges from 0 to 86.72 meters, with shallow aquifers tending to lie in lacustrine deposits (Figure 2.10). Deep aquifers are generally found in glacial moraine deposits, with dry cells in the Quaternary layer most commonly in till areas. Errors in head are slightly heteroskedastic, with a larger range of simulated head values at high head (Figure 2.11). 61 Figure 2.10. Depth to water. Modeled depth to water in the Quaternary aquifer layer in meters, ranging from 0 to 86.72 meters. Deep aquifers are most often found under glacial moraines, with dry cells occurring primarily in areas with till. Flooded cells lie in stream beds, with an additional flooded area in the northernmost point of the domain and one just northwest of the former Great Black Swamp. 62 Figure 2.11. Flow modeling validation. Simulated versus observed head elevation at wells across the study area (Hamlin et al., 2020). Residuals are slightly heteroskedastic, with larger error in areas of high head. 4.2 Nutrient Concentrations Model outputs provide maps of nutrient concentrations in solution and in sorption as well as mass fluxes from groundwater for each simulated time period. The spatial distributions of phosphorus in solution and attached to sediments across the study area are similar for each adsorption capacity tested (Figure 2.12). In the Quaternary aquifer, these patterns appear to be largely dependent on the patterns of input phosphorus concentrations in recharge and groundwater travel time. In simulations with higher adsorption capacities (more sorption), concentrations appear to be more dependent on input concentrations in recharge. However, at high sorption capacities, concentrations are smoothed with little variability. In simulations with 63 lower adsorption capacities (less sorption), concentrations are higher in areas with low groundwater travel times. Patterns in the bedrock aquifer appear to be opposite those of the Quaternary layer for simulations with low adsorption capacity, but similar to patterns in Quaternary concentrations for high adsorption capacity. Bedrock concentrations in sorption were consistently low. With the highest adsorption capacity tested (10-4 mg P/mg matrix), phosphorus concentrations in solution range from 0.0163 to 0.613 mg/L, while concentrations in sorption range from 0.00548 to 0.0395 mg/L. With the lowest adsorption capacity tested (10-8 mg/mg), phosphorus concentrations in solution range from 3.86x10-6 to 7.51 mg/L, while concentrations in sorption range from 1.03x10-9 to 9.30x10-6 mg/L. A B Figure 2.12. Phosphorus concentrations. Spatial distribution of TP concentrations in groundwater solution for the Quaternary aquifer. White areas are dry cells or lie outside of the domain. A) Adsorption capacity of 10-4 mg/mg. Concentrations range from 0.0163 to 0.613 mg/L. Concentrations have low variability and are smoothed across the domain. B) Adsorption capacity of 10-6 mg/mg. Concentrations range from 1.56x10-4 to 7.01 mg/L. C) Adsorption capacity of 10-8 mg/mg. Concentrations range from 3.86x10-6 to 7.51 mg/L. Concentrations generally follow patterns in groundwater travel time, with high concentrations in areas of low groundwater travel times. 64 Figure 2.12 (cont’d) C Patterns in nitrogen distribution are similar to those of phosphorus modeled with low adsorption capacity, again following patterns in groundwater travel time (Figure 2.13). Areas of low groundwater travel time had the highest modeled concentrations in the Quaternary aquifer. In the bedrock aquifer, again concentrations are roughly opposite those in the Quaternary deposits, with lowest concentrations in areas of low groundwater travel time. Highest concentrations appear to be in relatively isolated plumes throughout the domain, with some in the northeast spreading from areas of low concentration and some in the southwest beyond St. Marys River. Nitrogen concentrations range from 0.00156 to 192 mg/L. 65 Figure 2.13. Nitrogen concentrations. Spatial distribution of TN concentrations in groundwater solution for the Quaternary aquifer. White areas within the domain indicate dry cells. Concentrations range from 0.00156 to 192 mg/L. Concentrations generally follow patterns in groundwater travel time, with high concentrations in areas of low groundwater travel time. 4.4 Nutrient Loads Results indicate that a significant proportion of nutrient loads from the Maumee River watershed can be attributed to indirect groundwater discharge to streams and tile drains, while only a very small percentage of loads are attributed to direct groundwater discharge to Lake Erie. Total phosphorus loading to Lake Erie from groundwater ranges from 110.5 to 279.2 kg/day, depending on adsorption capacity. Total nitrogen loading to Lake Erie from groundwater is estimated to be 4,457.1 kg/day. Streams export the largest loads of phosphorus and nitrogen from groundwater within the study area, exporting about 55.0–65.6% and 50.6% of total phosphorus and nitrogen loads from groundwater respectively. Tile drains also have high loads, estimated to be about 33.5–43.7 % of TP from groundwater and 48.0% of TN from groundwater. Highest loads of TP from groundwater to tile drains are simulated at intermediate values of adsorption 66 capacity (10-6 mg/mg), while highest loads to streams are simulated at high and low values of adsorption capacity. Surface ditches and canals had the lowest loads via indirect pathways, contributing about 0.94–1.3% of TP from groundwater and 1.4% of TN from groundwater. Direct loading from groundwater was much lower, contributing less than 0.0001% of TP and TN loads from groundwater in all simulations. Due to groundwater transport delays and desorption of stored nutrients, continued nutrient discharge from groundwater is predicted to contribute significant loads for decades to come, even if all nutrient leaching to groundwater was eliminated today. Average groundwater travel time in the model domain is 181 years, while median groundwater travel time is 110 years. Combined with decades-long legacies of nutrient input to the landscape, released via desorption when recharge concentrations are low relative to sorbed concentrations, delays in groundwater export of nutrients can lead to future loads that are similar or even higher those observed today. Model simulations indicate that after 100 years of near-complete nutrient leaching cessation (except atmospheric deposition), loads from groundwater in the Maumee River watershed may be as high as 110.0–247.5 kilograms per day for phosphorus (99.5% and 88.6% of estimated loads today) and 3,678.6 kilograms per day for nitrogen (82.5% of estimated loads today). Future nutrient delivery with reduced leaching to groundwater shifts towards streams, which receive >60% of groundwater nutrients in all simulations due to the relatively shorter time scales of tile drainage. 4.5 Nutrient Validation Comparison of nutrient modeling results with upstream point loads to concentrations and loads observed in baseflow reveal a similar range of values despite low correlations (Figure 2.14). Simple linear regressions performed on simulated TP concentrations with point loads 67 against observed concentrations have low R2 values, ranging from 0.0749 to 0.0893 (p < 0.01). Lower adsorption capacities have higher R2 values. Linear regressions on TP loads have R2 values ranging from 0.333 to 0.341. For TN, R2 values were higher, at 0.264 for concentrations and 0.606 for loads. However, errors in concentration and yield plotted against observation site watershed size show that the largest errors occur in small watersheds, where local variation is less likely to be captured by the model and has a higher influence on observed values. A B C D Figure 2.14. Nutrient validation. A) Simulated and observed TP loads for an adsorption capacity of 10-6 mg/mg. Simple linear regression shows a significant R2 value of 0.339. Loads show similar magnitudes in simulated and observed values. B) Residual TP concentrations for an adsorption capacity of 10-6 mg/mg shown against the size of each observation site’s watershed. Largest errors are in small watersheds. C) Simulated and observed TN loads. Linear regression shows a significant R2 value of 0.606. Again, simulated and observed loads have a similar magnitude. D) Residual TN concentrations versus the size of each observation point’s watershed. Again, the largest errors tend to be in the smallest watersheds. 68 5. Discussion 5.1 Contribution of Groundwater Pathways to Maumee Nutrient Loads Estimated groundwater nutrient loads make up a notable portion of total nutrient loads from the Maumee River watershed to Lake Erie. Groundwater loads are often lumped in total accounts of nutrient loading with tributary (as baseflow or input to tile drains) and nonpoint source loads. This practice is highly problematic for nutrient management efforts due to the time lags associated with groundwater nutrient transport (Meals et al., 2010; Wang et al., 2013; Martin et al., 2017, 2021; Van Meter et al., 2018, 2021; Vero et al., 2018). Based on reported measures of average TP loading from the Maumee River watershed (Maccoux et al., 2016), our simulated groundwater loads (including both direct and indirect loads via baseflow, tile drains, and ditches) make up 1.5–3.9% of total loads, with the low estimate based on high adsorption capacity and high estimate based on low adsorption capacity. However, given that groundwater is largely unable to transport large amounts of particulate phosphorus, we can reasonably assume that all phosphorus exported from groundwater is soluble reactive phosphorus. If we make this assumption, simulated groundwater loads make up 7.2–18.2% of soluble reactive phosphorus loadings from the Maumee River watershed (Maccoux et al., 2016). Based on simulated SPARROW model TP loads from the Maumee River watershed (Robertson et al., 2019a), our simulated groundwater loads make up 1.7–4.3%. Reported time-averaged estimates of TN loading from the Maumee River watershed were not found in the literature. However, utilizing smoothed TN trends from 2000-2010 (Stow et al., 2015), we estimate recent TN loads to be about 3,150 metric tonnes per month, or 37,800 metric tonnes annually. Of this, our simulated groundwater loads contribute 4.3%. Compared to simulated SPARROW model TN loads (Robertson et al., 2019a), our estimates of groundwater loads contribute 4.0% of loads from the 69 Maumee River watershed, which support assertions that groundwater can be a significant contributor of nutrients to lakes and large tributaries (e.g., Brock et al., 1982; Belanger et al., 1985; Oberdorfer et al., 1990; Valiela et al., 1990; Hagerthey and Kerfoot, 1998; Kang et al., 2005; Holman et al., 2008; Nakayama and Watanabe, 2008; Meinikmann et al., 2013, 2015; Kidmose et al., 2015; Lewandowski et al., 2015; McDowell et al., 2015; Rosenberry et al., 2015; Robinson, 2015; Schilling et al., 2016; Martin et al., 2017; Kazmierczak et al., 2020; Förster et al., 2021; Brookfield et al., 2021). Our results also directly refute findings that sorption and other natural subsurface attenuation processes negate groundwater nutrient transport, especially for phosphorus (Logan et al., 1980; Edwards and Withers, 2007; King et al., 2015). Instead, our simulations are consistent with other evidence that sorption is instead a temporary storage of phosphorus, leading to desorption upon saturation of the sorption capacity or reversal of diffusion gradients (Sharpley et al., 2013; Van Meter et al., 2021). 5.2 Legacy Nutrient Loads While reducing nutrient inputs today is essential to reduce the formation of large harmful algal blooms in the future, we must also acknowledge the legacy impacts of past inputs to advise landscape management and meet water quality goals (Meals et al., 2010; Wang et al., 2013, 2016; Van Meter et al., 2016, 2018, 2021; Martin et al., 2017, 2021; Vero et al., 2018; Liu et al., 2021). The presence of legacy nutrients in the landscape leads to a lag time in land use changes or management practices and the effects those changes have on nutrient concentrations in surface waters such as Lake Erie (Meals et al., 2010; Wang et al., 2013, 2016; Van Meter et al., 2016, 2018, 2021; Martin et al., 2017, 2021; Vero et al., 2018; Liu et al., 2021). Nutrients have two types of legacies: a hydrologic legacy, whose lag time is dependent on the travel times of groundwater carrying the nutrients; and a biogeochemical legacy, whose lag time is dependent 70 on accumulated subsurface nutrient mass and rates of biogeochemical transformation (e.g., sorption/desorption, mineralization, microbial degradation) (Van Meter et al., 2016). The combination of these legacies compound lag times, so it is essential for nutrient reduction efforts to quantify the amount of nutrients present as legacies and the associated potential lag times. Here, we quantify nutrient legacies by modeling nutrient loads 100 years into the future with no inputs during that time, except for input from atmospheric deposition. Simulations show that the vast majority of nutrients present in groundwater today will still be present 100 years in the future. Without concerted remediation efforts on nutrients already in the landscape, loads from groundwater will still be 82–99% of today’s loads from groundwater. Loads will diminish gradually over time, but the timescales on which reductions naturally happen are not feasible for current stakeholders. 5.3 Nitrogen Contamination in Drinking Water The World Health Organization sets the limit for nitrate as nitrogen at 11.3 mg/L due to concerns about methemoglobinemia (2011). Simulated total nitrogen concentrations exceed this limit for 264 km2 of the study area, or roughly 1.4%. Most nitrogen in groundwater occurs as nitrate (Lewandowski et al., 2015), so the potential for high nitrate concentrations is dangerous for those in the study area who rely on groundwater as their source of potable water. Due to the hydrologic legacy, these concentrations are likely to remain high for decades, perpetuating this water quality issue. After 100 years of no nutrient leaching, except for atmospheric deposition, simulations show that concentrations will still exceed nitrogen limits for 0.23% of the domain. Natural attenuation and reduced inputs will not be sufficient to protect future drinking water quality in all groundwater-reliant communities; instead, active groundwater remediation is necessary. 71 5.4 Model Validation Validation of simulated nutrient concentrations and loads showed disparities in simulated versus observed values in baseflow, but overall magnitudes were similar. Variations in exact values should be expected in this case due to the nature of the observations used. For TP, observed soluble reactive phosphorus concentrations and loads may have been more appropriate for comparison to simulations, as groundwater is largely unable to transport significant amounts of particulate matter. Unfortunately, very few observations for soluble reactive phosphorus were available in the study area. For TN and TP, nutrient loading to baseflow can come from many sources outside of groundwater, such as point loading (e.g., wastewater treatment plant effluent), which we addressed based on available data. Other sources, such as runoff from dry season irrigated cropland, overland flow from dry season precipitation, or erosion of riverbed sediments could not be accounted for. We did not model in-stream processing of nutrients, which can also affect baseflow nutrient concentrations. In addition, the largest errors are located in the smallest watersheds, where local fluctuations are more likely to have a large influence on observed concentrations. Simulation resolution may not be able to accurately characterize concentrations in these watersheds. 5.5 Tile Drainage In highly agricultural areas, tile drainage is an especially important consideration for nutrient modeling and management decisions, as it can both increase loads in overland runoff (Blann et al., 2009) and decrease groundwater travel times, short-circuiting natural subsurface attenuation mechanisms (Schilling et al., 2015). Without accurate representation of tile drains in nutrient transport models, contributions from groundwater and other subsurface transport pathways may be vastly underestimated. Modeling results indicate that tile drains in the Maumee 72 River watershed receive a high percentage of groundwater nutrient loads, emphasizing their importance to groundwater management considerations. These results support previous findings that tile drains carry a large proportion of total nutrient losses (Smith et al., 2015). To reduce future loads from groundwater and legacy nutrient stores, stakeholders should consider BMPs that target tile drainage, including drainage water management (Williams et al., 2015) and reduction of drainage density, as well as those that target agriculture. As the largest source of both nitrogen and phosphorus to Lake Erie (Hamlin et al., 2020), agricultural source reduction is essential. Reduced application of fertilizer may be a feasible solution, as previous studies have shown that legacy nutrient stores can sustain crop yields for some years with decreased fertilizer applications (Dalgaard et al., 2014; Rudolph et al., 2015; Muenich et al., 2016). Improved understanding of optimal fertilizer application to individual crops is needed to both increase crop yield and protect future water quality (Cammarano et al., 2021). In addition, manure application to fields is a significant contributor to nutrient leaching (Hamlin et al., 2020; Basso and Ritchie, 2005). Reductions in manure application to fields also has the potential to improve water quality and future nutrient loads. 5.6 Limitations As a process-based model, this study is reliant on an array of outside datasets, each of which have their own limitations and errors. We assumed legacy nutrient inputs to be constant from 1900-1952, with no nutrient input to groundwater before 1900. This time constraint may not allow groundwater nutrients sufficient time to travel to their current positions. Our groundwater flow model was built as a steady-state model, and does not simulate transient changes in flow conditions. With effects of climate change readily becoming apparent, this assumption becomes less valid with time, especially for future projections. Improvements could 73 be made to the model by incorporating spatially explicit estimates of hydraulic conductivity rather than relying on general assumptions about geological deposits. Additionally, more intensive calibration may improve comparisons between modeled results and observation points. Observed groundwater concentrations and loads, in addition to data on soluble reactive phosphorus, may have reduced much uncertainty in model performance. Arguably the most significant limitation was a lack of data on sorption isotherm parameters in aquifers of the area, leading to highly variable estimates of phosphorus loading and concentrations. Developing a fully coupled groundwater-surface water model was outside the scope of this study, but would likely lead to improved simulations and provide additional information about interactions of surface flow routing with tile drains, in-stream nutrient processing, vadose and hyporheic zone dynamics, and the effectiveness of surface-based best management practices. Here, we lay the foundation for a coupled groundwater-surface water model by providing focus and context to the understudied groundwater component of transport. 6. Conclusions In this study, we used MODFLOW and MT3D-USGS to build a spatially explicit, process-based model of groundwater nutrient transport from the Maumee River watershed to Lake Erie. We found that groundwater contributes as much as 4% of total phosphorus and nitrogen loads from the Maumee River watershed, or as much as 18% of soluble reactive phosphorus loads if we assume all phosphorus from groundwater is in a soluble reactive form. Most of these loads are delivered to Lake Erie as baseflow in streams (51–66%) or via tile drains (34–48%). We find that even if all nutrient leaching to groundwater were to stop today, except for atmospheric deposition, loads would be sustained at 82–99% of current loads for at least 100 years. To improve future nutrient loading, we recommend both that surface loading, as the 74 largest nutrient pathway, and groundwater legacies be reduced. Reducing nutrient inputs to the landscape, for example via decreased fertilizer application, is essential for decreased loading from all pathways and reduced future accumulation of subsurface legacies. However, groundwater legacies have the potential to sustain loads in the long-term when other sources are lowered. We argue that reducing the legacy store of subsurface nutrients, for example via tile drainage management or groundwater remediation, is also an important step to improve water quality. Simulations indicate that 264 km2 of the domain likely has groundwater nitrogen concentrations exceeding limits set by the World Health Organization. This is a significant concern for current and future groundwater-reliant communities, as concentrations are expected to remain dangerously high for several decades in some areas given estimated groundwater travel times. Active groundwater remediation is necessary to protect these communities. Future studies should couple spatially explicit groundwater models, like the one built here, and surface water models for a more holistic picture of nutrient transport and BMPs. Models with more accurate sorption isotherm parameters would constrain simulated phosphorus concentrations. 7. Acknowledgments I would like to thank Quercus Hamlin, Luwen Wan, Sherry Martin, Anthony Kendall, Dave Hyndman, and others on the SENSEmap/SENSEflux team for providing baseline data for this chapter. Funding was provided by Michigan State University, USDA-National Institute of Food and Agriculture/National Science Foundation award number 2018-67003-27406, and NOAA grant NA12OAR4320071. 75 REFERENCES Apostel, A., Kalcic, M., Dagnew, A., Evenson, G., Kast, J., King, K., Martin, J., Muenich, R.L., and Scavia, D., 2021, Simulating internal watershed processes using multiple SWAT models: Science of The Total Environment, v. 759, article 143920, doi:10.1016/j.scitotenv.2020.143920. Arhonditsis, G.B., Neumann, A., Shimoda, Y., Kim, D.-K., Dong, F., Onandia, G., Yang, C., Javed, A., Brady, M., Visha, A., Ni, F., and Cheng, V., 2019, Castles built on sand or predictive limnology in action? Part A: Evaluation of an integrated modelling framework to guide adaptive management implementation in Lake Erie: Ecological Informatics, v. 53, article 100969, doi:10.1016/j.ecoinf.2019.05.015. Baker, D.B., Confesor, R., Ewing, D.E., Johnson, L.T., Kramer, J.W., and Merryfield, B.J., 2014, Phosphorus loading to Lake Erie from the Maumee, Sandusky and Cuyahoga rivers: The importance of bioavailability: Journal of Great Lakes Research, v. 40, p. 502–517, doi:10.1016/j.jglr.2014.05.001. Bakker, M. et al., 2021, FloPy v3.3.3 — release candidate: U.S. Geological Survey Software Release. Basso, B., and Ritchie, J.T., 2005, Impact of compost, manure and inorganic fertilizer on nitrate leaching and yield for a 6-year maize–alfalfa rotation in Michigan: Agriculture, Ecosystems & Environment, v. 108, p. 329–341, doi:10.1016/j.agee.2005.01.011. Bedekar, V., Morway, E.D., Langevin, C.D., and Tonkin, M., 2016, MT3D-USGS version 1.0.0: Groundwater Solute Transport Simulator for MODFLOW: U.S. Geological Survey Software Release, 30 September 2016, doi: 10.5066/F75T3HKD. Belanger, T.V., Mikutel, D.F., and Churchill, P.A., 1985, Groundwater seepage nutrient loading in a Florida Lake: Water Research, v. 19, p. 773–781, doi:10.1016/0043-1354(85)90126-5. Blann, K.L., Anderson, J.L., Sands, G.R., and Vondracek, B., 2009, Effects of agricultural drainage on aquatic ecosystems: A review: Critical Reviews in Environmental Science and Technology, v. 39, p. 909–1001, doi:10.1080/10643380801977966. Bosch, N.S., Allan, J.D., Dolan, D.M., Han, H., and Richards, R.P., 2011, Application of the Soil and Water Assessment Tool for six watersheds of Lake Erie: Model parameterization and calibration: Journal of Great Lakes Research, v. 37, p. 263–271, doi:10.1016/j.jglr.2011.03.004. Botts, L., and Krushelnicki, B., 1995, The Great Lakes: An Environmental Atlas and Resource Book (3rd ed.) (K. Fuller, H. Shear, & H. Wittig, Eds.): Toronto, Canada and Chicago, Illinois, Government of Canada and United States Environmental Protection Agency, 46 p. Brock, T.D., Lee, D.R., Janes, D., and Winek, D., 1982, Groundwater seepage as a nutrient source to a drainage lake; Lake Mendota, Wisconsin: Water Research, v. 16, p. 1255–1263, doi:10.1016/0043-1354(82)90144-0. 76 Brookfield, A.E., Hansen, A.T., Sullivan, P.L., Czuba, J.A., Kirk, M.F., Li, L., Newcomer, M.E., and Wilkinson, G., 2021, Predicting algal blooms: Are we overlooking groundwater? Science of The Total Environment, v. 769, 15 p., doi:10.1016/j.scitotenv.2020.144442. Brown, J.F., Tollerud, H.J., Barber, C.P., Zhou, Q., Dwyer, J.L., Vogelmann, J.E., Loveland, T.R., Woodcock, C.E., Stehman, S.V., Zhu, Z., Pengra, B.W., Smith, K., Horton, J.A., Xian, G., Auch, R.F., Sohl, T.L., Sayler, K.L., Gallant, A.L., Zelenak, D., Reker, R.R., and Rover, J., 2020, Lessons learned implementing an operational continuous United States national land change monitoring capability—The Land Change Monitoring, Assessment, and Projection (LCMAP) approach: Remote Sensing of Environment, v. 238, article 111356, doi: 10.1016/j.rse.2019.111356. Calhoun, F.G., Baker, D.B., and Slater, B.K., 2002, Soils, water quality, and watershed size: Interactions in the Maumee and Sandusky River basins of northwestern Ohio: Journal of Environmental Quality, v. 31, p. 47–53, doi:10.2134/jeq2002.4700a. Cammarano, D., Basso, B., Holland, J., Gianinetti, A., Baronchelli, M., and Ronga, D., 2021, Modeling spatial and temporal optimal N fertilizer rates to reduce nitrate leaching while improving grain yield and quality in malting barley: Computers and Electronics in Agriculture, v. 182, article 105997, doi:10.1016/j.compag.2021.105997. Chaffin, J.D., Bridgeman, T.B., Bade, D.L., and Mobilian, C.N., 2014, Summer phytoplankton nutrient limitation in Maumee Bay of Lake Erie during high-flow and low-flow years: Journal of Great Lakes Research, v. 40, p. 524–531, doi:10.1016/j.jglr.2014.04.009. Chapuis, R.P., and Gill, D.E., 1989, Hydraulic anisotropy of homogeneous soils and rocks: Influence of the densification process: Bulletin of the International Association of Engineering Geology, v. 39, p. 75–86. Cheng, X., Zeng, Y., Guo, Z., and Zhu, L., 2014, Diffusion of nitrogen and phosphorus across the sediment-water interface and in seawater at aquaculture areas of Daya Bay, China: International Journal of Environmental Research and Public Health, v. 11, p. 1557–1572, doi:10.3390/ijerph110201557. Choquette, A.F., Hirsch, R.M., Murphy, J.C., Johnson, L.T., and Confesor, R.B., 2019, Tracking changes in nutrient delivery to western Lake Erie: Approaches to compensate for variability and trends in streamflow: Journal of Great Lakes Research, v. 45, p. 21–39, doi:10.1016/j.jglr.2018.11.012. Cousino, L.K., Becker, R.H., and Zmijewski, K.A., 2015, Modeling the effects of climate change on water, sediment, and nutrient yields from the Maumee River watershed: Journal of Hydrology: Regional Studies, v. 4, p. 762–775, doi:10.1016/j.ejrh.2015.06.017. Cronan, C.S., and Aiken, G.R., 1985, Chemistry and transport of soluble humic substances in forested watersheds of the Adirondack Park, New York: Geochimica et Cosmochimica Acta, v. 49, p. 1697–1705, doi:10.1016/0016-7037(85)90140-1. Culbertson, A.M., Martin, J.F., Aloysius, N., and Ludsin, S.A., 2016, Anticipated impacts of 77 climate change on 21st century Maumee River discharge and nutrient loads: Journal of Great Lakes Research, v. 42, p. 1332–1342, doi:10.1016/j.jglr.2016.08.008. Cussler, E.L., 1997, Mass Transfer in Fluid Systems (2nd ed.): Cambridge, United Kingdom, Cambridge University Press, 580 p. Dagnew, A., Scavia, D., Wang, Y.-C., Muenich, R., and Kalcic, M., 2019, Modeling phosphorus reduction strategies from the international St. Clair-Detroit River system watershed: Journal of Great Lakes Research, v. 45, p. 742–751, doi:10.1016/j.jglr.2019.04.005. Dalgaard, T., Hansen, B., Hasler, B., Hertel, O., Hutchings, N.J., Jacobsen, B.H., Jensen, L.S., Kronvang, B., Olesen, J.E., Schjørring, J.K., Kristensen, I.S., Graversgaard, M., Termansen, M., and Vejre, H., 2014, Policies for agricultural nitrogen management—trends, challenges and prospects for improved efficiency in Denmark: Environmental Research Letters, v. 9, article 115002, doi:10.1088/1748-9326/9/11/115002. Daloğlu, I., Cho, K.H., and Scavia, D., 2012, Evaluating causes of trends in long-term dissolved reactive phosphorus loads to Lake Erie: Environmental Science & Technology, v. 46, p. 10660–10666, doi:10.1021/es302315d. Dewitz, J., 2019, National Land Cover Database (NLCD) 2019 Land Cover Conterminous United States: United States Geological Survey. Dobiesz, N.E., Hecky, R.E., Johnson, T.B., Sarvala, J., Dettmers, J.M., Lehtiniemi, M., Rudstam, L.G., Madenjian, C.P., and Witte, F., 2010, Metrics of ecosystem status for large aquatic systems – A global comparison: Journal of Great Lakes Research, v. 36, p. 123–138, doi:10.1016/j.jglr.2009.11.003. Domenico, P.A., and Schwartz, F.W., 1990, Physical and Chemical Hydrogeology: New York, John Wiley and Sons, 824 p. Dubrovsky, N.M., and Hamilton, P.A., 2010, Nutrients in the nation’s streams and groundwater: National findings and implications: U.S. Geological Survey, v. USGS Fact Sheet 2010– 3078, p. 6. Edwards, A.C., and Withers, P.J.A., 2007, Linking phosphorus sources to impacts in different types of water body: Soil Use and Management, v. 23, p. 133–143, doi:10.1111/j.1475- 2743.2007.00110.x. Falcone, J.A., 2021, Tabular county-level nitrogen and phosphorus estimates from fertilizer and manure for approximately 5-year periods from 1950 to 2017: U.S. Geological Survey data release, https://doi.org/10.5066/P9VSQN3C. Fang, S., Del Giudice, D., Scavia, D., Binding, C.E., Bridgeman, T.B., Chaffin, J.D., Evans, M.A., Guinness, J., Johengen, T.H., and Obenour, D.R., 2019, A space-time geostatistical model for probabilistic estimation of harmful algal bloom biomass and areal extent: Science of the Total Environment, v. 695, article 133776, doi:10.1016/j.scitotenv.2019.133776. 78 Förster, W., Scholten, J.C., Schubert, M., Knoeller, K., Classen, N., Lechelt, M., Richard, J.-H., Rohweder, U., Zunker, I., and Wanner, S.C., 2021, Phosphorous supply to a eutrophic artificial lake: Sedimentary versus groundwater sources: Water, v. 13, 20 p., doi:10.3390/w13040563. Forster, D.L., Smith, E.C., and Hite, D., 2000, A bioeconomic model of farm management practices and environmental effluents in the western Lake Erie Basin: Journal of Soil and Water Conservation, v. 55, p. 177–182. Freeze, R.A., and Cherry, J.A., 1979, Groundwater: Englewood Cliffs, New Jersey, Prentice- Hall, 604 p. Geologic Survey Division of the Michigan Department of Environmental Quality (MDEQ), 1987, Bedrock Geology: Michigan Department of Environmental Quality. Governments of Canada and the U.S.A., 2012, Great Lakes Water Quality Agreement: International Joint Commission, p. 56. Gray, H.H., Ault, C., Keller, S., and Harper, D., 2010, BEDROCK_GEOL_MM48_IN: Bedrock Geology of Indiana (Indiana Geological Survey, 1:500,000, Polygon Shapefile): Indiana Geological Survey. Gray, H.H., and Walls, C., 2002, SURFICIAL_GEOL_MM49_IN: Quaternary Geologic Map of Indiana (Indiana Geological Survey, 1:500,000, Polygon Shapefile): Indiana Geological Survey. Gunn, K.M., Baule, W.J., Frankenberger, J.R., Gamble, D.L., Allred, B.J., Andresen, J.A., and Brown, L.C., 2018, Modeled climate change impacts on subirrigated maize relative yield in northwest Ohio: Agricultural Water Management, v. 206, p. 56–66, doi:10.1016/j.agwat.2018.04.034. Guo, T., Confesor, R., Saleh, A., and King, K., 2020, Crop growth, hydrology, and water quality dynamics in agricultural fields across the Western Lake Erie Basin: Multi-site verification of the Nutrient Tracking Tool (NTT): Science of The Total Environment, v. 726, article 138485, doi:10.1016/j.scitotenv.2020.138485. Hagerthey, S.E., and Kerfoot, W.C., 1998, Groundwater flow influences the biomass and nutrient ratios of epibenthic algae in a north temperate seepage lake: Limnology and Oceanography, v. 43, p. 1227–1242, doi:10.4319/lo.1998.43.6.1227. Haines, M., Fishback, P., and Rhode, P., 2018, United States Agriculture Data, 1840 – 2012: Inter-university Consortium for Political and Social Research [distributor], https://doi.org/10.3886/ICPSR35206.v4. Hamlin, Q.F., Kendall, A.D., Martin, S.L., Whitenack, H.D., Roush, J.A., Hannah, B.A., and Hyndman, D.W., 2020, Quantifying landscape nutrient inputs with spatially explicit nutrient source estimate maps: Journal of Geophysical Research: Biogeosciences, v. 125, article e2019JG005134, doi:10.1029/2019JG005134. 79 Hannah, B.A., Kendall, A.D., Martin, S.L., and Hyndman, D.W., 2020, Quantifying linkages between watershed factors and coastal wetland plant invasion in the US Great Lakes: Landscape Ecology, v. 35, p. 2843–2861, doi:10.1007/s10980-020-01124-3. Harbaugh, A.W., Langevin, C.D., Hughes, J.D., Niswonger, R.N., and Konikow, L.F., 2017, MODFLOW-2005 version 1.12.00, the U.S. Geological Survey modular groundwater model: U.S. Geological Survey Software Release, http://dx.doi.org/10.5066/F7RF5S7G. Hartig, J.H., and Wallen, D.G., 1984, Seasonal variation of nutrient limitation in western Lake Erie: Journal of Great Lakes Research, v. 10, p. 449–460, doi:10.1016/S0380- 1330(84)71862-4. Hegglin, M. I., D. Kinnison, D. Plummer, et al., Historical and future ozone database (1850- 2100) in support of CMIP6, GMD, in preparation. Holman, I.P., Whelan, M.J., Howden, N.J.K., Bellamy, P.H., Willby, N.J., Rivas-Casado, M., and McConvey, P., 2008, Phosphorus in groundwater—an overlooked contributor to eutrophication? Hydrological Processes, v. 22, p. 5121–5127, doi:10.1002/hyp.7198. Hu, Y., Long, C.M., Wang, Y.-C., Kerkez, B., and Scavia, D., 2019, Urban total phosphorus loads to the St. Clair-Detroit River System: Journal of Great Lakes Research, v. 45, p. 1142–1149, doi:10.1016/j.jglr.2019.09.009. Hyndman, D.W., Kendall, A.D., and Welty, N.R.H., 2007, Evaluating temporal and spatial variations in recharge and streamflow using the Integrated Landscape Hydrology Model (ILHM), in Hyndman, D.W., Day-Lewis, F.D., and Singha, K. eds., Geophysical Monograph Series, Washington, D. C., American Geophysical Union, v. 171, p. 121–141, doi:10.1029/171GM11. Igras, J.D., and Creed, I.F., 2020, Uncertainty analysis of the performance of a management system for achieving phosphorus load reduction to surface waters: Journal of Environmental Management, v. 276, article 111217, doi:10.1016/j.jenvman.2020.111217. Jankowiak, J., Hattenrath-Lehmann, T., Kramer, B.J., Ladds, M., and Gobler, C.J., 2019, Deciphering the effects of nitrogen, phosphorus, and temperature on cyanobacterial bloom intensification, diversity, and toxicity in western Lake Erie: Limnology and Oceanography, v. 64, p. 1347–1370, doi:10.1002/lno.11120. Jarvie, H.P., Johnson, L.T., Sharpley, A.N., Smith, D.R., Baker, D.B., Bruulsema, T.W., and Confesor, R., 2017, Increased soluble phosphorus loads to Lake Erie: Unintended consequences of conservation practices? Journal of Environmental Quality, v. 46, p. 123– 132, doi:10.2134/jeq2016.07.0248. Jetoo, S., Grover, V.I., and Krantzberg, G., 2015, The Toledo drinking water advisory: Suggested application of the water safety planning approach: Sustainability, v. 7, p. 9787– 9808, doi:10.3390/su7089787. Kaatz, M.R., 1955, The Black Swamp: A study in historical geography: Annals of the 80 Association of American Geographers, v. 45, p. 35, doi:10.1111/j.1467- 8306.1955.tb01481.x. Kane, D.D., Conroy, J.D., Peter Richards, R., Baker, D.B., and Culver, D.A., 2014, Re- eutrophication of Lake Erie: Correlations between tributary nutrient loads and phytoplankton biomass: Journal of Great Lakes Research, v. 40, p. 496–501, doi:10.1016/j.jglr.2014.04.004. Kang, W.-J., Kolasa, K.V., and Rials, M.W., 2005, Groundwater inflow and associated transport of phosphorus to a hypereutrophic lake: Environmental Geology, v. 47, p. 565–575, doi:10.1007/s00254-004-1180-3. Kazmierczak, J., Postma, D., Müller, S., Jessen, S., Nilsson, B., Czekaj, J., and Engesgaard, P., 2020, Groundwater-controlled phosphorus release and transport from sandy aquifer into lake: Limnology and Oceanography, v. 65, p. 2188–2204, doi:10.1002/lno.11447. Kidmose, J., Engesgaard, P., Ommen, D.A.O., Nilsson, B., Flindt, M.R., and Andersen, F.Ø., 2015, The role of groundwater for lake-water quality and quantification of N seepage: Groundwater, v. 53, p. 709–721, doi:10.1111/gwat.12281. King, K.W., Williams, M.R., Johnson, L.T., Smith, D.R., LaBarge, G.A., and Fausey, N.R., 2017, Phosphorus availability in Western Lake Erie Basin drainage waters: Legacy evidence across spatial scales: Journal of Environmental Quality, v. 46, p. 466–469, doi:10.2134/jeq2016.11.0434. King, K.W., Williams, M.R., Macrae, M.L., Fausey, N.R., Frankenberger, J., Smith, D.R., Kleinman, P.J.A., and Brown, L.C., 2015, Phosphorus transport in agricultural subsurface drainage: A review: Journal of Environmental Quality, v. 44, p. 467–485, doi:10.2134/jeq2014.04.0163. Knights, D., Parks, K.C., Sawyer, A.H., David, C.H., Browning, T.N., Danner, K.M., and Wallace, C.D., 2017, Direct groundwater discharge and vulnerability to hidden nutrient loads along the Great Lakes coast of the United States: Journal of Hydrology, v. 554, p. 331–341, doi:10.1016/j.jhydrol.2017.09.001. Kroeger, K.D., Cole, M.L., and Valiela, I., 2006, Groundwater-transported dissolved organic nitrogen exports from coastal watersheds: Limnology and Oceanography, v. 51, p. 2248– 2261, doi:10.4319/lo.2006.51.5.2248. Kujawa, H., Kalcic, M., Martin, J., Aloysius, N., Apostel, A., Kast, J., Murumkar, A., Evenson, G., Becker, R., Boles, C., Confesor, R., Dagnew, A., Guo, T., Muenich, R.L., Redder, T., Scavia, D., and Wang, Y.-C., 2020, The hydrologic model as a source of nutrient loading uncertainty in a future climate: Science of The Total Environment, v. 724, article 138004, doi:10.1016/j.scitotenv.2020.138004. Lewandowski, J., Meinikmann, K., Nützmann, G., and Rosenberry, D.O., 2015, Groundwater – the disregarded component in lake water and nutrient budgets. Part 2: effects of groundwater on nutrients: Hydrological Processes, v. 29, p. 2922–2955, 81 doi:10.1002/hyp.10384. Liu, Y., Engel, B.A., Collingsworth, P.D., and Pijanowski, B.C., 2017, Optimal implementation of green infrastructure practices to minimize influences of land use change and climate change on hydrology and water quality: Case study in Spy Run Creek watershed, Indiana: Science of The Total Environment, v. 601–602, p. 1400–1411, doi:10.1016/j.scitotenv.2017.06.015. Liu, J., Meter, K.J.V., McLeod, M.M., and Basu, N.B., 2021, Checkered landscapes: hydrologic and biogeochemical nitrogen legacies along the river continuum: Environmental Research Letters, v. 16, article 115006, doi:10.1088/1748-9326/ac243c. Liu, Y., Wang, R., Guo, T., Engel, B.A., Flanagan, D.C., Lee, J.G., Li, S., Pijanowski, B.C., Collingsworth, P.D., and Wallace, C.W., 2019, Evaluating efficiencies and cost- effectiveness of best management practices in improving agricultural water quality using integrated SWAT and cost evaluation tool: Journal of Hydrology, v. 577, 16 p., doi:10.1016/j.jhydrol.2019.123965. Logan, T.J., 1987, Diffuse (non-point) source loading of chemicals to Lake Erie: Journal of Great Lakes Research, v. 13, p. 649–658, doi:10.1016/S0380-1330(87)71679-7. Logan, T.J., Randall, G.W., and Timmons, D.R., 1980, Nutrient content of tile drainage from cropland in the north central region: Ohio Agricultural Research and Development Center Research Bulletin 1119, North Central Regional Research Publication 268, 19 p. Luscz, E.C., Kendall, A.D., and Hyndman, D.W., 2017, A spatially explicit statistical model to quantify nutrient sources, pathways, and delivery at the regional scale: Biogeochemistry, v. 133, p. 37–57, doi:10.1007/s10533-017-0305-1. Luscz, E.C., Kendall, A.D., and Hyndman, D.W., 2015, High resolution spatially explicit nutrient source models for the Lower Peninsula of Michigan: Journal of Great Lakes Research, v. 41, p. 618–629, doi:10.1016/j.jglr.2015.02.004. Maccoux, M.J., Dove, A., Backus, S.M., and Dolan, D.M., 2016, Total and soluble reactive phosphorus loadings to Lake Erie: A detailed accounting by year, basin, country, and tributary: Journal of Great Lakes Research, v. 42, p. 1151–1165, doi:10.1016/j.jglr.2016.08.005. Manson, S., Schroeder, J., Van Riper, D., Kugler, T., and Ruggles, S., 2021, IPUMS National Historical Geographic Information System: Version 16.0 [dataset]: Minneapolis, MN, IPUMS, http://doi.org/10.18128/D050.V16.0. Martin, S.L., Hamlin, Q.F., Kendall, A.D., Wan, L., and Hyndman, D.W., 2021, The land use legacy effect: looking back to see a path forward to improve management: Environmental Research Letters, v. 16, 11 p., doi:10.1088/1748-9326/abe14c. Martin, S.L., Hayes, D.B., Kendall, A.D., and Hyndman, D.W., 2017, The land-use legacy effect: Towards a mechanistic understanding of time-lagged water quality responses to land 82 use/cover: Science of the Total Environment, v. 579, p. 1794–1803, doi:10.1016/j.scitotenv.2016.11.158. McCallister, D.L., and Logan, T.J., 1978, Phosphate adsorption-desorption characteristics of soils and bottom sediments in the Maumee River Basin of Ohio: Journal of Environmental Quality, v. 7, p. 87–92, doi:10.2134/jeq1978.00472425000700010018x. McDowell, R.W., Cox, N., Daughney, C.J., Wheeler, D., and Moreau, M., 2015, A national assessment of the potential linkage between soil, and surface and groundwater concentrations of phosphorus: JAWRA Journal of the American Water Resources Association, v. 51, p. 992–1002, doi:10.1111/1752-1688.12337. Meals, D.W., Dressing, S.A., and Davenport, T.E., 2010, Lag time in water quality response to best management practices: A review: Journal of Environmental Quality, v. 39, p. 85–96, doi:10.2134/jeq2009.0108. Mehan, S., Aggarwal, R., Gitau, M.W., Flanagan, D.C., Wallace, C.W., and Frankenberger, J.R., 2019, Assessment of hydrology and nutrient losses in a changing climate in a subsurface- drained watershed: Science of The Total Environment, v. 688, p. 1236–1251, doi:10.1016/j.scitotenv.2019.06.314. Meinikmann, K., Hupfer, M., and Lewandowski, J., 2015, Phosphorus in groundwater discharge – A potential source for lake eutrophication: Journal of Hydrology, v. 524, p. 214–226, doi:10.1016/j.jhydrol.2015.02.031. Meinikmann, K., Lewandowski, J., and Nützmann, G., 2013, Lacustrine groundwater discharge: Combined determination of volumes and spatial patterns: Journal of Hydrology, v. 502, p. 202–211, doi:10.1016/j.jhydrol.2013.08.021. Merriman, K.R., Daggupati, P., Srinivasan, R., Toussant, C., Russell, A.M., and Hayhurst, B., 2018, Assessing the impact of site-specific BMPs using a spatially explicit, field-scale SWAT model with edge-of-field and tile hydrology and water-quality data in the Eagle Creek Watershed, Ohio: Water, v. 10, article 1299, doi:10.3390/w10101299. Michigan Natural Features Inventory, 1999, Quaternary Geology of Michigan: Michigan Department of Natural Resources. Mitsch, W.J., 2017, Solving Lake Erie’s harmful algal blooms by restoring the Great Black Swamp in Ohio: Ecological Engineering, v. 108, p. 406–413, doi:10.1016/j.ecoleng.2017.08.040. Molina-Navarro, E., Bailey, R.T., Andersen, H.E., Thodsen, H., Nielsen, A., Park, S., Jensen, J.S., Jensen, J.B., and Trolle, D., 2019, Comparison of abstraction scenarios simulated by SWAT and SWAT-MODFLOW: Hydrological Sciences Journal, v. 64, p. 434–454, doi:10.1080/02626667.2019.1590583. Mooney, R.J., Stanley, E.H., Rosenthal, W.C., Esselman, P.C., Kendall, A.D., and McIntyre, P.B., 2020, Outsized nutrient contributions from small tributaries to a Great Lake: 83 Proceedings of the National Academy of Sciences, v. 117, p. 28175–28182, doi:10.1073/pnas.2001376117. Muenich, R.L., Kalcic, M., and Scavia, D., 2016, Evaluating the impact of legacy P and agricultural conservation practices on nutrient loads from the Maumee River Watershed: Environmental Science & Technology, v. 50, p. 8146–8154, doi:10.1021/acs.est.6b01421. Murray, C.J., Müller-Karulis, B., Carstensen, J., Conley, D.J., Gustafsson, B.G., and Andersen, J.H., 2019, Past, present and future eutrophication status of the Baltic Sea: Frontiers in Marine Science, v. 6, 12 p. Nakagaki, N. and Wieczorek, M.E., 2016, Estimates of subsurface tile drainage extent for 12 Midwest states, 2012: U.S. Geological Survey data release, doi: 10.5066/F7W37TDP. Nakayama, T., and Watanabe, M., 2008, Missing role of groundwater in water and nutrient cycles in the shallow eutrophic Lake Kasumigaura, Japan: Hydrological Processes, v. 22, p. 1150–1172, doi:10.1002/hyp.6684. National Agricultural Statistics Service (NASS), 2021, Cropland Data Layer: U.S. Department of Agriculture (USDA). National Oceanic and Atmospheric Administration (NOAA), Tides and Currents: Water Level Stations: https://tidesandcurrents.noaa.gov/stations.html?type=Water+Levels (accessed August 2021). Neitsch, S.L., Arnold, J.G., Kiniry, J.R., and Williams, J.R., 2011, Soil and Water Assessment Tool Theoretical Documentation: Version 2009: Texas Water Resources Institute, report no. 406, 618 p. Nelson, A.M., Moriasi, D.N., Talebizadeh, M., Steiner, J.L., Gowda, P.H., Starks, P.J., and Tadesse, H.K., 2018, Use of soft data for multicriteria calibration and validation of Agricultural Policy Environmental eXtender: Impact on model simulations: Journal of Soil and Water Conservation, v. 73, p. 623–636, doi:10.2489/jswc.73.6.623. Newell, S.E., Davis, T.W., Johengen, T.H., Gossiaux, D., Burtner, A., Palladino, D., and McCarthy, M.J., 2019, Reduced forms of nitrogen are a driver of non-nitrogen-fixing harmful cyanobacterial blooms and toxicity in Lake Erie: Harmful Algae, v. 81, p. 86–93, doi:10.1016/j.hal.2018.11.003. Niswonger, R.G., Panday, S., and Ibaraki, M., 2011, MODFLOW-NWT, A Newton formulation for MODFLOW-2005: U.S. Geological Survey Techniques and Methods 6-A37, 44 p., https://doi.org/10.3133/tm6A37. Nolan, B.T., and Hitt, K.J., 2006, Vulnerability of shallow groundwater and drinking-water wells to nitrate in the United States: Environmental Science & Technology, v. 40, p. 7834–7840, doi:10.1021/es060911u. Oberdorfer, J.A., Valentino, M.A., and Smith, S.V., 1990, Groundwater contribution to the 84 nutrient budget of Tomales Bay, California: Biogeochemistry, v. 10, p. 199–216, doi:10.1007/BF00003144. Oelsner, G.P., and Stets, E.G., 2019, Recent trends in nutrient and sediment loading to coastal areas of the conterminous U.S.: Insights and global context: Science of The Total Environment, v. 654, p. 1225–1240, doi:10.1016/j.scitotenv.2018.10.437. Ohio Environmental Protection Agency (Ohio EPA), 2014, Major aquifers in Ohio and associated water quality: Technical Series on Ground Water Quality, 28 p. Oldfield, L., Rakhimbekova, S., Roy, J.W., and Robinson, C.E., 2020, Estimation of phosphorus loads from septic systems to tributaries in the Canadian Lake Erie Basin: Journal of Great Lakes Research, v. 46, p. 1559–1569, doi:10.1016/j.jglr.2020.08.021. Pabich, W.J., Valiela, I., and Hemond, H.F., 2001, Relationship between DOC concentration and vadose zone thickness and depth below water table in groundwater of Cape Cod, U.S.A.: Biogeochemistry, v. 55, p. 247–268, doi:10.1023/A:1011842918260. Paerl, H.W., and Otten, T.G., 2013, Harmful cyanobacterial blooms: Causes, consequences, and controls: Microbial Ecology, v. 65, p. 995–1010, doi:10.1007/s00248-012-0159-y. Pavey, R. R., Goldthwait, R. P., Brockman, C. S., Hull, D. N., Swinford, E. M., and Van Horn, R. G., 1999, Quaternary geology of Ohio: Ohio Division of Geological Survey Map M-2, 1:500,000-scale map and 1:250,000-scale GIS files. Pollock, D.W., 2017, MODPATH v7.2.01: A particle-tracking model for MODFLOW: U.S. Geological Survey Software Release. Rakhimbekova, S., O’Carroll, D.M., Oldfield, L.E., Ptacek, C.J., and Robinson, C.E., 2021, Spatiotemporal controls on septic system derived nutrients in a nearshore aquifer and their discharge to a large lake: Science of The Total Environment, v. 752, article 141262, doi:10.1016/j.scitotenv.2020.141262. Read, E.K., Carr, L., De Cicco, L., Dugan, H.A., Hanson, P.C., Hart, J.A., Kreft, J., Read, J.S., and Winslow, L.A., 2017, Water quality data for national‐scale aquatic research: The Water Quality Portal: Water Resources Research, v. 53, p. 1735-1745. River, M., and Richardson, C.J., 2019, Dissolved reactive phosphorus loads to western Lake Erie: The hidden influence of nanoparticles: Journal of Environmental Quality, v. 48, p. 645–653, doi:10.2134/jeq2018.05.0178. Robertson, D.M., and Saad, D.A., 2011, Nutrient inputs to the Laurentian Great Lakes by source and watershed estimated using SPARROW watershed models: JAWRA Journal of the American Water Resources Association, v. 47, p. 1011–1033, doi:10.1111/j.1752- 1688.2011.00574.x. Robertson, D.M., Saad, D.A., Benoy, G.A., Vouk, I., Schwarz, G.E., and Laitta, M.T., 2019a, Phosphorus and nitrogen transport in the binational Great Lakes Basin estimated using 85 SPARROW watershed models: JAWRA Journal of the American Water Resources Association, v. 55, p. 1401–1424, doi:10.1111/1752-1688.12792. Robertson, W.D., Van Stempvoort, D.R., and Schiff, S.L., 2019b, Review of phosphorus attenuation in groundwater plumes from 24 septic systems: Science of The Total Environment, v. 692, p. 640–652, doi:10.1016/j.scitotenv.2019.07.198. Robinson, C., 2015, Review on groundwater as a source of nutrients to the Great Lakes and their tributaries: Journal of Great Lakes Research, v. 41, p. 941–950, doi:10.1016/j.jglr.2015.08.001. Rosenberry, D.O., Lewandowski, J., Meinikmann, K., and Nützmann, G., 2015, Groundwater - the disregarded component in lake water and nutrient budgets. Part 1: effects of groundwater on hydrology: Hydrological Processes, v. 29, p. 2895–2921, doi:10.1002/hyp.10403. Rowland, F.E., Stow, C.A., Johengen, T.H., Burtner, A.M., Palladino, D., Gossiaux, D.C., Davis, T.W., Johnson, L.T., and Ruberg, S., 2020, Recent patterns in Lake Erie phosphorus and chlorophyll a concentrations in response to changing loads: Environmental Science & Technology, v. 54, p. 835–841, doi:10.1021/acs.est.9b05326. Rudolph, D. l., Devlin, J. f., and Bekeris, L., 2015, Challenges and a strategy for agricultural BMP monitoring and remediation of nitrate contamination in unconsolidated aquifers: Groundwater Monitoring & Remediation, v. 35, p. 97–109, doi:10.1111/gwmr.12103. Sayers, M.J., Grimm, A.G., Shuchman, R.A., Bosse, K.R., Fahnenstiel, G.L., Ruberg, S.A., and Leshkevich, G.A., 2019, Satellite monitoring of harmful algal blooms in the Western Basin of Lake Erie: A 20-year time-series: Journal of Great Lakes Research, v. 45, p. 508–521, doi:10.1016/j.jglr.2019.01.005. Schilling, K.E., Streeter, M.T., Quade, D., and Skopec, M., 2016, Groundwater loading of nitrate-nitrogen and phosphorus from watershed source areas to an Iowa Great Lake: Journal of Great Lakes Research, v. 42, p. 588–598, doi:10.1016/j.jglr.2016.03.015. Schilling, K.E., Wolter, C.F., Isenhart, T.M., and Schultz, R.C., 2015, Tile drainage density reduces groundwater travel times and compromises riparian buffer effectiveness: Journal of Environmental Quality, v. 44, p. 1754–1763, doi:10.2134/jeq2015.02.0105. Shao, G., Zhang, D., Guan, Y., Xie, Y., and Huang, F., 2019, Application of SWAT model with a modified groundwater module to the semi-arid Hailiutu River Catchment, Northwest China: Sustainability, v. 11, 20 p., doi:10.3390/su11072031. Sharma, S., Nalley, D., and Subedi, N., 2018, Characterization of temporal and spatial variability of phosphorus loading to Lake Erie from the western basin using wavelet transform methods: Hydrology, v. 5, 50 p., doi:10.3390/hydrology5030050. Sharpley, A., Jarvie, H.P., Buda, A., May, L., Spears, B., and Kleinman, P., 2013, Phosphorus legacy: Overcoming the effects of past management practices to mitigate future water 86 quality impairment: Journal of Environmental Quality, v. 42, p. 1308–1326, doi:10.2134/jeq2013.03.0098. Sinsabaugh, R.L., Findlay, S., Franchini, P., and Fischer, D., 1997, Enzymatic analysis of riverine bacterioplankton production: Limnology and Oceanography, v. 42, p. 29–38, doi:10.4319/lo.1997.42.1.0029. Sinsabaugh, R.L., and Shah, J.J.F., 2010, Integrating resource utilization and temperature in metabolic scaling of riverine bacterial production: Ecology, v. 91, p. 1455–1465, doi:10.1890/08-2192.1. Slucher, E.R., (principal compiler), Swinford, E.M., Larsen, G.E., and others, with GIS production and cartography by Powers, D.M., 2006, Bedrock geologic map of Ohio: Ohio Division of Geological Survey Digital Map Series BG-1, version 6.0, 1:500,000 scale. Smith, D.R., King, K.W., Johnson, L., Francesconi, W., Richards, P., Baker, D., and Sharpley, A.N., 2015, Surface runoff and tile drainage transport of phosphorus in the midwestern United States: Journal of Environmental Quality, v. 44, p. 495–502, doi:10.2134/jeq2014.04.0176. Soil Survey Staff, 2019, Gridded Soil Survey Geographic (gSSURGO) Database for the Conterminous United States (FY2019 official release): United States Department of Agriculture, Natural Resources Conservation Service (https://gdg.sc.egov.usda.gov/). Soller, D.R., and Garrity, C.P., 2018, Quaternary sediment thickness and bedrock topography of the glaciated United States east of the Rocky Mountains: U.S. Geological Survey Scientific Investigations Map 3392, 2 sheets, doi: 10.3133/sim3392. Sonzogni, W.C., Chapra, S.C., Armstrong, D.E., and Logan, T.J., 1982, Bioavailability of phosphorus inputs to lakes: Journal of Environmental Quality, v. 11, p. 555–563, doi:10.2134/jeq1982.00472425001100040001x. Stow, C.A., Cha, Y., Johnson, L.T., Confesor, R., and Richards, R.P., 2015, Long-term and seasonal trend decomposition of Maumee River nutrient inputs to western Lake Erie: Environmental Science & Technology, v. 49, p. 3392–3400, doi:10.1021/es5062648. Thomas, K.E., Lazor, R., Chambers, P.A., and Yates, A.G., 2018, Land-use practices influence nutrient concentrations of southwestern Ontario streams: Canadian Water Resources Journal / Revue canadienne des ressources hydriques, v. 43, p. 2–17, doi:10.1080/07011784.2017.1411211. Ullman, W.J., and Aller, R.C., 1982, Diffusion coefficients in nearshore marine sediments: Limnology and Oceanography, v. 27, p. 552–556, doi:10.4319/lo.1982.27.3.0552. U.S. Environmental Protection Agency (EPA), 2018, U.S. Action Plan for Lake Erie: U.S. Environmental Protection Agency, 114 p. U.S. Geological Survey (USGS), 2018, National Elevation Dataset. 87 U.S. Geological Survey (USGS) and National Geospatial Program, 2020, USGS National Hydrography in FileGDB 10.1 format (published 20201127). Valiela, I., Costa, J., Foreman, K., Teal, J.M., Howes, B., and Aubrey, D., 1990, Transport of groundwater-borne nutrients from watersheds and their effects on coastal waters: Biogeochemistry, v. 10, p. 177–197, doi:10.1007/BF00003143. Van Meter, K.J., Basu, N.B., Veenstra, J.J., and Burras, C.L., 2016, The nitrogen legacy: emerging evidence of nitrogen accumulation in anthropogenic landscapes: Environmental Research Letters, v. 11, article 035014, doi:10.1088/1748-9326/11/3/035014. Van Meter, K.J., McLeod, M.M., Liu, J., Tenkouano, G.T., Hall, R.I., Van Cappellen, P., and Basu, N.B., 2021, Beyond the mass balance: Watershed phosphorus legacies and the evolution of the current water quality policy challenge: Water Resources Research, v. 57, article e2020WR029316, doi:10.1029/2020WR029316. Van Meter, K.J., Van Cappellen, P., and Basu, N.B., 2018, Legacy nitrogen may prevent achievement of water quality goals in the Gulf of Mexico: Science, v. 360, p. 427–430, doi:10.1126/science.aar4462. Verhamme, E.M., Redder, T.M., Schlea, D.A., Grush, J., Bratton, J.F., and DePinto, J.V., 2016, Development of the Western Lake Erie Ecosystem Model (WLEEM): Application to connect phosphorus loads to cyanobacteria biomass: Journal of Great Lakes Research, v. 42, p. 1193–1205, doi:10.1016/j.jglr.2016.09.006. Verma, S., Bhattarai, R., Bosch, N.S., Cooke, R.C., Kalita, P.K., and Markus, M., 2015, Climate change impacts on flow, sediment and nutrient export in a Great Lakes watershed using SWAT: CLEAN – Soil, Air, Water, v. 43, p. 1464–1474, doi:10.1002/clen.201400724. Vero, S.E., Basu, N.B., Van Meter, K., Richards, K.G., Mellander, P.-E., Healy, M.G., and Fenton, O., 2018, Review: The environmental status and implications of the nitrate time lag in Europe and North America: Hydrogeology Journal, v. 26, p. 7–22, doi:10.1007/s10040- 017-1650-9. Wallace, C.W., Flanagan, D.C., and Engel, B.A., 2017, Quantifying the effects of future climate conditions on runoff, sediment, and chemical losses at different watershed sizes: Transactions of the ASABE, v. 60, p. 915–929, doi:10.13031/trans.12094. Wang, L., Butcher, A.S., Stuart, M.E., Gooddy, D.C., and Bloomfield, J.P., 2013, The nitrate time bomb: A numerical way to investigate nitrate storage and lag time in the unsaturated zone: Environmental Geochemistry and Health, v. 35, p. 667–681, doi:10.1007/s10653-013- 9550-y. Wang, L., Stuart, M.E., Lewis, M.A., Ward, R.S., Skirvin, D., Naden, P.S., Collins, A.L., and Ascott, M.J., 2016, The changing trend in nitrate concentrations in major aquifers due to historical nitrate loading from agricultural land across England and Wales from 1925 to 2150: Science of the Total Environment, v. 542, p. 694–705, doi:10.1016/j.scitotenv.2015.10.127. 88 Wang, Q., and Boegman, L., 2021, Multi-year simulation of western Lake Erie hydrodynamics and biogeochemistry to evaluate nutrient management scenarios: Sustainability, v. 13, article 7516, doi:10.3390/su13147516. Wang, Z., Zhang, T.Q., Tan, C.S., Xue, L., Bukovsky, M., and Qi, Z.M., 2021, Modeling impacts of climate change on crop yield and phosphorus loss in a subsurface drained field of Lake Erie region, Canada: Agricultural Systems, v. 190, article 103110, doi:10.1016/j.agsy.2021.103110. Watson, S.B., Miller, C., Arhonditsis, G., Boyer, G.L., Carmichael, W., Charlton, M.N., Confesor, R., Depew, D.C., Höök, T.O., Ludsin, S.A., Matisoff, G., McElmurry, S.P., Murray, M.W., Richards, R.P., Rao, Y.R., Steffen, M.M., and Wilhelm, S.W., 2016, The re- eutrophication of Lake Erie: Harmful algal blooms and hypoxia: Harmful Algae, v. 56, p. 44–66, doi:10.1016/j.hal.2016.04.010. Wei, X., and Bailey, R.T., 2021, Evaluating nitrate and phosphorus remediation in intensively irrigated stream-aquifer systems using a coupled flow and reactive transport model: Journal of Hydrology, v. 598, 21 p., doi:10.1016/j.jhydrol.2021.126304. Welch, H.L., Kingsbury, J.A., and Coupe, R.H., 2010, Occurrence of phosphorus in groundwater and surface water of northwestern Mississippi: Mississippi Water Resources Conference, p. 142–155. Wellogic, 2005, Wellogic -Statewide Wells: Michigan Department of Environmental Quality, Statewide Groundwater Database (https://www.egle.state.mi.us/wellogic/Login.aspx). Williams, M.R., King, K.W., and Fausey, N.R., 2015, Drainage water management effects on tile discharge and water quality: Agricultural Water Management, v. 148, p. 43–51, doi:10.1016/j.agwat.2014.09.017. World Health Organization, 2011, Guidelines for Drinking-water Quality (4th ed.): Geneva, Switzerland, World Health Organization, 541 p. Wurtsbaugh, W.A., Paerl, H.W., and Dodds, W.K., 2019, Nutrients, eutrophication and harmful algal blooms along the freshwater to marine continuum: WIREs Water, v. 6, 27 p., doi:10.1002/wat2.1373. Wynne, T.T., and Stumpf, R.P., 2015, Spatial and temporal patterns in the seasonal distribution of toxic cyanobacteria in western Lake Erie from 2002–2014: Toxins, v. 7, p. 1649–1663, doi:10.3390/toxins7051649. Yuan, Y., Bingner, R.L., Locke, M.A., Theurer, F.D., and Stafford, J., 2011, Assessment of subsurface drainage management practices to reduce nitrogen loadings using AnnAGNPS: Applied Engineering in Agriculture, v. 27, p. 335–344, doi:10.13031/2013.37075. Zhang, H., Culver, D.A., and Boegman, L., 2008, A two-dimensional ecological model of Lake Erie: Application to estimate dreissenid impacts on large lake plankton populations: Ecological Modelling, v. 214, p. 219–241, doi:10.1016/j.ecolmodel.2008.02.005. 89 Zhang, W., Wilson, R.S., Burnett, E., Irwin, E.G., and Martin, J.F., 2016, What motivates farmers to apply phosphorus at the “right” time? Survey evidence from the Western Lake Erie Basin: Journal of Great Lakes Research, v. 42, p. 1343–1356, doi:10.1016/j.jglr.2016.08.007. 90 APPENDIX Table of all 24 articles examined in the first phase of Chapter 1 (Lake Erie nutrient models). Articles are color-coded by depth, with surface studies in grey and subsurface studies in blue. Studies are described by study area, nutrient form examined, model used, source analyzed, and nutrient quantifications. Abbreviations used: LE = Lake Erie, SWAT = Soil and Water Assessment Tool, ELEMeNT = Exploration of Long-tErM Nutrient Trajectories, OLS = ordinary least squares, SENSEmap = Spatially Explicit Nutrient Source Estimate map, TP = total phosphorus, TDP = total dissolved phosphorus, SRP = soluble reactive phosphorus, DRP = dissolved reactive phosphorus, SP = soluble phosphorus, TN = total nitrogen, DKN = dissolved Kjeldahl nitrogen, DON = dissolved organic nitrogen, ON = organic nitrogen, SN = soluble nitrogen, DIN = dissolved inorganic nitrogen, TDN = total dissolved nitrogen. Table A1. Lake Erie articles in Chapter 1. Table of articles examined in the first phase (Lake Erie) of Chapter 1. Surface or Region/Study Paper Nutrients Model Source/Quantification subsurface Area Forster et al., Subsurface Western LE TP, NO , ON 3 Erosion Productivity No-till nutrient losses compared to 2000 basin Impact Calculator (EPIC) conventional tillage model NO : 5-14% increase in runoff 3 ON: 44-47% decrease in soil TP: 45-46% decrease in soil Liu et al., 2021 Subsurface Grand River TN ELEMeNT-N Estimated N inputs to watershed Watershed Manure: ~23% Atmospheric deposition: ~11% Inorganic fertilizer: ~25% Biological fixation: ~35% Domestic waste: ~6% Fates of cumulative N surplus Exported through stream: 27% 91 Table A1 (cont’d) Denitrification in subsurface: 35% WWTP removal: 10% Legacy N storage in groundwater: 3.9% Legacy N storage in soil: 26.1% Oldfield et al., Subsurface Canadian LE TP Developed geospatial Percent of TP loads from Canada 2020 basin model sourced in septic tanks 1.7 ± 0.8–5 ± 2.3% of P loads and predict future increases Rakhimbekova Subsurface Ontario LE PO , NO , 4 3 Developed geospatial Nutrient loads from septic systems et al., 2021 shoreline NH 3 model MODFLOW- within a distance of the LE NWT shoreline Within 100m: Up to 17 MT/yr of P, up to 110 MT/yr of N Within 250m: Up to 25 MT/yr of P, up to 149 MT/yr of N Van Meter et Subsurface Grand River TP ELEMeNT-P Legacy P accumulation in the al., 2021 Watershed subsurface since 1900 Net P retained: ~96%, with 45% as soil organic P, 34% in strongly sorbed forms, and 21% in weakly sorbed forms P exported to downstream waters: ~4.5% P inputs from manure: 86% from 1900-1945, 59% in 1977 Apostel et al., Surface Maumee River TP, DRP, SWAT Field-scale annual loads with 2021 Watershed TN, NO 3 standard deviations TP: Surface 2.68 ± 2.59 kg/ha, tile 0.04 ± 0.03 kg/ha DRP: Surface 0.29 ± 0.24 kg/ha, tile 0.04 ± 0.03 kg/ha 92 Table A1 (cont’d) TN: Surface 7.10 ± 5.57 kg/ha, tile 25.31 ± 27.04 kg/ha NO : Surface 2.82 ± 3.27 kg/ha, tile 3 25.31 ± 27.04 kg/ha Culbertson et Surface Maumee River TP, SRP SWAT Predicted responses to climate al., 2016 Watershed change by late-century holding fertilizer application rates constant Annual TP runoff: decrease by 7.8- 14.6% Annual SRP runoff: decrease by 8.6-11.8% Predicted responses to climate change by late-century increasing fertilizer application with plant uptake Annual TP runoff: increase by 9.9- 24.6% Annual SRP runoff: increase by 26.7-42.0% Dagnew et al., Surface St. Clair-Detroit TP, DRP SWAT Predicted reductions in response to 2019 River system best-performing management strategies (combo of cover crops, filter strips, and wetlands) 55% spatially random adoption: TP 20-61%, DRP 10-51% 55% targeted adoption: TP 14-37%, DRP 15-61% 100% adoption: TP 32-81%, DRP 23-83% 93 Table A1 (cont’d) Daloğlu et al., Surface Sandusky River DRP SWAT Estimated responses in DRP load 2012 Watershed from baseline to tillage practice changes Conventional: decrease by ~5-20% No-till: increase by ~20-25% Guo et al., Surface NW Ohio DRP, Agricultural Field-scale average annual P loads 2020 (within Maumee particulate P Policy/Environmental DRP in tile discharge: 0.14 kg/ha and Sandusky eXtender (APEX) model DRP in surface runoff: 0.12 kg/ha watersheds) and Nutrient Tracking Particulate P in surface runoff: 0.30 Tool kg/ha Hu et al., 2019 Surface St. Clair-Detroit TP EPA Storm Water Urban area TP load by source River system Management Model and percentage regression models Detroit WWTP: 56% Other point sources: 25% Runoff: 10% CSOs: 9% Detroit River TP load by source percentage Urban point sources: 19% Runoff: 2.5% CSOs: 2.5% Igras and Surface Grand River TP, DRP ISO 31010:2009 Bowtie Probability of achieving ≤ 40% Creed, 2020 Watershed Risk Analysis Tool with a load reduction (TP/DRP) under Bayesian relief network different BMP adoption rates model All BMPs: ≤ 20% adoption = 91.3/29.2% chance, 20-40% adoption = 99.9/76% chance, 40- 60% adoption = 100/96.1% chance, 60-80% adoption = 100/99.1% chance, >80% adoption = 100/99.8% chance 94 Table A1 (cont’d) Commonly used BMPs: ≤ 20% adoption = 17.5/2.1% chance, 20- 40% adoption = 82.2/10.5% chance, 40-60% adoption = 99.6/32.3% chance, 60-80% adoption = 100/50.4% chance, >80% adoption = 100/64.1% chance TP-effective BMPs: ≤ 20% adoption = 71.5/37.3% chance, 20- 40% adoption = 99.4/84% chance, 40-60% adoption = 100/98.2% chance, 60-80% adoption = 100/99.8% chance, >80% adoption = 100/100% chance DRP-effective BMPs: ≤ 20% adoption = 43/26.5% chance, 20- 40% adoption = 93.6/86.6% chance, 40-60% adoption = 100/99.9% chance, 60-80% adoption = 100/100% chance, >80% adoption = 100/100% chance Jarvie et al., Surface Sandusky, SRP Empirical regression Increase in SRP load source 2017 Raisin, and models attribution by percent Maumee Sandusky: 65% increased SRP watersheds delivery, 35% increased runoff Raisin: 63% increased SRP delivery, 37% increased runoff Maumee: 66% increased SRP delivery, 34% increased runoff 95 Table A1 (cont’d) Kujawa et al., Surface Maumee River TP, DRP, TN SWAT Predicted response of loads to 2020 Watershed climate change based on several models’ average TP: 7% decrease, range of 29% decrease to 20% increase DRP: 5% increase, range of 25% decrease to 38% increase TN: 1% decrease, range of 31% decrease to 22% increase Luscz et al., Surface Michigan lower TP, TN SENSEmap Mean loading rates by land use 2015 peninsula Agriculture: TN 118 kg/ha/yr, TP 17 kg/ha/yr Pasture: TN 36 kg/ha/yr, TP 4 kg/ha/yr Urban: TN 31 kg/ha/yr, TP 2 kg/ha/yr Unmanaged: TN 17 kg/ha/yr, TP 0.2 kg/ha/yr Proportion of nutrient loading by source Chemical agricultural fertilizer: TN 43%, TP 59% Manure: TN ~17%, TP ~28% Atmospheric deposition: TN ~33%, TP ~2% Non-agricultural fertilizer: TN ~2%, TP ~2% Septic tanks: TN ~2%, TP ~5% Point sources: TN ~3%, TP ~4% 96 Table A1 (cont’d) Mehan et al., Surface Mason Ditch TP, SP, TN, SWAT Responses of annual nutrient loads 2019 Watershed NO 3 to climate change by the end of the (within Maumee century watershed) Organic P: increase by 15-75% Soluble P: decrease of 60% to increase of 75% Subsurface drain soluble P: decrease by 35-60% Subsurface drain NO : decrease by 3 25-75% Surface runoff NO : increase by 50- 3 100% Merriman et Surface Eagle Creek TP, DRP, SWAT Percent reduction in field-scale al., 2018 watershed TN, NO 3 average annual nutrient surface (within Maumee runoff loads with best-performing watershed) management (conservation cover and filter strips) TP: 94% DRP: 60% Percent reduction in watershed- scale average annual nutrient loads under high BMP implementation TP: 23% DRP: 9.9% NO : ~28% 3 TN: ~32% Nelson et al., Surface Rock Creek TP, TN Agricultural Nutrient load ranges under 2018 Watershed Policy/Environmental differing water management (within eXtender (APEX) model scenarios Sandusky With tile drainage: TP range = watershed) 2.48-4.04 kg/ha, TN range = 26.03- 55.10 kg/ha 97 Table A1 (cont’d) Without tile drainage: TP range = 0.49-2.37 kg/ha, TN range = 5.36- 10.04 kg/ha Thomas et al., Surface 29 SW Ontario TP, TN, OLS regression model Most important land use descriptor 2018 tributaries to LE TDP, SRP, variables in predicting nutrient and Lake Huron DKN, NH , 3 loads and standard coefficients NO /NO , 3 2 DON: WWTP population density, DON, DIN, 0.8688 TDN NH3: Watershed impervious %, 0.7549 DKN: WWTP population density, 0.8461 TDN: Headwater buffer pervious %, −0.9788 TN: Headwater buffer pervious %, −1.0745 NO3/NO2: Stream buffer forage %, −0.5636 DP: 600m buffer forage %, 0.7514 SRP: WWTP population density, 0.8026 TP: WWTP population density, 0.9225 Wallace et al., Surface Cedar Creek TP, SP, TN, SWAT Responses of annual nutrient loads 2017 Watershed SN to climate change by the end of the (within Maumee century watershed) Total P: No significant changes Soluble P: No significant changes to decrease of 25.5% Total N: No significant changes to increase of 19.8% Soluble N: No significant changes 98 Table A1 (cont’d) Wang and Surface Maumee River TP, PO 4 General Lake Model- Load changes under various BMP Boegman, 2021 (for nutrient Aquatic Ecosystem scenarios loading) Dynamics No-till: TP +2.46%, PO −1.54% 4 Cover crops: TP −2.37%, PO 4 −1.66% Filter strips: TP −2.97%, PO 4 −3.08% Combo spatially random: TP −0.88%, PO −4.41% 4 Combo in high source areas: TP −6.92%, PO −8.04% 4 Combo near river mouth: TP −0.86%, PO −0.18% 4 Wang et al., Surface Woodslee, ON DRP Environmental Policy Average responses of annual DRP 2021 (within St. Clair- Integrated Climate loads to climate change by mid- Detroit (EPIC) model century watershed) Total DRP loss: decrease by 11%, DRP in surface runoff: increase by 4% DRP in subsurface drainage: decreased by 22% Wang et al., Surface Woodslee, ON DRP Environmental Policy Proportion of DRP loss sources in 2018 (within St. Clair- Integrated Climate surface runoff Detroit (EPIC) model and Solid manure plots: 19% from watershed) SurPhos model manure, 81% from soil Liquid manure plots: 4.5% from manure, 95.5% from soil Zhang et al., Surface Lake Erie TP, TN EcoLE Percentages of TP and TN in the 2008 water column excreted by dreissenid mussels 1997 P: 23.2% western basin, 1.6% central basin, 1.2% eastern basin 99 Table A1 (cont’d) 1997 N: 12.5% WB, 3.8% CB, 2.8% EB 1998 P: 19.4% WB, 1.1% CB, 0.8% EB 1998 N: 11.8% WB, 3.7% CB, 3.4% EB 1999 P: 55.5% WB, 1.7% CB, 0.8% EB 1999 N: 14.6% WB, 4.3% CB, 2.5% EB 100