Lake hydrodynamics intensify the potential impact of watershed pollutants on coastal ecosystem services

Watersheds deliver numerous pollutants to the coastline of oceans and lakes, thereby jeopardizing ecosystem services. Regulatory frameworks for stressors often focus on loading rates without accounting for the physical dynamics of the receiving water body. Here, we use a three-dimensional hydrodynamic model to simulate the transport of a generic tributary-delivered anthropogenic pollutant within Lake Michigan based on the location and timing of loading. Simulating pollutant plumes from 11 rivers, and their intersections with coastal ecosystem services, reveals strong mediation of potential impacts by lake physics. Trapped pollutants accumulate in nearshore waters during spring peak flows, and become diluted by spreading offshore during the summer. The threat to coastal ecosystem services posed by pollutant loading differs sharply among rivers; high potential impact arises from the spatiotemporal coincidence of tributary input rates, lake mixing dynamics, and multiple human uses of the shoreline. Simultaneous pollution from multiple rivers yields overlapping plumes, creating a second way in which lake hydrodynamics can amplify potential impacts on coastal ecosystem services. Our simulations demonstrate that the physical dynamics of large water bodies can create a dynamic stressor landscape arising from multiple independent sources of non-point-source pollution. The design and implementation of pollution regulations rarely account for spatial and temporal complexities of load processing in receiving waters, yet the resulting variability is likely to strongly mediate impacts on society. As hydrodynamic models improve, our analytical framework could be applied to a wide range of pollutants and waterbodies to enhance the sustainable use of coastal ecosystems.


