Management Transition to the Great Lakes Nearshore: Insights from Hydrodynamic Modeling

<jats:p>The emerging shift in Great Lakes management from offshore to nearshore waters will require attention to complexities of coastal hydrodynamics and biogeochemical transformations. Emphasizing hydrodynamics, this work resolves transport processes in quantifying discharge plume and pollutant of concern (POC) footprint dimensions, the latter being the portion of the plume where water quality standards are not met. A generic approach, isolated from pollutant-specific biokinetics, provides first-approximation estimates of the footprint area. A high-resolution, linked hydrodynamic-tracer model is applied at a site in the Greater Toronto Area on Lake Ontario. Model results agree with observed meteorological and hydrodynamic conditions and satisfactorily simulate plume dimensions. Footprints are examined in the context of guidelines for regulatory mixing zone size and attendant loss of beneficial use. We demonstrate that the ratio of the water quality standard to the POC concentration at discharge is a key determinant of footprint dimensions. Footprint size for traditional pollutants (ammonia, total phosphorus) meets regulatory guidelines; however, that for soluble reactive phosphorus, a presently unattended pollutant, is ~1–2 orders of magnitude larger. This suggests that it may be necessary to upgrade treatment technologies to maintain consistency with regulatory guidelines and mitigate manifestations of the eutrophication-related soluble reactive phosphorus POC.</jats:p>

section_content: Located at the watershed-waterbody interface, coastal margins are both the primary standard for stakeholder perception of water quality and a bellwether of future conflict in ecosystem services. The Laurentian Great Lakes include 9783 miles of coastal margin in the United States (US) and Canada, and collectively have long been referred to as "the inland seas" [1]. North America's inland seas are likened to ocean systems, as coastal dynamics in both are driven by similar physical processes (e. g., coastal jets [2] [3] [4] ). Unlike the coastal oceans, however, the Great Lakes are closed systems where the influence of hydrodynamic phenomena on ecosystem function is magnified [4]. For example, discharges to nearshore waters may be entrained by coastal jets, extending the region of influence along the broader coastal zone beyond the point of entry and into the larger lake [4] [5] [6]. Similarities between the coastal margin of marine and large lake systems are such that similar challenges are faced [7] and common hydrodynamic modeling platforms are often utilized in simulating mass transport (e. g., Princeton Ocean Model (POM) [8], Finite Volume Community Ocean Model (FVCOM) [9, 10], Estuary and Lake Computer Model (ELCOM) [11] ). In regulatory applications, engineering design guidance for ocean outfalls and offshore discharges to freshwaters share a common manual [12]. As a result of these similarities, contributions to the body of knowledge on physical processes in large lakes merit a place in the coastal ocean literature (historically [3] and more recently [13, 14] ).
section_name: Introduction

section_content: Management of pollutants of concern (POCs) in the nearshore calls for an approach that is radically different from that applied to open lakes and offshore waters. The open lake is often conceptualized as a completely mixed system where, upon input, individual loads are mixed instantaneously across the entire water body, resulting in only a modest increase in lakewide pollutant concentration and attendant impact (Figure 1a ). However, this conceptualization translates poorly to the nearshore. Coastal zone management recognizes the system as being incompletely mixed as tributary and point source inputs are received at locations proximate to the land-water boundary. These inputs are retained in the nearshore for extended periods, potentially becoming juxtaposed with receptors (e. g., beaches, algal habitat, water intakes), while undergoing dilution and biokinetic transformation prior to transport to the open lake (Figure 1b ). This conceptualization of the nearshore highlights the importance of nutrient loading reductions and separation of nutrient-rich point source effluent plumes from algal habitat in ecosystem management, with common approaches including advanced point source treatment and extension of offshore discharges (see Reference [15] ). and (c ) an effluent plume, over which a dilution gradient (dashed lines) is manifested, as well as the footprint for a pollutant of concern (gray shading). The extent of the plume is set arbitrarily at the 99. 5% dilution isopleth and represents the confluence of water containing effluent and the unimpacted boundary condition. The extent of the footprint is set at the 98% dilution isopleth for this example.
section_name: Science, Engineering, and Management Context

section_content: The western shoreline of Lake Ontario is characterized as "an almost unbroken urban and (c ) an effluent plume, over which a dilution gradient (dashed lines) is manifested, as well as the footprint for a pollutant of concern (gray shading). The extent of the plume is set arbitrarily at the 99. 5% dilution isopleth and represents the confluence of water containing effluent and the unimpacted boundary condition. The extent of the footprint is set at the 98% dilution isopleth for this example. 

Point source engineering design for achieving the required separation between nutrient-rich wastewater effluent plumes (Figure 1b, c ) and receptors begins with the definition of discharge plume dimensions. The plume boundary marks the location at which constituent concentrations in plume waters approach those of ambient waters, illustrated in Figure 1c as the 99. 5% dilution isopleth. The footprint is located within the plume and represents the region in which the POC concentration exceeds the water quality standard (illustrated in Figure 1c as the 98% dilution isopleth). In regulatory applications, a mixing zone is defined as a region where compliance with water quality standards is not required [16]. However, regulatory guidelines state that mixing zones are to be kept as small as possible, not interfere with beneficial uses, and not be used as an alternative to reasonable and practicable treatment [16]. Where a POC footprint conforms to these guidelines, regulatory objectives are met, water quality management goals are realized, and the terms footprint and mixing zone may be considered synonymous. Where guidelines are not met, the discharge must be moved further offshore (enhancing separation from receptors) or additional treatment (e. g., reducing discharge effluent concentration) is required. The focus of this work is to characterize plume and footprint dimensions for tributary and point source inputs and better understand their impacts on coastal zone management.
section_name: Study Site

section_content: The western shoreline of Lake Ontario is characterized as "an almost unbroken urban landscape" extending from the mouth of the Niagara River, across the Greater Toronto Area (GTA), to Oshawa, Ontario (ON) [17] (Figure 2 ) and is today home to a population approaching seven million. Eastward along the north shore, land use is dominated by agriculture with some forest cover. A gradient in phosphorus enrichment and attendant water quality parallels the gradient in urbanization [18] [19] [20] [21]. 

In the lesser developed region to the east of the GTA, that gradient is overlain by a mosaic of discharges from generally small tributaries and wastewater treatment plants [21]. While consensus was built regarding the significance of effluent plumes, a critical knowledge gap remains in establishing plume dimensions (i. e., longshore and offshore extent) and their residence time in the nearshore. 