Introduction
The coastlines of oceans and large lakes are often severely stressed (Millennium Ecosystem Assessment 2005, Diaz and Rosenberg 2008, Halpern et al 2008, Allan et al 2013, in part due to chemical pollution from tributary watersheds. In parallel, human use of large water bodies is concentrated along the coasts, especially near inflowing rivers (Allan et al 2015). The spatial coincidence of tributary-delivered pollution and human reliance on coastlines raises the possibility of widespread impacts on key ecosystem services. Moreover, coastal ecosystems have less dilution capacity than deeper offshore waters, so pollution from watersheds can become concentrated when trapped or transported along the coast (Hoffman and Hittinger 2017). These mixing dynamics of coastal waters are highly dynamic in space and time (Holland andKay 2003, Rao andSchwab 2007), as are the concentration and timing of pollutant inputs delivered by tributary rivers (Dolan and Chapra 2012).
Ideally, limits on pollutant loading to coastal waters should account for the spatial and temporal dynamics of both inputs from watersheds and processing within receiving water bodies. Under the United States' Clean Water Act, limits on pollutant loading are expressed as Total Maximum Daily Loads (TMDLs). While thousands of TMDLs have been developed in the United States (EPA 2019), few are developed for large lakes or marine coastal areas with complex hydrodynamics, with some exceptions such as Chesapeake Bay (EPA 2010). Many of these waters are under the jurisdiction of more than one state, which complicates selection of targets, allocation of pollutant limits, and coordination of technical development and implementation through regulatory programs. In addition, the complexity of modeling pollutant transport and its effect on ecosystem processes and human uses may delay or deter integrating it into regulatory frameworks. For example, the 2012 Great Lakes Water Quality Agreement (United States and Canada 2012) between the United States and Canada sets forth a goal of reviewing and updating the phosphorus loading targets for each Great Lake for the purpose of meeting multiple ecosystem objectives, most of which are related to eutrophication. As of 2019, targets have only been updated for Lake Erie, and the parties have recognized gaps in knowledge of the distribution and movement of nutrients between nearshore and offshore zones in Lakes Ontario, Huron, and Michigan (United States and Canada 2019a).
The Great Lakes of North America offer a particularly ideal setting for evaluating the implications of the physical dynamics of large receiving waters for TMDL-type regulatory paradigms. Along their extensive coastlines, there is enormous spatiotemporal variation in tributary pollutant inputs, hydrodynamics of the receiving water body, and reliance on ecosystem services. The stakes are high for managing pollution of these massive lakes, which hold 84% of North America's surface fresh water (21% of the global total), support nearly 1.5 million jobs, and engender $8B in annual recreational expenditures (Vaccaro and Read 2011). For local communities, including many major cities, these coastlines offer potable water, fisheries, cooling water for power plants, and a wide range of aesthetic and recreational benefits (Allan et al 2015). Threats to these ecosystem services from harmful algal blooms, hypoxia, toxic chemicals, and species invasions inspired the ongoing Great Lakes Restoration Initiative (GLRI) in the United States, which spent $2.3B from 2010-2017 to enhance ecosystem health (EPA 2017). Nonetheless, closures of swimming beaches and municipal drinking water intakes remain a regular occurrence due to pollution from industrial facilities, agricultural landscapes, and urban centers.
TMDL-type regulatory paradigms often focus on the sources and dynamics of pollution, but they also have the potential to incorporate understanding of how lake physics mediate the propagation of pollutants in space and time within large aquatic ecosystems. The depth and thermal structure of lakes creates three-dimensional heterogeneity in water density that can foster either rapid dilution by mixing or the trapping of a pollutant entering from a particular location (Auer and Gatzke 2004, Makarewicz et al 2012, Yurista et al 2015. Each of the Great Lakes shows a pronounced annual cycle in its thermal structure (Boyce et al 1989, Beletsky andSchwab 2001), changing from vertically mixed in winter to vertically stratified and laterally mixed in summer. Prior to the onset of summer stratification, spring surface warming of the shallow water column along the coast separates it from cooler water offshore. The convergence zone between the inshore and offshore water masses occurs at the temperature of maximum density (~4 • C), and this 'thermal bar' (Holland andKay 2003, Bai et al 2013) acts as a barrier to their mixing. As the surface continues to absorb radiation and warms, the thermal bar dissipates, and is replaced by stable vertical stratification across the lake during the summer and early fall. These strong vertical versus lateral gradients of water temperature-and thus density-modify circulation patterns within each Great Lake (Gbah and Murthy 1998, Auer and Bub 2004, Rao and Schwab 2007, Makarewicz et al 2012 and can influence ecological functioning (Scavia and Bennett 1980). Both observations and models have illustrated the importance of these hydrodynamic patterns for coastal trapping of nutrients and contaminants (Spain et al 1976, Boyce et al 1989, Auer and Bub 2004, Auer and Gatzke 2004. The Great Lakes have been the subject of relatively comprehensive analyses of the geography of both stressors and ecosystem services. Examples of tributary derived stressors include non-point sources such as suspended sediments, nutrients, and combined sewer overflow events or chemical pollution such as mercury, copper, and polychlorinated biphenyl compounds (PCBs). A synthesis of 34 anthropogenic stressors shows that cumulative stress is highest in coastal areas with high levels of human activity on shore and in upstream watersheds (Allan et al 2013). Similarly, Great Lakes coastlines are the focal point for many aspects of human use, representing high but spatially variable reliance on ecosystem services (Allan et al 2015). Thus, the Great Lakes offer a unique opportunity to integrate realistic spatiotemporal dynamics of pollutant inputs and in-lake transport with ecosystem service distributions to evaluate whether potential impacts on ecosystem services are closely related to watershed loading rates. This integrative approach could enable pollution control programs to direct their efforts toward watersheds where loading patterns and lake hydrodynamics create the highest potential societal benefits from reducing coastal ecosystem degradation. Here, we simulate the fate of watershed-derived pollutants in Lake Michigan. Lake Michigan has a drainage area of 118 000 km 2 , a retention time of 99 years, an average depth of 85 m, and a volume of 4920 km 3. We couple a three-dimensional advection-transport model of Lake Michigan forced with realistic meteorology with input time series of a passive, generic tracer from 11 different tributaries (figure 1). Tracer input time series are based on phosphorus, an important anthropogenic pollutant. However, we are not attempting to simulate phosphorus concentrations which are, in the real lake, modulated by a variety of biological processes. Instead, we use this time history to represent a generic tributary-derived pollutant that would subsequently be acted on by lake hydrodynamics. Within the lake, our generic pollutant is modeled as a conservative tracer that could represent any contaminant entering the lake via time-varying river flows. Recent large-scale stressor assessments in the Great Lakes have overlooked the spatial and temporal variability of watershed-derived stressors, despite their contribution to aggregate ecosystem stress status (Allan et al 2013). After developing realistic watershed input estimates, we evaluate whether the seasonal hydrodynamics of a large receiving body lead to trapping, dilution, or transport of pollution. In addition to testing how the capacity of the lake to attenuate pollutant concentrations varies with the location of the river mouth, we overlay the pollution plumes from all 11 watersheds. Finally, we intersect the aggregate pollution distribution with maps of summertime use of various coastal ecosystem services. These spatial comparisons enable us to quantify the cumulative stress to services based on the duration and spatial extent of overlap between elevated pollutant concentrations and coastal ecosystem service locations, thereby focusing on the potential implications of pollution rather than its intensity perse.

Hydrodynamic model description
We used the Massachusetts Institute of Technology general circulation model (MITgcm) (Marshall et al 1997a(Marshall et al , 1997b configured to the bathymetry of Lake Michigan (National Geophysical Data Center 1996). For computational efficiency, the lake is modeled as a closed basin with a horizontal resolution of 1 min (approximately 2 km) and 28 vertical levels. The thickness of vertical levels increases from 5 m at the surface to 31 m at the bottom. Sub-grid-scale processes are simulated via the K-Profile Parameterization (KPP) vertical mixing scheme (Large et al 1994) and the Smagorinsky horizontal diffusivity scheme (Smagorinsky 1963). Further discussion of the physical model and its validation can be found elsewhere (Bennington et al 2010, Pilcher et al 2015. The model solves the tracer transport tendency equation at each grid cell (equation (1)) where (U · ∇) C represents advection of a tracer volume concentration (C) by the three-dimensional flow field (U), S includes convective mixing and sub-grid parameterizations of mixing, and I represents time-varying tracer input from 11 tributaries (figure 1). The tracer is neutrally buoyant, biologically inert, and non-interactive with the atmosphere, suspended particles, and bottom substrate. This approach allows us to focus on pollutant trapping and redistribution via hydrodynamics. The model is forced at the surface using meteorological and radiative fields with a temporal resolution of 3 h. Atmospheric downward shortwave and longwave radiation, air temperature, specific humidity, and wind stress are imposed at surface from the North American Regional Reanalysis Project (NARR) (Mesinger et al 2006). Daily lake ice fractions are taken from the U.S. National Ice Center (U.S. National Ice Center 2010). There is an established warm bias in the physical model due to both the NARR forcing and the ice mask (Bennington et al 2010, Pilcher et al 2015. Despite this known bias, NARR is selected due to its higher spatial resolution relative to other choices. The model is spun up for two years by repeating 2007 forcing to achieve a repeating annual cycle (Pilcher et al 2015), after which it is run for a four-year period (2007)(2008)(2009)(2010)(2011). We focus our analysis on surface concentrations between May and August (figure 2).
Input of our generic tracer is based on estimated daily histories of total phosphorus export for 11 rivers (supplementary figure 1 is available online at '(stacks.iop.org/ERL/15/064028/mmedia)') from 2007-2010. These estimates were derived from the weighted regressions of concentrations on time, discharge and season method (WRTDS, supplementary table 1) (Hirsch et al 2010). To validate our input physics, we assess correlation to bi-weekly observations of surface phosphorus concentrations near Milwaukee, and find that the model reasonably estimates the temporal variation of in-lake concentrations at this location (r = 0.7; supplementary figure 2). This indicates modeled physical dynamics for input and near-shore redistribution are realistic. Our generic pollutant is conservative, in contrast to phosphorus that is quite reactive. By making this comparison to phosphorus data collected close to the input location, the effect of our neglect of reactions is minimized. In the model, tracer concentrations in the lake are reset to zero on 1 January of each year to reflect overturn during the cold season and because our focus is to assess the annual pattern of in-lake plume formation and dispersion. For our results, pollutant concentrations arising from tributary loading were scaled by dividing by the expected background concentration from fully mixing the total annual load from all rivers into the entire lake volume (~1 µg l −1 ). This allows regions above a well-mixed background concentration to be easily identified.

Ecosystem services
Maps of ecosystem service delivery on a 2 × 2 km grid were drawn from the Great Lakes Environmental Assessment Mapping (GLEAM) Project (Allan et al 2013(Allan et al , 2015. From the 34 stressors assessed in GLEAM, we select 3 that are most illustrative of the human uses of the coastline that will be impacted by plumes of tributary pollutants. These maps indicate the location-based usage of each service (municipal water intake, beaches, and marinas/boat launches) (figure 3), but do not represent the economic value or supply-demand balance of each service. Following Allan et al (2015), each usage metric is log 10 (x + 1) transformed, and then linearly normalized between the maximum and minimum value. These timeinvariant services maps are primarily representative of usage in the summer months (June-August).
In the GLEAM dataset, service location and usage are mapped separately for each service using publicly available data. Using data maintained by the Great Lakes Commission and individual state contacts, municipal water usage is quantified as the annual (2013 or 2014) water withdrawal amounts in millions of gallons per year from intake/outfall sites within 5 km of the shore. Beach locations are obtained from the US EPA BEACH Act Geospatial database. As a proxy for beach visitation, the InVEST model (Natural Capital Project 2013) counts the number of geo-tagged Flickr photos within a 500 m buffer of each beach location (see Wood et al (2013)). Marina/boat launch locations are identified from internet sites, governmental sources, tourist information publications, and from marinas that reported boat launch availability; locations are confirmed using Google Earth. Usage for marina and boat launches is obtained from marina websites, agency sources, and Google Earth imagery. Marina and boat launch usage is quantified as the number of boat slips and boat launch parking spaces, respectively.

Cumulative Stress Days (CSD)
Cumulative Stress Days (CSD) quantify the intersection of service usage and tributary-delivered pollutant concentrations, and is calculated by multiplying the number of days exceeding the background concentration (Exceedance Days, ED) by normalized service usage at adjacent grid cells. Our CSD metric is designed to quantify the cumulative stress to services caused by local exposure to elevated pollutant concentrations. Since the background concentration is equal to the annual load being well mixed throughout the entire lake, this identifies times and locations where pollutant concentrations are elevated due to inlake hydrodynamics.
To calculate CSD, we first calculate ED separately for each river and for all the rivers grouped in the North, East, and West. ED is a proxy for potential stress to ecosystem services, and is defined as the period over which pollutant concentrations (C) exceed an estimated background concentration state (  C bg  ) . At each location, and for each day of the year, a value of 1 is assigned if C is greater than C bg and 0 otherwise. Then, for a particular nutrient load source, we tally these binary (0,1) scores for each grid cell over the M days in each month to get a monthly ED (equation (2)).

ED =
Finally, for each month, each service metric (Boating, Beaches, or Water) is weighted by ED at each of the N locations and summed to yield the CSD for that month and for the respective region or river (equation (3)).

CSD =
Large CSD values can arise by persistently high tracer concentrations within the region, many days of residence of the tracer plume on the shore, or the density of ecosystem services. Changes in service supply or demand associated with an excessive concentration are not accounted for with CSD. Other metrics that weigh the absolute concentration of nutrients from tributary loading more than persistence of the threshold were considered. Since both types of metrics give qualitatively similar results, we chose the simpler ED and CSD metrics.

Seasonal barriers to mixing
Our simulations of time-varying tributary loads and lake mixing reveal that both strongly affect the spatial (figure 1) and temporal (supplementary movie) propagation of pollutants arriving from each watershed. The strong lateral temperature gradient-the 'thermal bar'-in May acts as a barrier to horizontal transport (figure 2(A)), thereby trapping pollution near the coast and acting as a conduit for spreading along the coast ( figure 2(B)). As spring progresses, the lake surface warms everywhere and the lateral barrier is diminished. In its place, vertical thermal stratification develops (supplementary figure 3). The progression toward a laterally-uniform surface temperature (figure 2(C)) allows the pollutant to spread offshore and become diluted (figure 2(D)). During the late summer and fall, vertical stratification traps the pollution above the thermocline without constraining its lateral spread (supplementary figure 3).
When all 11 tributary plumes are overlaid, we find a highly heterogeneous distribution of aggregate pollutant concentrations (figures 1, 2(B), and (D)). Peak concentrations occur in regions where inputs from multiple tributaries coincide. Monthly average pollutant concentration is particularly elevated along the eastern coastline and within Green Bay  (figures 2(A) and (B), supplementary figure 4). Peak coastal concentrations occur in May due to the influence of both the thermal bar and high input rates associated with annual maximum discharge from each tributary ( figure 2(B)). Interannual variation in wind patterns creates substantial differences in the dispersal of tributary plumes from year to year (supplementary figure 4).
To represent the point beyond which aggregate pollutant concentrations might begin to threaten ecosystem services from the lake, we developed an Exceedance Days metric (ED, see Methods). ED is a tally of the number of days per year when local accumulation of the conservative pollutant leads to concentrations that exceed expectations from mixing the same total load throughout the lake. We find that ED is particularly high within Green Bay, but is also consistently higher in the nearshore zone than in open water. There is also substantial temporal variability in ED arising from the dynamics of plume location and concentration. For example, ED near the Sheboygan River is elevated in June as its pollution load remains trapped along the shore (figure 3(D)), then drops sharply in July and August as the plume mixes laterally (figures 3(F) and (H)).

Potential stress on coastal ecosystem services
The distribution of coastal ecosystem service utilization (figure 3) is no less spatially heterogeneous than that of tributary-derived pollution (figure 1). Elevated pollutant concentrations intersect with many important ecosystem service sites during the summer (supplementary figure 5). To integrate across boating, beaches, and water intakes as potentially impacted services, we calculated a Cumulative Stress Days metric (CSD, see Methods) that sums the ED for each tributary plume that directly contacts coastal sites where these three ecosystem services are documented.
Ecosystem services in the North and East regions of Lake Michigan are exposed to more tributary-derived pollution than services in the West (figure 4(A)). High values of CSD can arise from either of two distinct patterns: sustained high pollutant concentrations that overlap with a modest number of services, or intermittently high concentrations intersecting with a larger number of services. For instance, although Green Bay's load from the Fox River drives average summer pollutant concentrations higher in the North (5.1 µg l −1 ) compared to the East (0.7-1.1 µg l −1 ), the river plumes in the East overlap spatially with twice as many services (supplementary figure 5). These intermittent high concentrations at sites of multiple services leads to peak CSD in the East ( figure 4(A)). In contrast, high CSD values in the North are due to prolonged high concentration of pollution at a few service locations (supplementary figure 5).
Within each region, there is variability in CSD across the summer months that reflects the dynamic extent of river plumes in the vicinity of service sites (figure 3). For example, low ED in June leads to minimal CSD at Washington Island (figure 3(E)), while CSD generally increases later in the summer as plumes expand across the North region (figures 3(G), (I) and 4(A)). In contrast, pollution tends to be trapped along the coast in the East (figure 3(A)) and West (figure 3(D)) during June, thereby creating elevated CSD at service sites ( figure 4(A)). In July and August, ED values are reduced (figure 4(A)) in the East (figures 3(B) and (C)) and West (figures 3(F) and (H)) as river plumes are diluted by mixing with open water. However, high ED continues and expands within Green Bay throughout the summer due to its lack of dilution capacity and limited exchange with the open lake (figures 3(E) and (H)).
Calculating CSD separately for each river plume underscores the roles of both loading rates and ecosystem service distributions in governing potential impacts on societal use of Lake Michigan. CSD in the West region is predominantly driven by the Milwaukee River. However, in June the Sheboygan River plume spreads further along the shore (figure 3(D)), and makes a greater contribution to CSD (supplementary figure 5). The number and spatial distribution of services is limited in the North compared to other regions (figure 3, supplementary figure 5), but CSD is high due to high ED from persistent inputs of pollution from the Fox River (Dolan and Chapra 2012) into shallow Green Bay, which mixes slowly with the open lake. In the East, ED is dominated by the Grand, Kalamazoo, and St. Joseph Rivers (figure 4(A)), each of which drives high CSD values as its plume overlaps with many service sites along the coastline ( figure 3, supplementary figure 5). These three river plumes spread along the shoreline in June due to the thermal bar ( figure 3(A)), giving rise to high CSD values (figure 4(A), supplementary figure 5). Moreover, elevated CSD values to the north of the Muskegon River reflect overlapping loading from four rivers that each discharge to a stretch of the Lake Michigan shoreline where coastal currents carry pollution to areas that people use regularly for boating, beaches, and water supplies (supplementary figure 6). The intensity of this joint pollution is alleviated by offshore mixing later in the summer, but not entirely eliminated ( figure 4(A)).

Discussion
Merging realistic riverine pollutant loading and three-dimensional lake hydrodynamics with ecosystem service maps offers a new perspective on pollution control strategies for large water bodies, which can lead to more effective coastal management approaches. Physical mixing within Lake Michigan yields extensive redistribution of watershed-derived pollution during our study period (2007)(2008)(2009)(2010). Strong offshore temperature gradients, such as the thermal bar, in spring and early summer inhibit offshore mixing, thereby trapping pollution against the shore ( figure 2). Notably, the resulting high expected pollutant concentrations at that time are primarily an effect of in-lake physics rather than high loading rates, even though nutrient control measures often focus on load dynamics alone. We find that, regardless of the instantaneous loading rate, the thermal bar not only causes incoming pollution to accumulate in the vicinity of the river mouth, but also promotes its spread along the shoreline, where societal reliance on ecosystem services is greatest (Allan et al 2015). This pattern is readily visualized through an animation of the model results (see supplementary movie), including the transition from pollution being trapped by the thermal bar (April 1st) to spreading widely under vertical stratification as the lake surface becomes isothermal during the summer (August 1st).
Our findings illustrate the insights that can emerge from unifying spatial analyses of stressors and ecosystem service use with models of the hydrodynamics of the receiving water body. As mixing dynamics cause pollution to be trapped and spread throughout the nearshore zone (figure 2), the risk that watershed-derived pollution will affect coastal ecosystem services becomes amplified. The effectiveness of this trapping, when it coincides with the locations of population centers and ecosystem service utilization, could lead to problems with harmful algal blooms, bacterial contamination, and other water quality problems even when the load arriving from the nearest tributary is not particularly high ( figure 4(B)). Indeed, the East region of Lake Michigan exemplifies this scenario; the likelihood of impacts on services arises from the overlap of multiple river plumes, each of which could only modestly increase pollutant concentrations in lake waters when considered in isolation.
Our incorporation of lake hydrodynamics into analysis of the spatial relationships between pollutant inputs and human use of Lake Michigan's shoreline is merely a start toward quantifying the impact of pollution on coastal ecosystem services (Dodds et al 2008, Wolf et al 2017, but it has clear implications for regulations in the Great Lakes. The first phosphorus loading goals for the Laurentian Great Lakes were expressed as total loads (International Joint Commission 1978), an approach that implicitly treated all sources and locations equally. Though overall phosphorus loads to Lake Michigan have decreased below the target concentration under the existing binational agreement (Dolan and Chapra 2012), problems with nuisance algae remain in many embayments and near-shore areas (Brooks et al 2015) due to elevated coastal phosphorus (Yurista et al 2015). Algal growth in the nearshore zone is exacerbated by invasive zebra and quagga mussels that capture planktonic nutrients and deposit them on the substrate (Hecky et al 2004, Mosley and Bootsma 2015, Pilcher et al 2017, Rowe et al 2017, but lake hydrodynamics plays a mediating role even in that process (Pilcher et al 2015). Indeed, the ironic state of Lake Michigan is that invasive filter feeders and reduced phosphorus loading have diminished open-water phytoplankton production to the point of undercutting the energetic basis for important pelagic fisheries (Rowe et al 2017), yet coastal ecosystems still suffer from widespread overabundance of benthic and planktonic algae (Brooks et al 2015). The contribution of temperature-driven trapping of watershed nutrient loads along the shoreline has received little attention, but our simulations suggest that it may be an important factor.
The continuing evolution of regulatory frameworks that set pollutant loading targets offers opportunities to incorporate the hydrodynamics of receiving bodies. A focus on total loads, even when specific to a particular tributary, can inadvertently overlook both pollutant accumulation in the nearshore (Yurista et al 2015) and advection to far-away sites that are critical for coastal ecosystem services (Hoffman and Hittinger 2017). An approach like ours can be used to estimate the spatial and temporal variation in the effective load to the coastal zone under a range of watershed management and climate scenarios, thereby enabling the next generation of regulatory load targets to account for seasonality of both inputs and in-lake processing. As the United States and Canada launch efforts to update phosphorus load targets for Lake Michigan and other Great Lakes in the next few years (United States and Canada 2019b), we hope that the growing evidence that lake hydrodynamics either trap or dilute tributary loads will be considered. While it may be too cumbersome to expect regulatory agencies to run sophisticated hydrodynamic models for every receiving body, our findings suggest that setting targets without accounting for mixing dynamics of large lakes and coastal oceans would undercut the objective of protecting public health and environmental quality.
Regulatory mandates are often framed around avoiding the impairment of human use of water bodies, and our spatial analysis of the intersection of aggregate pollution patterns with coastal ecosystem services represents a feasible approach for tackling this issue explicitly. However, we recognize that our Cumulative Stress Days (CSD) metric is a blunt instrument for translating pollutant concentrations into potential impact on services. Its strengths include integrating through time and across multiple services, but it focuses on exceedance of an estimated background concentration in a well-mixed system as a proxy for impact. Refining this approach could start with integrating the sensitivity of ecosystem service levels to the relative dosage of key pollutants, as well as differentiating between reactive and conservative behavior of pollutants. It also would be ideal to account for potential synergistic effects of multiple pollutants entering from the same tributary (Smith et al 2019).
As eutrophication and hypoxia continue to expand along the world's coastlines (Smith 2003, Diaz andRosenberg 2008), it is critical to bear in mind both successes and limitations in managing loads of nutrients and other pollutants. The current paradigm in the United States of regulating loads from each watershed independently has unquestionably had positive effects over the last 50 years, yet hypoxia continues to plague commercial fisheries in the Gulf of Mexico (Scavia et al 2017) and nuisance algal blooms are expanding in lakes of all sizes (Dolan andChapra 2012, Brooks et al 2015). These forms of environmental degradation threaten commerce and quality of life for coastal communities (Dodds et al 2008, Allan et al 2015, Wolf et al 2017, and partly reflect a failure to recognize that the internal dynamics of large water bodies often redirects pollutants. Our simulations of Lake Michigan illustrate the insights into nutrient management that can come from jointly accounting for spatiotemporal patterns of pollutant inputs, large-scale hydrodynamics, and human use of coastal ecosystems. This analytical approach is transferable to any large water body, and incorporating it into the next generation of regulatory loading targets is likely to enhance the protection of critical coastal ecosystem services.

Author contributions
L G, G A M, R J M, and P B M designed the study. L G and G A M adapted the hydrodynamic model for this study. J D A and P B M guided use of the ecosystem service maps. M W D and R J M estimated the tributary phosphorus flux histories. L G performed all simulations and analyses, and drafted the manuscript. All authors discussed the results and refined the manuscript.