The study site is located along a 25-km stretch of the Lake Ontario nearshore from Pickering to Whitby, ON, centered on Ajax, ON and extending offshore for a distance of ~3 km (4. 5 km at two stations; Figure 2 ). Although positioned near the eastern boundary of the GTA, the study site remains highly developed. Maximum depths range from 20 to 35 m, the approximate point of separation for nearshore and offshore waters in Lake Ontario [22] (see also References [5, 6] ). Three tributaries (the Rouge River, and Duffins and Carruthers Creeks), three wastewater treatment plants (WWTPs; Highland Creek, Duffin Creek, and Corbett Creek), and the Pickering Nuclear Generating Station (PNGS; cooling water intake and discharge) are located within the study site (Figure 2 ; Figure 3 ; Table 1 ). Based on total phosphorus (TP) load [19], the Duffin Creek, Highland Creek, and Corbett Creek WWTPs are, respectively, the second, seventh, and ninth largest of the 38 plants discharging to the Canadian waters of Lake Ontario [19]. The Pickering Nuclear Generating Station withdrawal and return represents use for once-through cooling water. While this intake and discharge impacts local flow conditions, it has no impact on water quality constituents. Management of POCs in the Pickering → Ajax → Whitby nearshore is a source of public concern [23]. For example, this region experiences nuisance growth of the attached, filamentous, green alga Cladophora with episodic beach fouling leading to loss of beneficial uses [24]. Nuisance growth is stimulated by sources of soluble reactive phosphorus (SRP), the freely and fully available form of the nutrient [24, 25]. The western portion of our site was the subject of field and hydrodynamic water quality modeling studies in 2007 and 2008 [11]. The authors reported that the outfall from the Duffin Creek WWTP was responsible for elevated SRP levels noted in both field monitoring and modeling results. It was concluded that, if model-predicted episodic elevations in SRP (some not observed in monitoring data) were real, they could have a significant stimulatory impact on Cladophora. We believe that the episodic events were caused by incursions of SRP-rich WWTP effluent across their study area, reflecting the highly dynamic nature of the nearshore environment. Thus, our interest is in more fully capturing plume dimensions and characterizing pollutant of concern footprints for the tributary and point source discharges to the Pickering → Ajax → Whitby nearshore. This location was also the focus of studies quantifying phosphorus loads [26] and characterizing the Cladophora-phosphorus dynamic through monitoring [11, 21, [27] [28] [29] [30] [31]. 

Table 1. Tributary and point source discharges included in the study. Tributary flows are summer (May-September) averages, while those for WWTPs are annual averages for 2014. Data sources are footnoted here and described in detail below. Flows provided here facilitate characterization of study site hydrology, while those for specific periods used in model runs are provided in Table 2.
section_name: Study Site

section_name: Flow (m

section_content: More than 15 years ago, Rao et al. [17] warned that the technologies being utilized at wastewater treatment plants discharging to western Lake Ontario were nearing the limits of their design and noted that "the demand for clean water and the need for suitable waste disposal facilities continue to rise at an ever-increasing rate". Since that time, the population estimated at four million by Rao et al. [17] rose to almost seven million and is projected to approach 10 million by 2041 (Statistics Canada, 2017; https://www12. statcan. gc. ca/census-recensement). Yet, the treatment technologies applied across the region (e. g., chemically enhanced secondary treatment for phosphorus) remain those introduced in the 1970s and 1980s. From this, one may posit that a shift in management focus to the Lake Ontario nearshore will reveal a necessity to reassess the impact of POC plumes with respect to regulatory guidelines for coastal waters. We hypothesize that the footprints for previously unattended and emerging POCs may be much larger than traditional POCs and need to be effectively reduced to maintain compliance in a regulatory context. 

We test this hypothesis through application of a hydrodynamic model facilitating comparison of effluent and tributary plume and POC footprint dimensions for traditional (e. g., total ammonia nitrogen (TAN) and total phosphorus; see [17] ) and a previously unattended POC (soluble reactive phosphorus). Management of POCs requires attention to both mass transport (hydrodynamic processes, unrelated to the type of POC) and biogeochemical (biokinetics, POC-specific) transformations. This work focuses on the first of these, mass transport, with the intent to characterize POC footprint maxima. Therefore, the target POC is treated as a conservative substance, an approach appropriate for determination of discharge plume dimensions and POC footprint maxima solely as mediated by mass transport, i. e., in isolation from biokinetic processes. 

The hydrodynamic model was developed, and simulations were performed and applied in three steps as follows:

1. 

Test model performance by comparing output to field measurements to establish the credibility and reliability required to support management decisions; 2. 

Characterize dynamic plume behavior (e. g., variation in the direction of alongshore movement) as a determinant of plume dimensions; 3. 

Evaluate POC footprint dimensions for their consistency with mixing zone guidelines for nearshore water quality management. 

In testing our hypothesis, we seek to provide guidance for the transition from open lake to nearshore management of previously unattended and/or emerging POCs in the Great Lakes. We believe that our approach and our findings will find application in other large lake and coastal ocean environments exploring the need for next-generation effluent management strategies.
section_name: Objectives and Approach

section_name: Methods

section_content: Nearshore and coastal waters are highly dynamic, thus posing a challenge with respect to the development of hydrodynamic models. For example, the ability to simulate transport of effluent plumes is highly dependent on the accuracy of modeled current flow. An acoustic Doppler current profiler (ADCP) providing continuous current flow measurements was deployed from July to September 2014 (see Figure 3 for location). ADCP data were processed and analyzed along with meteorological measurements obtained from National Data Buoy Center (NDBC) Buoy 45159 (Figure 3 ) for validation of hydrodynamic model input forcing and simulation results. In addition, to characterize the spatial variability of effluent plumes and provide ground truth for model calibration and confirmation, we also conducted water quality sampling on 19 July and 9 September 2014. We sampled along 14 transects (86 stations) in July and 12 transects (71 stations) in September, spanning the 25-km breadth of the study site (Figure 3 ). A YSI 6600 V2-4 multi-parameter water quality sonde and an ISUS V3 sonde were used to obtain surface-to-bottom (1-20 m) profiles of temperature ( • C), conductivity (µS/cm), and nitrate (as nitrogen, NO 3 -N) at each station.
section_name: Field Sampling and Laboratory Methods

section_content: The research described here utilizes the Finite Volume Community Ocean Model (FVCOM), a platform widely applied in estuarine and coastal ocean environments [32] [33] [34] and emerging as a commonly used tool for application to the Great Lakes [9, [35] [36] [37] [38] in both stand-alone hydrodynamic modeling and coupled with a water quality [39] or regional climate model [10, 40]. In this section, we describe the methods underlying linkage of a hydrodynamic model with a soluble tracer tracking model. FVCOM, like many other hydrodynamic frameworks used for the Great Lakes (e. g., References [41] [42] [43] [44] [45] ), is a three-dimensional (3D) free-surface, primitive-equation model that solves the momentum, continuity, temperature, salinity, and density equations, and is closed physically and mathematically using the 2. 5 level turbulence submodel [46] for vertical mixing and the Smagorinsky formulation [47] for horizontal diffusion. FVCOM is solved numerically using the finite volume method in the integral form of the primitive equations over an unstructured triangular grid mesh and vertical sigma layers. This approach combines the advantage of an unstructured grid for shoreline fitting and the flexibility of local mesh refinements (similar to finite element methods), as well as numerical efficiency and code simplicity (similar to finite difference methods). These characteristics of FVCOM make it particularly attractive in applications to coastal waters [33], where strong nonlinear interactions occur due to local geometric complexity and highly variable flow patterns [4, 48]. These characteristics make the model well suited for the study of Lake Ontario described here. 

Although this study focuses on the Pickering → Ajax → Whitby nearshore (Figure 2 ), FVCOM is configured to simulate physical dynamics for all of Lake Ontario in a manner similar to that applied by Xue et al. [14], thus providing reliable representation of large-scale background circulation and the role of remote (offshore) forcing in driving nearshore water movement; additionally, this configuration avoids the impact of setting an artificial numerical boundary condition for our target region. The advantage of an unstructured grid is that model resolution varies from 2-4 km (coarse) in the open lake to 50-150 m (fine) in the targeted nearshore, affording a high degree of resolution across the 25-km study site and adequately resolving the geographic complexity and coastal hydrodynamic conditions of that system (Figure 3 ). The model configuration yields a total of 41,442 grid elements (cells) in the horizontal plane. Vertically, 20 evenly spaced sigma coordinate layers are implemented with a resolution of <1 m in the nearshore region of interest and 5-10 m in the open lake. The bathymetric data used in the hydrodynamic model are derived from a three-arc-second (~90-m cell size), grid-based dataset from https://www. ngdc. noaa. gov/mgg/greatlakes/greatlakes. html. The Great Lakes bathymetric data were compiled by the National Oceanographic and Atmospheric Administration (NOAA) and made available to the public for Great Lakes science. method in the integral form of the primitive equations over an unstructured triangular grid mesh and vertical sigma layers. This approach combines the advantage of an unstructured grid for shoreline fitting and the flexibility of local mesh refinements (similar to finite element methods), as well as numerical efficiency and code simplicity (similar to finite difference methods). These characteristics of FVCOM make it particularly attractive in applications to coastal waters [33], where strong nonlinear interactions occur due to local geometric complexity and highly variable flow patterns [4, 48]. These characteristics make the model well suited for the study of Lake Ontario described here. Although this study focuses on the Pickering → Ajax → Whitby nearshore (Figure 2 ), FVCOM is configured to simulate physical dynamics for all of Lake Ontario in a manner similar to that applied by Xue et al. [14], thus providing reliable representation of large-scale background circulation and the role of remote (offshore) forcing in driving nearshore water movement; additionally, this configuration avoids the impact of setting an artificial numerical boundary condition for our target region. The advantage of an unstructured grid is that model resolution varies from 2-4 km (coarse) in the open lake to 50-150 m (fine) in the targeted nearshore, affording a high degree of resolution across the 25-km study site and adequately resolving the geographic complexity and coastal hydrodynamic conditions of that system (Figure 3 ). The model configuration yields a total of 41,442 grid elements (cells) in the horizontal plane. Vertically, 20 evenly spaced sigma coordinate layers are implemented with a resolution of <1 m in the nearshore region of interest and 5-10 m in the open FVCOM requires input of environmental (meteorological) forcing conditions including downward solar radiation, wind speed and direction, air temperature, atmospheric pressure, relative humidity, and cloud cover. Forcing conditions were obtained as hourly data from the National Centers for Environmental Prediction (NCEP)'s Climate Forecast System Reanalysis (CFSR), a resource that provides meteorological data beginning in 1979, with enhanced horizontal resolution (0. 2 degrees, ~20 km) available after 2011. A recent study [37] showed that CFSR reliably represents the wind patterns and radiation fields required to support robust simulations of circulation patterns and thermal structure in the Great Lakes. The CFSR wind data are also compared to local meteorological measurements from NDBC Buoy 45159 (Figure 3 ) to ensure accuracy of that input forcing condition at a local scale.
section_name: Hydrodynamic Modeling Methods

section_content: A soluble tracer-based model is used to simulate the distribution of the effluent plume. The transport and mixing of the plume with a concentration field (C) is calculated using the mass conservation equation as follows:

where H is the total water depth, and u, v, and w are the horizontal and vertical components of the water velocity in the sigma (σ) coordinate system [32]. The σ coordinate transformation is commonly used in the vertical component in order to obtain a smooth representation of irregular variable bottom topography. K h is the vertical thermal diffusion coefficient, F c is the horizontal diffusion term, and C source and C sink represent sources and sinks of C. For the conservative tracer modeling used in this study, no sink/source term is included. In our model, the tracer model is embedded within the hydrodynamic model so that we can conduct "online" tracking, meaning that the tracer model can use instantaneous results of water flow and mixing from every eight-second time step of the hydrodynamic simulation to model the distribution of the effluent plume. This approach provides more accurate results but requires more computational power than "offline" tracking [49], which reads archived hydrodynamic model output (often in an hourly time interval) as input to drive the tracer model. A conservative tracer is required in simulating plume hydrodynamics so that mass transport and kinetic influences on plume distribution may be separated; conductivity is traditionally applied for this purpose. Here, we utilize nitrate (NO 3 -N) as the tracer due to its degree of enrichment proximate to WWTP outfalls relative to ambient lake water (42% and 23% in July and September 2014, respectively, for nitrate and 5% and 2% in July and September 2014, respectively, for conductivity; Figure 4 ). Additionally, NO 3 -N may be considered as approximately conservative over the time interval of measurement used here to characterize plume dimensions. Seasonal surface water nitrate depletion due to phytoplankton nutrient uptake is regularly observed in the Great Lakes [18, [50] [51] [52]. The rate of nitrate depletion at the study site, determined from contemporary measurements made on four dates in 2014 (5 June, 6 August, 3 September, and 17 September) by the Toronto and Region Conservation Authority, http://arcgis01. trca. on. ca/nswq/,), was 1. 2 µg NO 3 -N•L -1 •day -1. This rate is comparable to that observed in our measurements (1. 6 µg NO 3 -N•L -1 •d -1 and is <1% of the range in nitrate levels measured in our single day sampling events. The near-conservative nature of nitrate is also evidenced by the linearity of its relationship with paired measurements of conductivity (Figure 4 ). 

Both field observations and modeled NO 3 -N concentrations help us to better define tributary and WWTP plume dimensions. Tributary nitrate concentrations and temperatures were obtained from the Toronto and Region Conservation Authority (TRCA, trca. on. ca). Mean daily flow for the Rouge River and Duffins Creek were obtained from the Environment Canada Water Office (wateroffice. ec. gc. ca), and the flow of Carruthers Creek was estimated by applying a watershed area ratio obtained from the TRCA to the flow data for Duffins Creek. Tributary flows were scaled up from the point of measurement to represent entire watersheds. Time series of mean WWTP effluent flow, temperature, and NO 3 -N concentration were obtained from annual wastewater treatment plant reports published by the City of Toronto (web. toronto. ca) and through special requests to the Region of Durham. The Pickering Nuclear Generating Station was included as a point source of effluent (spent cooling water) with data provided by the Region of Durham, ON. The NO 3 -N concentration in the two PNGS discharges (PNGS A and PNGS B) was assumed to be equivalent to the mean lake surface (1 m) concentration measured at the two sampling locations nearest to the PNGS outfall. 

For the modeling application presented here, NO 3 -N baseline concentrations of 355 and 270 µg•L -1 were applied lakewide (for July and September, respectively). These values were chosen based on concentrations representing the lower fifth percentile of our datasets, i. e., those uninfluenced by tributary and point source discharges. It was recognized that higher constituent concentrations could be present near the western boundaries of the study system, i. e., those potentially impacted by discharges to the west within the GTA. Inspection of field measurements revealed that antecedent meteorological and hydrodynamic conditions (longshore transport) resulted in nitrate levels to the west of our study site well in excess of the baseline condition. This was accommodated by including the NO 3 -N discharge of the Ashbridges Bay WWTP (west of the study site; Figure 3 inset) in our model configuration. The plant's flow rate and NO 3 -N concentration were obtained from the plant's 2016 Annual Report (published by the City of Toronto; Tables 1 and 2 ). 

flow, temperature, and NO3-N concentration were obtained from annual wastewater treatment plant reports published by the City of Toronto (web. toronto. ca) and through special requests to the Region of Durham. The Pickering Nuclear Generating Station was included as a point source of effluent (spent cooling water) with data provided by the Region of Durham, ON. The NO3-N concentration in the two PNGS discharges (PNGS A and PNGS B) was assumed to be equivalent to the mean lake surface (1 m) concentration measured at the two sampling locations nearest to the PNGS outfall. 

For the modeling application presented here, NO3-N baseline concentrations of 355 and 270 μg•L -1 were applied lakewide (for July and September, respectively). These values were chosen based on concentrations representing the lower fifth percentile of our datasets, i. e., those uninfluenced by tributary and point source discharges. It was recognized that higher constituent concentrations could be present near the western boundaries of the study system, i. e., those potentially impacted by discharges to the west within the GTA. Inspection of field measurements revealed that antecedent meteorological and hydrodynamic conditions (longshore transport) resulted in nitrate levels to the west of our study site well in excess of the baseline condition. This was accommodated by including the NO3-N discharge of the Ashbridges Bay WWTP (west of the study site; Figure 3 inset) in our model configuration. The plant's flow rate and NO3-N concentration were obtained from the plant's 2016 Annual Report (published by the City of Toronto; Tables 1 and 2 ).
section_name: Soluble Tracer Sub-Model

section_content: FVCOM is applied in conducting three numerical experiments, corresponding to the steps outlined in Objectives and Approach. While the FVCOM/tracer-submodel framework remains constant in all cases, specifics of model application vary and, thus, are described here for each case.
section_name: Design of Numerical Experiments

section_name: Testing Model Performance by Comparison with Field Measurements

section_content: The design for numerical experiments here includes a run (run#1) for the full project period (210 days; 1 March to 30 September 2014) utilizing all tributary and WWTP loads and applied to calibration/ confirmation of meteorological data, hydrodynamic results, and tracer concentrations. Simulations were initialized for March 2014 with an initial water temperature of 2 • C (similar to Xue et al. [37] ) and run for a spin-up period through June, when summer stratification typically develops [53]. In testing model performance, meteorological conditions and the results of hydrodynamic simulations are validated with available data from buoy and ADCP measurements, and tracer model results are evaluated against data from two sampling events occurring on 19 July and 9 September 2014 (Section 2). Nitrate loads for tributary and point source inputs were calculated as the product of their individual flow and nitrate concentration (Table 2 ). As conditions of flow and concentration are relatively stable at wastewater treatment plants, we utilized mean monthly flow and mean monthly concentration. For tributaries, where conditions are more variable, we used mean daily flows. Sampling of tributary concentrations (Toronto Region Conservation Authority) was conducted on only two dates (one in July and one in September); thus, single measurements were utilized there. Model inputs and environmental forcing conditions are summarized in Table 2. FVCOM is applied here with the linked tracer model to define the plume originating from the three tributary and three point source inputs discharging directly to the study system (Tables 1 and 2, and Figure 2 ). A generic, conservative tracer is employed so that the impacts of source hydrology (rate of discharge) and receiving water hydrodynamics may be isolated from those associated with biokinetic processes. Initial and boundary conditions assume that no tracer is present in the receiving water. Tracer loads are calculated as the product of the flow for each source and the tracer concentration. In this section, concentration can be normalized equal to one, so that model output is presented as a percentage of the tracer concentration at discharge, an expression that describes the role of physical processes (hydrology and hydrodynamics) in a general fashion subsequently applicable to a broad range of POCs. This mode of expression also facilitates both bounding of the plume and characterization of tracer gradients within the plume. Simulations are performed for stand-alone and multiple sources: the Duffin Creek WWTP to examine variability in the plant's plume position (run#2: 13 days, 9-21 July); Duffins Creek to compare that tributary's footprint during high-and low-flow periods (run#3: high flow, 30 days, 21 March-21 April; run#4: low flow, 30 days, 30 June-30 July) and to examine the nearshore flushing rate (run#5: flushing, 20 days, 10-30 April) and all three WWTPs inputs to quantify plume overlap and the system-wide extent of the combined plumes (run#6: 13 days, 9-21 July).
section_name: Testing Model Performance by Comparison with Field Measurements

section_content: Finally, the model is applied to define the POC footprint of the discharge, i. e., the portion of the plume occupied in meeting the water quality standard, and the degree to which the guidelines for a regulatory mixing zone are met. The analysis is limited to a single discharge, that of the Duffin Creek WWTP, and is applied for three POCs (ammonia nitrogen, total phosphorus, and soluble reactive phosphorus). The resulting footprints are then examined with respect to regulatory guidelines to keep the mixing zone as small as possible and avoid interference with beneficial uses.
section_name: Evaluating POC Footprint Dimensions in a Regulatory Context

section_name: Results and Discussion

section_name: Testing Model Performance

section_content: Wind-driven flow plays an important role in governing mass transport in coastal regions. Although the CFSR meteorological forcing used here is considered suitable for lake-wide simulations [37], our first step for this study was to ensure that wind conditions measured at a local buoy were consistent with CFSR output for the coastal region of interest. A comparison of eastward and northward wind components at the study site (Figure 5, upper panel) demonstrates that CFSR-derived data are consistent with buoy-measured winds. Notice that buoy winds were measured at the 5-m level, which was adjusted to wind at the 10-m level based on the formula of log wind profile [54, 55],

where the reference wind speed u(z 1 ) is measured at height z 1 = 5 m, u(z 2 ) is the wind speed at height z 2 = 10 m, z 0 is the roughness length (z 0 = 0. 0002 m for the water surface), and d is the zero-plane displacement (the height in meters above the ground at which zero wind speed is achieved, assumed to be zero for the water surface). Wind and currents in the nearshore of northern Lake Ontario are well correlated. For example, a comparison of the longshore wind component derived from buoy measurements and the alongshore currents measured by the ADCP yielded a correlation coefficient of 0. 61, even though the buoy and the ADCP were not at exactly the same locations (the buoy is located farther offshore; Figure 3 ). In addition, the longshore wind component obtained from the CFSR data and alongshore currents simulated by FVCOM were well correlated (r = 0. 76); here, both results are taken at the ADCP location. These results suggest that wind patterns significantly affect flow conditions in this region, and that hydrodynamic models will capture realistic flow conditions if quality meteorological forcing is used to drive the model. This conclusion is confirmed by a comparison of flow conditions measured by the ADCP and simulated by FVCOM (Figure 6 ). Model-simulated current conditions agree well with ADCP measurements for both alongshore and cross-shore flow components (Figure 6 ). The In July, the prevailing wind (>50% of the time) was from the southwest (a longshore direction) with a daily averaged speed ranging from 2 to 6 m•s -1 (Figure 5, lower panel). Northeasterly winds (also longshore) occurred from time to time (10%), while cross-shore winds were observed less frequently. Southwesterly winds remained predominant in September; however, wind direction was generally more random at that time. The wind was generally stronger during September with a typical speed ranging from 4 to 8 m•s -1 (Figure 5, lower panel). 

Wind and currents in the nearshore of northern Lake Ontario are well correlated. For example, a comparison of the longshore wind component derived from buoy measurements and the alongshore currents measured by the ADCP yielded a correlation coefficient of 0. 61, even though the buoy and the ADCP were not at exactly the same locations (the buoy is located farther offshore; Figure 3 ). In addition, the longshore wind component obtained from the CFSR data and alongshore currents simulated by FVCOM were well correlated (r = 0. 76); here, both results are taken at the ADCP location. These results suggest that wind patterns significantly affect flow conditions in this region, and that hydrodynamic models will capture realistic flow conditions if quality meteorological forcing is used to drive the model. This conclusion is confirmed by a comparison of flow conditions measured by the ADCP and simulated by FVCOM (Figure 6 ). Model-simulated current conditions agree well with ADCP measurements for both alongshore and cross-shore flow components (Figure 6 ). The longshore component varied between -10 cm•s -1 and 10 cm•s -1 and was much stronger than the cross-shore flow (<4 cm•s -1 ) at the near-surface layer (2 m). Synoptic events causing a strong longshore flow were well captured by the model (e. g., 28 July, 12-14 August, 25 August, and 5 September 2014; Figure 6 ). Longshore flow at the near-surface layer (2 m) shows the largest variability and decreases from surface to bottom due to bottom friction. Cross-shore flows were weak throughout the water column and no significant vertical variability was observed. This result suggests that transport in the longshore direction was much more significant than for the cross-shore. 

As thermal conditions also influence flow patterns, we compared model-simulated lake surface temperature (LST) with buoy measurements (Figure 6 ). Both model results and buoy data show that the LST increased from 9 • C to 20 • C during June. Two cooling events, causing an abrupt and significant drop in temperature in July and August, were captured as well. Throughout September, both observations and model results show that the LST decreases gradually from 20 • C to 15 • C, reflecting typical seasonal cooling. This provides us with confidence in model performance at the local scale of the region of interest, which is critical to the determination of pollutant transport and characterization of plume dynamics. Finally, we note that a similar hydrodynamic model configuration was applied for lake-wide simulations of the Great Lakes from 2002-2014; model results validated at the basin scale are presented in Xue et al. [10], targeting a coupled lake-atmosphere simulation. 

J. Mar. Sci. Eng. 2019, 2, x 13 of 27 longshore component varied between -10 cm•s -1 and 10 cm•s -1 and was much stronger than the crossshore flow (<4 cm•s -1 ) at the near-surface layer (2 m). Synoptic events causing a strong longshore flow were well captured by the model (e. g., 28 July, 12-14 August, 25 August, and 5 September 2014; Figure 6 ). Longshore flow at the near-surface layer (2 m) shows the largest variability and decreases from surface to bottom due to bottom friction. Cross-shore flows were weak throughout the water column and no significant vertical variability was observed. This result suggests that transport in the longshore direction was much more significant than for the cross-shore. As thermal conditions also influence flow patterns, we compared model-simulated lake surface temperature (LST) with buoy measurements (Figure 6 ). Both model results and buoy data show that the LST increased from 9 °C to 20 °C during June. Two cooling events, causing an abrupt and significant drop in temperature in July and August, were captured as well. Throughout September, both observations and model results show that the LST decreases gradually from 20 °C to 15 °C, reflecting typical seasonal cooling. This provides us with confidence in model performance at the local scale of the region of interest, which is critical to the determination of pollutant transport and characterization of plume dynamics. Finally, we note that a similar hydrodynamic model configuration was applied for lake-wide simulations of the Great Lakes from 2002-2014; model results validated at the basin scale are presented in Xue et al. [10], targeting a coupled lakeatmosphere simulation.
section_name: Simulating Meteorological and Hydrodynamic Conditions

section_name: Simulating Tracer Plume Dimensions

section_content: Having demonstrated the efficacy of FVCOM in representing coastal hydrodynamic conditions, attention now turns to the ability of the model to simulate plume dimensions across the study site. Here, we link FVCOM with a soluble tracer model to track the plume associated with tributary and point source discharges. Model performance was first assessed by comparing measured and modeled NO 3 -N concentrations (Figure 7 ) for each sampling transect (see Figure 2 ). The transect-average concentrations presented here are calculated as the mean of the station-average concentrations making up that transect. Station averages are the means of the vertical profiles measured for each station. The model-data comparison shows that, for the July sampling, the model successfully predicts elevations in NO 3 -N levels proximate to the Highland Creek (>450 µg NO 3 -N•L -1 ) and Duffin Creek WWTP (~500 µg NO 3 -N•L -1 ) outfalls and constant, near-baseline levels to the east (~350 µg NO 3 -N•L -1 ) toward the Corbett Creek WWTP outfall. On this occasion, the discharge from the Duffin Creek WWTP is seen to impact a longshore distance of 5 km (reaching 2 km east and 3 km west of the outfall) as evidenced by lower levels of NO 3 -N between adjacent areas of WWTP influence (troughs in Figure 7 ). A similar, yet less clear pattern is seen in results from the September survey, where much less variability in NO 3 -N levels is observed across the study site. The model results agree well with observations at sampling locations to the west of the Duffin Creek WWTP outfall except at the westernmost transect, where the model under-predicts the NO 3 -N concentration by ~30 µg•L -1 (<10%) in both July and September. To the east of the Duffin Creek WWTP outfall, model results agree well with observations of near-baseline NO 3 -N levels with some underestimation at a single transect 4 km east of the Duffin Creek WWTP outfall in September. An objective comparison of modeled and measured transect mean NO 3 -N concentrations indicates a good model fit for both the calibration (r 2 = 0. 87; normalized root-mean-square error (RMSE) = 4%) and confirmation (r 2 = 0. 70; normalized RMSE = 5%) cases (Figure 7, lower panel). 

Creek WWTP (~500 μg NO3-N•L ) outfalls and constant, near-baseline levels to the east (~350 μg NO3-N•L -1 ) toward the Corbett Creek WWTP outfall. On this occasion, the discharge from the Duffin Creek WWTP is seen to impact a longshore distance of 5 km (reaching 2 km east and 3 km west of the outfall) as evidenced by lower levels of NO3-N between adjacent areas of WWTP influence (troughs in Figure 7 ). A similar, yet less clear pattern is seen in results from the September survey, where much less variability in NO3-N levels is observed across the study site. The model results agree well with observations at sampling locations to the west of the Duffin Creek WWTP outfall except at the westernmost transect, where the model under-predicts the NO3-N concentration by ~30 μg•L -1 (<10%) in both July and September. To the east of the Duffin Creek WWTP outfall, model results agree well with observations of near-baseline NO3-N levels with some underestimation at a single transect 4 km east of the Duffin Creek WWTP outfall in September. An objective comparison of modeled and measured transect mean NO3-N concentrations indicates a good model fit for both the calibration (r 2 = 0. 87; normalized root-mean-square error (RMSE) = 4%) and confirmation (r 2 = 0. 70; normalized RMSE = 5%) cases (Figure 7, lower panel). The transect-average measurements presented above serve well in defining the alongshore extent of point source effluent plumes and the juxtaposition of regions of impact with adjoining regions having constituent concentrations characteristic of the offshore. Further insight regarding plume distributions may be gained through a two-dimensional (2D) spatial representation (Figure 8 ). It is challenging to capture the details of transient patterns at every station. However, model output and data agree well in defining effluent patterns in both the alongshore and onshore-offshore The transect-average measurements presented above serve well in defining the alongshore extent of point source effluent plumes and the juxtaposition of regions of impact with adjoining regions having constituent concentrations characteristic of the offshore. Further insight regarding plume distributions may be gained through a two-dimensional (2D) spatial representation (Figure 8 ). It is challenging to capture the details of transient patterns at every station. However, model output and data agree well in defining effluent patterns in both the alongshore and onshore-offshore dimensions. On the dates of simulation, the Duffin Creek and Highland Creek WWTP plumes extend to the west and are gradually attenuated. In the case of the July survey, the model successfully predicts the strong gradient in effluent NO 3 -N concentrations near the Duffin Creek WWTP outfall. For a plume represented by NO 3 -N concentrations >480 µg•L -1, the effluent is seen to spread ~2 km to the west and in the offshore direction; concentrations immediately to the east of the outfall are much lower, however, resulting in a strong gradient near the point of discharge. The plume may be seen to extend up to ~4 km to the west if defined by a nitrate concentration >435 µg NO 3 -N•L -1. The plume extent was slightly less for the Highland Creek WWTP discharge, with a longshore influence of ~1 km based on plume concentrations >480 µg•L -1 and 4 km for concentrations >435 µg NO 3 -N•L -1. With loads much lower than those of the other study site WWTPs (Table 2 ), the Corbett Creek WWTP has a much smaller effect on tracer concentrations across the adjoining nearshore. The Corbett Creek WWTP plume is barely evident in the July survey and of modest dimensions in September (Figure 8 ), following the shoreline to the north and west. 

We emphasize that model performance in matching transient data on a point-by-point basis is particularly challenging, as the spatiotemporal variability of effluent constituent concentrations is extremely high in the vicinity of WWTP discharges. For example, plume concentrations proximate to the Duffin Creek WWTP outfall could vary from more than 500 µg NO 3 -N•L -1 to levels approaching the lake baseline over a 36-h period (see Figure S1, Supplementary Materials) depending on conditions of advective transport (currents). The agreement of model output and field observations demonstrates the ability of the linked FVCOM-tracer model to predict the extent of plume influence. Specifically, the simulations effectively locate and bound the plumes and characterize the transition regions between discharges, i. e., high concentrations becoming moderate between the Duffin Creek and Highland Creek WWTP outfalls and moderate concentrations becoming low between the Duffin Creek and Corbett Creek WWTP discharges. With loads much lower than those of the other study site WWTPs (Table 2 ), the Corbett Creek WWTP has a much smaller effect on tracer concentrations across the adjoining nearshore. The Corbett Creek WWTP plume is barely evident in the July survey and of modest dimensions in September (Figure 8 ), following the shoreline to the north and west. We emphasize that model performance in matching transient data on a point-by-point basis is particularly challenging, as the spatiotemporal variability of effluent constituent concentrations is extremely high in the vicinity of WWTP discharges. For example, plume concentrations proximate to the Duffin Creek WWTP outfall could vary from more than 500 μg NO3-N•L -1 to levels approaching the lake baseline over a 36-hour period (see Figure S1, Supplementary Materials) depending on conditions of advective transport (currents). The agreement of model output and field observations
section_name: Simulating Tracer Plume Dimensions

section_content: Characterization of tributary and point source plumes is complicated by the fact these discharges to the nearshore are continuously present and highly dynamic with respect to position: horizontally as driven by constantly shifting coastal currents and vertically as discharges seek a position of density equilibrium in a seasonally warming and cooling water column. Monitoring and modeling results presented to this point (Figures 7 and 8 ) represent a plume as only a snapshot in time. Over a time frame of days, tributary and point source plumes distribute their constituents across the nearshore in a manner that may be likened to a stationary source (gas stack) discharging to an atmosphere with variable wind speed and direction. The dynamic nature of plumes may be illustrated by performing simulations over a period longer than the snapshot (here, three simulation results within seven days are selected to illustrate strong variability in plume position), thus incorporating variability in current speed and direction. An FVCOM simulation of this phenomenon is performed using the Duffin Creek WWTP outflow for July 2014 (Table 2 ). 

Simulation results illustrate the dynamic behavior of the Duffin Creek WWTP discharge upon entry to the lake (Figure 9 ), with plume location reversing direction in response to changes in current structure. Over the period of simulation, the WWTP plume tracks east → west on 12 July, is centrally located on 15 July, and tracks west → east on 19 July; Figure 9 ). This dynamic behavior applies to both tributary and point source discharges and persists throughout the year. Thus, from an ecosystem management perspective, quantification of plume dimensions (mixing zones in a regulatory context) must accommodate this variability, with the resulting plume describing the entire region influenced by the discharge. J. Mar. Sci. Eng. 2019, 2, x
section_name: Characterizing Dynamic Plume Behavior

section_content: The annual total phosphorus loads for the three tributaries discharging to the study site, the Rouge River, Duffins Creek, and Carruthers Creek, rank 10th, 16th, and 51st, respectively, among the 95 Canadian tributaries to Lake Ontario tabulated for 2008 by Makarewicz et al. [19]. Annual total phosphorus loads for Duffins Creek were estimated to range from approximately 10 to 90 metric tons per year [26], with the upper bound approaching the annual load for the Duffin Creek WWTP (mean
section_name: Tributary Inputs

section_content: The annual total phosphorus loads for the three tributaries discharging to the study site, the Rouge River, Duffins Creek, and Carruthers Creek, rank 10th, 16th, and 51st, respectively, among the 95 Canadian tributaries to Lake Ontario tabulated for 2008 by Makarewicz et al. [19]. Annual total phosphorus loads for Duffins Creek were estimated to range from approximately 10 to 90 metric tons per year [26], with the upper bound approaching the annual load for the Duffin Creek WWTP (mean of 138 metric tons per year, 2012-2016 [56] ). This may be interpreted as meaning that Duffins Creek and the Duffin Creek WWTP are of comparable importance in mediating local nearshore nutrient concentrations. However, tributary loads exhibit significant intra-annual variability, reflecting seasonality in constituent concentration and, particularly, flow (Figure 10a ).
section_name: Tributary Inputs

section_content: The annual total phosphorus loads for the WWTP effluents received by the study site (Highland Creek, Duffin Creek and Corbett Creek) rank second, seventh, and ninth among the 38 Canadian Two FVCOM simulations (run#3 and#4) are performed here for Duffins Creek representing early spring (high flow) and late spring and summer (low flow) conditions. It is evident from these simulations that the seasonal decline in tributary flow (Figure 10a ) and attendant reduction in load (Table 2 ) result in plume areas (calculated for concentrations C ≥ 1% of the discharge level) that differ by an order of magnitude for the two periods: ~4800 ha for early spring, high-flow (March-April) conditions (Figure 10b ) and ~370 ha for late spring-summer, low-flow (May-July) conditions (Figure 10c ). Model runs were then continued for 20 days, reducing the tributary load to zero, to quantify the purging time of constituent concentrations from the nearshore over the high-flow → low-flow transition (run#5). The simulation demonstrates that nearshore constituent concentrations characteristic of the high-flow period (Figure 10a ) are purged (<1% remaining) from the system in 6-7 days (Figure 10d ) and are, thus, not a critical driver of the nutrient-alga (e. g., phosphorus-Cladophora) dynamic during the May-June growth season. This analysis may be extended to tributary wet-weather discharges during the growing season, where purging would limit the residence of their impact to several days. We conclude that the influence of tributaries on nutrient conditions during the interval of Cladophora growth is limited to the extreme nearshore at the point where the tributary enters the lake. It is noted, however, that discharges from larger tributaries to Lake Ontario (e. g., the Genesee River at Rochester, New York (NY)) likely maintain a significant plume (and are, thus, a driving force) even in the summer.
section_name: Point Source Inputs

section_content: The annual total phosphorus loads for the WWTP effluents received by the study site (Highland Creek, Duffin Creek and Corbett Creek) rank second, seventh, and ninth among the 38 Canadian plants discharging to Lake Ontario, as tabulated for 2008 by Makarewicz et al. [19]. Unlike tributaries, effluent discharges from WWTPs exhibit relatively constant flows over the year; however, like all discharges to the nearshore, plume position varies dynamically under the influence of coastal jets (Figure 9 ). WWTP plume dimensions and the interaction of juxtaposed effluent discharges were explored by performing a 13-day FVCOM run (run#6) using the meteorological and mass transport conditions of 9-21 July 2014. WWTP effluent loads were calculated as the product of flow (Table 2 ) and the conservative tracer concentration (assumed equal to one and identical for each input). Model output is presented as a percentage of the tracer concentration at discharge, thus representing the interplay of discharge rate and ambient mass transport in determining plume distribution (Figure 11 ). 

plants discharging to Lake Ontario, as tabulated for 2008 by Makarewicz et al. [19]. Unlike tributaries, effluent discharges from WWTPs exhibit relatively constant flows over the year; however, like all discharges to the nearshore, plume position varies dynamically under the influence of coastal jets (Figure 9 ). WWTP plume dimensions and the interaction of juxtaposed effluent discharges were explored by performing a 13-day FVCOM run (run#6) using the meteorological and mass transport conditions of 9-21 July 2014. WWTP effluent loads were calculated as the product of flow (Table 2 ) and the conservative tracer concentration (assumed equal to one and identical for each input). Model output is presented as a percentage of the tracer concentration at discharge, thus representing the interplay of discharge rate and ambient mass transport in determining plume distribution (Figure 11 ). As anticipated, the dimensions of the WWTP plumes are proportional to their loads: Duffin Creek > Highland Creek > Corbett Creek (Figure 11 ). The Duffin Creek WWTP plume, defined by C ≥ 2%, extends ~5 km to the east and ~3 km to the west of the outfall. The C ≥ 2% plume for the Highland Creek WWTP discharge extends ~4 km to the east and ~1 km to the west, a reduction in plume dimensions, compared with the Duffin Creek WWTP, consistent with the difference in plant As anticipated, the dimensions of the WWTP plumes are proportional to their loads: Duffin Creek > Highland Creek > Corbett Creek (Figure 11 ). The Duffin Creek WWTP plume, defined by C ≥ 2%, extends ~5 km to the east and ~3 km to the west of the outfall. The C ≥ 2% plume for the Highland Creek WWTP discharge extends ~4 km to the east and ~1 km to the west, a reduction in plume dimensions, compared with the Duffin Creek WWTP, consistent with the difference in plant flows. Similarly, the C ≥ 2% plume for the considerably smaller Corbett Creek WWTP extends ~1 km to both the east and west. The C ≥ 2% plume for the Duffin Creek and Highland Creek WWTPs extends across a significant portion of the Cladophora growth habitat (20 m; lakeward boundary defined by the green line in Figure 11 ). A value of C ≥ 1% is taken here to represent the outer boundary of the plume, i. e., the confluence of the plume with the offshore boundary concentration (here, 0%). Plume dimensions are much larger for the C ≥ 1% plume boundary as ambient turbulence works to further dilute the point source discharges. The margin of the Duffin Creek WWTP plume defined by C ≥ 1% extends ~10 km to the east and likely a similar distance to the west, although this westward extent is difficult to isolate in a multi-source model run due to the juxtaposition of the Duffin Creek and Highland Creek WWTP plumes (Figure 11 ). The C ≥ 1% plume for the Highland Creek discharge extends ~3 km to the west, while its eastward extent is not discernable, again due to confluence with the Duffin Creek WWTP plume. The Corbett Creek WWTP plume extends 2 km to the west and ~2 km to the east, a smaller extent, again due to its smaller discharge. Perhaps most striking is the result demonstrating that the C ≥ 1% boundary for the two largest WWTPs extends lakeward, fully across the nearshore zone to the offshore edge of Cladophora growth habitat. This, considered in the context of the "unbroken urban landscape" of Rao et al. [17], begins to bring the implications of point source management for control of nuisance growth of Cladophora into clearer focus.
section_name: Point Source Inputs

section_content: The offshore discharge of WWTP effluents in the Great Lakes is regulated through the concept of a mixing zone. In this portion of the nearshore environment, segregated from the general region protected for unlimited use, water quality criteria may be exceeded [12, 16]. The adoption of a mixing zone concept recognizes the limitations of mandating that wastes be treated to the level of water quality standards, but also that enlargement of the mixing zone not be used as an alternative to reasonable and practical treatment [16]. 

Evaluation of a discharge to the nearshore with respect to regulatory mixing zone guidelines requires that the dimensions of the effluent plume and POC-specific footprint(s) be established. Here, we describe a generic means for characterizing plume and footprint dimensions, yielding a first approximation of compliance with regulatory guidelines. As previously defined, the POC footprint is the portion of the plume in which the water column concentration exceeds the water quality standard. The footprint boundary may be determined through calculation of D, the dilution (as a percentage of the effluent concentration) required to meet a POC-specific water quality standard through ambient mixing,

where C effluent is the POC concentration at the point of discharge, and C standard is the POC-specific water quality standard. Note that the determination of D does not involve the FVCOM-tracer model, but rather performs the calculation directly based on two user-specified variables (C effluent and C standard ). 

Values of C effluent and C standard are specified and values of D calculated (Table 3 ) for three POCs common to wastewater treatment facilities such as the Duffin Creek WWTP. Two of these, total ammonia nitrogen (TAN) and total phosphorus (TP), are traditional pollutants presently regulated in a mixing zone context for point source discharges to Great Lakes waters [12, 16]. The third, soluble reactive phosphorus, is a presently unattended POC mandated for management of nuisance algal growth under the Great Lakes Water Quality Agreement. Values of C standard for TAN (500 µg N/L) and TP (20 µg P/L) are those specified by the Ministry of the Environment (MOE) [16] (Provincial Water Quality Objectives). No regulatory standard is currently established for SRP. A value of C standard for SRP (1. 5 µg P/L) was calculated using the Great Lakes Cladophora Model v2 [57] and represents the upper boundary for P-limited growth of Cladophora. This SRP concentration is consistent with estimates offered previously by Canale and Auer [58] (<2 µg SRP/L), Tomlinson et al. [57] (<1 µg SRP/L), and the Targets Task Team [59] (0. 9 µg SRP/L). Values for C effluent for TAN are those published in the WWTP's Annual Performance Report for 2014 [60] and those for TP and SRP are the means of average annual concentrations for 2012-2016 published in the WWTP's Phosphorus Reduction Action Plan [56]. Recalling that the footprint is defined as the region where the POC concentration exceeds the water quality standard, it follows that the value of D circumscribes the footprint boundary within the plume dilution field (see Figure 1 ). It, thus, remains to establish the spatial scale over which the required dilution occurs, i. e., what is the size of the POC footprint? The role of the FVCOM-tracer model is to provide a map of the plume dilution field (see Figure 1 ), where 0% dilution occurs at the point of effluent discharge. The value of D is then mapped onto the FVCOM-tracer model plume dilution field to define the three POC footprints (Figure 12 ). For a value of D → 0, little (or no) dilution is required to constrain the footprint to dimensions consistent with regulatory mixing zone guidelines. As the difference between C standard and C effluent increases, D becomes larger, more dilution is required, and the POC-specific footprint comes to occupy a larger portion of the plume (Table 3 and Figure 12 ). the regulatory mixing zone guidance to minimize interference with lost beneficial uses. The next largest footprint is that for TP (<2 km × 2 km; ~342 ha), which has an area more than 85 times that for TAN. With dimensions of ~17 km × 2. 5 km and an area of ~3708 ha, the SRP footprint is even larger (more than 900 times that for TAN) and extends significantly into the far-field (Figure 12 ). Thus, the footprint area for SRP approaches one (TP) to two (TAN) orders of magnitude greater than mixing zones presently considered to be acceptable in a regulatory context (TP and TAN [23] ). Over much of its extent, particularly in the Ajax nearshore, the footprint fully engages and essentially overlays habitat colonized by Cladophora (Figure 12 ). 3 ), circumscribe the area where in-lake POC concentrations exceed the water quality standard. The green line is the 20-m depth representing the lakeward extent of colonization by Cladophora [24] and approximating the accepted separation of nearshore and offshore waters [5, 6, 22]. Model results are from run#2.
section_name: Evaluating POC Footprint Dimensions in a Regulatory Context

section_content: In this study, we developed a linked hydrodynamic tracer model for Lake Ontario with a high degree of resolution for the region of interest in the northern coastal zone of the lake. Model results show a reasonable agreement with observed hydrodynamic conditions and satisfactorily simulate plume dimensions observed during sampling events. Using this validated model, we developed a 3 ), circumscribe the area where in-lake POC concentrations exceed the water quality standard. The green line is the 20-m depth representing the lakeward extent of colonization by Cladophora [24] and approximating the accepted separation of nearshore and offshore waters [5, 6, 22]. Model results are from run#2. 

The smallest footprint is that for TAN (<0. 2 km × 0. 2 km; ~4 ha). While requirements relating to mixing zones are developed on a case-by-case basis [16], these dimensions are likely consistent with the regulatory mixing zone guidance to minimize interference with lost beneficial uses. The next largest footprint is that for TP (<2 km × 2 km; ~342 ha), which has an area more than 85 times that for TAN. With dimensions of ~17 km × 2. 5 km and an area of ~3708 ha, the SRP footprint is even larger (more than 900 times that for TAN) and extends significantly into the far-field (Figure 12 ). Thus, the footprint area for SRP approaches one (TP) to two (TAN) orders of magnitude greater than mixing zones presently considered to be acceptable in a regulatory context (TP and TAN [23] ). Over much of its extent, particularly in the Ajax nearshore, the footprint fully engages and essentially overlays habitat colonized by Cladophora (Figure 12 ).
section_name: Conclusions: Management Insights

section_content: In this study, we developed a linked hydrodynamic tracer model for Lake Ontario with a high degree of resolution for the region of interest in the northern coastal zone of the lake. Model results show a reasonable agreement with observed hydrodynamic conditions and satisfactorily simulate plume dimensions observed during sampling events. Using this validated model, we developed a generic approach for providing estimates of selected POC footprint dimensions defined by the ratio of water quality standards to effluent discharge concentrations. The generic approach has the advantage of being broadly applicable to a variety of POCs, but is a first-approximation tool offering less accuracy than a POC-specific biokinetic model. We consider these to be estimates to a first approximation as we utilize a conservative tracer to identify the POC footprint for a given source hydrology (WWTP rate of discharge) and receiving water hydrodynamic conditions, isolated from attendant biokinetic processes. If the dimensions of the footprint do not conform to the regulatory guidance for an initial mixing zone, it indicates additional efforts (e. g., extension of the discharge further offshore or reduction of effluent POC concentrations) may be needed to reduce the footprint. 

These results, treating POCs as conservative substances to target mass transport processes, indicate that the previously unattended SRP POC will not meet the guidelines for mixing zone-based regulation with current treatment technologies. Since there is a much greater difference between C standard and C effluent for SRP than for TAN or TP (Equation (3), Table 3 ), the footprint for SRP is much larger than that for TAN and TP (Figure 12 ) and a greater reduction in C effluent (D = 99. 4%) is necessary. These findings confirm, as a first approximation, our hypothesis that a shift in management focus to the Lake Ontario nearshore will reveal a necessity to reassess the impact of POC plumes with respect to regulatory guidelines for coastal waters. To achieve compliance with more stringent standards in nearshore applications, the footprints for previously unattended and emerging POCs may require advanced treatment to maintain compliance in a regulatory context. 

The variation in footprints dimensions presented here for TAN, TP, and SRP demonstrate that management challenges may lie ahead in meeting the "as small as possible" and avoidance of "interference with beneficial uses" features of a mixing zone strategy for the nearshore [12, 16]. For example, while the TAN and TP management objectives are likely achieved, that for SRP remains challenging in the context of controlling nuisance algal growth. The transition in focus from the open lake to the nearshore will require separation of effluent plumes from areas colonized by Cladophora, an alga stimulated by elevated levels of soluble reactive phosphorus and presently growing at nuisance levels in the study area. Simulations performed here demonstrate that the present engineering design of the Duffin Creek WWTP may not be capable of accomplishing that separation. 

The present WWTP discharge is located ~1 km offshore and has a footprint with a lakeward boundary reaching 2-3 km offshore (Figure 12 ). Across much of its longshore dimension, the footprint fully engages and even overlays colonizable Cladophora habitat (Figure 12 ). The outfall extension option would require that the discharge be sited such that its footprint becomes separated from Cladophora habitat (i. e., at depths greater than ~25 m, uncolonizable by the alga due to near-bottom light availability [24] ), with a margin of safety. The required offshore distance could be established using the modeling platform applied here for determining plume and footprint dimensions. Technologies capable of supporting the second option, enhanced treatment, are presently available. Tertiary treatment at the Metropolitan Syracuse (NY) Wastewater Treatment Plant, discharging to Lake Ontario via the Seneca and Oswego Rivers, maintained those exceptional levels of SRP removal for more than a decade [61]. 

The results presented here represent the critical first step in engineering design of a physical (hydrologic, hydrodynamic) foundation for modeling effluent plume and POC footprint dimensions. The second step in design is to apply a coupled hydrodynamic-biokinetic model to accommodate the non-conservative nature of a target POC (e. g., soluble reactive phosphorus with FVCOM-GEM [62, 63] ). The role of biokinetics is well recognized in mediating plume attenuation in traditional POCs (e. g., nitrification of ammonia nitrogen, cell mortality in Escherichia coli, and settling for total phosphorus). That role is less clear for the soluble reactive phosphorus POC where the interplay of biokinetic source (e. g., SRP excretion by invasive dreissenid mussels) and sink terms (e. g., SRP uptake by phytoplankton) may be appropriately considered. Field monitoring and mathematical modeling studies conducted along the Pickering → Ajax → Whitby nearshore revealed no apparent coupling of local phytoplankton growth with individual point source discharges of phosphorus [11, 27] and found Cladophora biomass density to be weakly correlated with dreissenid abundance [27]. Additionally, Leon et al. [11] reported that their model was able to predict major water quality variables related to nuisance algae growth without inclusion of dreissenid mussels and associated phosphorus cycling. Nonetheless, the significance of biokinetic processes is of sufficient potential in relation to plume/footprint dimensions and, thus, management approaches to merit further attention. We conclude that our results represent a first step, pointing to the potential need for a step change in management practices addressing lost beneficial use as the technologies of the 1970s reach their practical design limits along this "almost unbroken urban landscape" of the Lake Ontario nearshore [17].
section_name: Conclusions: Management Insights

section_content: Acknowledgments: This is contribution 41 of the Great Lakes Research Center at Michigan Technological University and contribution 338 of the Upstate Freshwater Institute. Superior, a high-performance computing infrastructure at Michigan Technological University, was used in obtaining results presented in this publication.

section_content: Funding: This research was funded by a grant to Michigan Technological University from the Town of Ajax, ON. The sponsor had no role in study design; collection, analysis and interpretation of data; in the writing of the manuscript, or in the decision to submit the article for publication.

section_name: Supplementary Materials:

section_content: The following are available online at http://www. mdpi. com/2077-1312/7/5/129/s1, Figure S1
section_name: Supplementary Materials:

section_content: The authors declare no conflicts of interest.
section_name: Conflicts of Interest: