Episodic nearshore-offshore exchanges of hypoxic waters along the north shore of Lake Erie

We investigate the nearshore-offshore exchange of hypoxic waters during episodic coastal upwelling events in the nearshore waters of northern Lake Erie using intensive ﬁeld observations and a validated hydrodynamic and water quality model. We observe wind-induced coastal upwelling events to be the dominant nearshore physical process in the lake which are energized every 5–10 days. When the winds were predominantly blowing from the west or south-west, epilimnetic waters were transported to the offshore bringing in hypolimnetic waters with low temperature (8–10 (cid:1) C), dissolved oxygen (DO: 0– 6 mg L (cid:1) 1 ) and pH (6–7) to the nearshore zones. During these events, vertical diffusivity coefﬁcients decreased from 10 (cid:1) 2 m 2 s (cid:1) 1 to values as low as ~ 10 (cid:1) 7 m 2 s (cid:1) 1 . In late summer, the coastal upwelling events in the nearshore waters lower the near bottom DO to hypoxic levels (DO < 2 mg L (cid:1) 1 ). Lake-wide observations of DO and pH show that they are positively and linearly correlated while in the near-shore DO and pH experience spatial and temporal variability where upwelling events were developed, which were further assessed using a three-dimensional model. The model accuracy to reproduce offshore hypoxia was ﬁrst assessed on a lake-wide basis using a coarse resolution model for a ﬁve-year period (2008–2012) and in nearshore waters using a higher resolution model for 2013. We use the model results to delineate the near bottom areas experiencing hypoxia at time scales longer than 48 h. (cid:3) 2021 The Author(s). Published by Elsevier B.V. on behalf of International Association for Great Lakes Research. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/ 4.0/).


Introduction
Dissolved oxygen (DO) is a critical water quality parameter that significantly impacts the suitability of physical habitat in aquatic systems (Wetzel, 2001).All temperate lakes that stratify during summer will have some degree of DO depletion in the hypolimnion; and, in some cases, this can lead to hypoxic conditions that can negatively affect biota (Watson et al., 2016;Del Giudice et al., 2018).Understanding the physical and biochemical processes maintaining DO concentrations above a critical hypoxic level (DO ! 2 mg L À1 ) has been a management priority for decades (Hawley et al., 2006;Rao et al., 2008).DO in waters is supplied mainly from atmospheric oxygen and/or aquatic photosynthesis (Charlton, 1980;Chapra, 1997).Seasonal stratification can limit the vertical distribution of DO by preventing its transport into deeper waters (Rao et al., 2008;Fischer et al., 1979) while respiration of organic matter on the lake bottom or in the overlaying waters consume DO, causing the DO concentrations in the hypolimnion to gradually decrease throughout the summer.
Over 12.5 million people in the United States and Canada rely on Lake Erie for drinking water (US EPA, 2018).Decades of binational efforts through Great Lakes Water Quality Agreements (GLWQA, 1972(GLWQA, -2015) ) between the governments of United States and Canada have been aimed to reduce phosphorus loads to Lake Erie in order to limit the seasonal offshore hypoxia extent to ~2 Â 10 3 km 2 .Despite many successes in reducing loads, the seasonal extent of offshore hypoxia in the central basin has significantly expanded during the last decade from ~2 Â 10 3 km 2 to ~9 Â 10 3 km 2 (Zhou et al., 2013 andZhou et al., 2015).In 2012, Annex 4 (nutrient annex) of the GLWQA revisited nutrient load targets for Lake Erie, and GLWQA 2015 recommended an ensemble of models be used to evaluate the extent and ecological consequences of hypoxia.Subsequently, renewed binational numerical modeling efforts began to further explore and simulate DO in the central basin, supported with field observations in the offshore (Bocaniov et al., 2016;Rucinski et al., 2016;Scavia et al., 2016) and the nearshore (Rao et al., 2014;Rowe et al., 2015Rowe et al., , 2019;;Nürnberg et al., 2019).
In Lake Erie, seasonal stratification and water circulation in the central basin are mainly driven by changes in air temperature, solar radiation and wind, whereas large riverine inflows, i.e., Detroit and Maumee Rivers, dominate circulation in the very weakly stratified west basin (e.g., Beletsky et al., 2013).In the offshore region, hourly and daily thermal structure undulations and circulations are influenced by wind driven basin-scale physical processes such as near-inertial waves and surface seiches (Bouffard et al., 2012;Rao et al., 2008;Valipour et al., 2015a,b), compared to nearshore waters which are predominantly under the influence of coastal upwelling and downwelling events in summer (Valipour et al., 2019).Coastal upwelling events are developed within the coastal boundary layer (comprising a zone within ~15 km from the shoreline) when wind events push the surface layer away from the shoreline through Ekman transport and cooler hypolimnetic water is transported toward the nearshore to conserve mass (Beletsky et al., 1997).These events can significantly impact the mixing and transport of waters in the nearshore (Rao andMurthy, 2001a, 2001b;Rao and Schwab, 2007).
In the central basin of Lake Erie, offshore hypoxia occurs during seasonal thermal stratification (Hawley et al., 2006;Scavia et al., 2014).Observations also show the development of seasonal low DO zones in the eastern basin (e.g., Boyce et al., 1980).During seasonal stratification when the hypolimnion is hypoxic, coastal upwelling events can transport hypolimnetic waters into the nearshore waters of the central basin in response to winds from the south-west or to the western basin in response to winds from the south (Rao et al., 2014;Rowe et al., 2015;Jabbari et al., 2019).The offshore hypoxic waters in the central basin are characterized with a yellow-color and are highly corrosive due to low pH levels with elevated levels of manganese (Ruberg et al., 2008).When transported to the nearshore via episodic coastal upwelling events, these events have the potential to release soluble reactive phosphorus (SRP) from fine sediment (Matisoff et al., 2016;Nürnberg et al., 2019).They can also reach certain drinking water intakes in Lake Erie, which creates a risk for source water protection and requires careful management including additional costly drinking water pre-treatments (Ruberg et al., 2008;Jabbari et al., 2019;Rowe et al., 2019).Nearshore hypoxia can also adversely impact ecosystems and fish habitant by abrupt lowering of DO concentrations.For example, thousands of deceased fish were observed along the northern shore of central Lake Erie in 2012 after hypoxic hypolimnetic waters had reached the nearshore (Rao et al., 2014).Hypoxic waters can also impact the taste and odor of lake water (imparting a musty smell) and other components of the food chain (e.g., decreased survival of bottom mayfly nymphs), Matisoff (2002).
The main objective of this study is to investigate the occurrence of hypoxic waters in the nearshore waters of Lake Erie with a focus on the northern shoreline of the central basin where recent observations show that episodic nearshore hypoxia is frequently observed and at times may result in fish kills (Rao et al., 2014).We use field observations and a validated hydrodynamic and water quality model to examine the development of coastal upwelling events in the lake and identify near bottom areas experiencing hypoxia at time scales longer than 48 h.

Study area
Lake Erie (Fig. 1; 388 km long and 92 km maximum width) is the shallowest of the Laurentian Great Lakes and is commonly treated as having three distinct basins: the western basin (WB), central basin (CB) and eastern basin (EB) with maximum depths of 11 m, 25 m, and 64 m, respectively.Lake Erie has the shortest hydraulic residence time of ~three years compared to the other Great Lakes (Chapra, 1997), and receives nutrient loads from one-third of the Great Lakes' population (Rao et al., 2008).Most of the inflow to Lake Erie is from the Detroit River (average discharge of ~6000 m 3 s À1 ) which enters the lake from the northwest.The main outflow is the Niagara River, located at the north-east end of the lake.

Measurements
As part of dedicated funding and increased focus by Environment and Climate Change Canada on Lake Erie nutrient status, sensor deployment and high-resolution modeling was started in 2012 with an aim to study nearshore nutrient dynamics in coordination with the Ontario Ministry of the Environment, Conservation and Parks (MECPO).In 2013, ECCC deployed seven YSI 6600 sondes (Yellow Springs Instruments, accuracy of 2% measured values) at stations 67,360,365,475,493,591 and 885 (Fig. 1) to study the lake-wide variability of bottom ( 1m above the bottom) DO, pH, turbidity and specific conductance (Table 1).In addition, MECPO established nearshore stations along the north shore of the CB to assess the spatiotemporal changes in nearshore DO and temperature (Table 1).The MECPO moorings comprised temperature sensors (Tidbit by Bourne, accuracy of ± 0.21 °C) at six stations (RE1, RE2, RW1, WT1, WT2, and EL1) to measure the profile of water temperature every 30 min, DO loggers (RBR-1050 with an accuracy of 0.4% saturation) at multiple depths at three stations (RE 1, RE2 and WT2) to measure DO concentrations every 10 min, and four Acoustic Doppler Current Profilers (ADCPs, Teledyne RD narrowband with an accuracy of ± 0.01 m s À1 ) to measure current profiles (three 600 kHz and one 1200 kHz) at four stations (RE1, RE2, WT1, and EL1).All ADCPs were bottom mounted to measure hourly velocities to the surface every 30 min.Further details on the ECCC field instrumentations can be found in Valipour et al (2016;2019)) and details of the MECPO moorings can be found in Nürnberg et al. (2019).
ECCC's Great Lakes Surveillance Program (GLSP) data for the period 2008-2016 were used to validate our findings on the relationships between DO and pH.The GLSP comprises the Canadian federal long-term water quality monitoring program conducted on the Great Lakes.Water quality samples were taken at discrete depths using either a vertical array of Niskin bottles or a single Van Dorn bottle.DO and pH were measured onboard, immediately after collection, following standard operating protocols (Dove et al., 2009).For DO, samples were transferred to BOD bottles under laminar flow to avoid the introduction of oxygen to the samples.Until 2016, DO was measured using Winkler titrations with an accuracy and precision of 0.1 mg L À1 and 0.01 mg L À1 , respectively.In 2016, after testing to ensure equivalency of the protocols, DO is measured using a Mettler-Toledo benchtop optical probe with equivalent accuracy and precision.pH was measured on samples transferred to stainless steel vessels using a Mettler-Toledo benchtop meter (an Oakton 510 m prior to 2005) with an accuracy and precision of 0.01 units.

Vertical diffusivity
In offshore waters of lakes, persistent net heat fluxes into the lake during the summer warm surface layers and wind-induced mixing within the surface layers result in the development of thermal stratification.In the nearshore waters and during coastal upwelling events, there is a shift of homogeneous turbulent conditions (mixed waters) in the water column to stratified conditions as the offshore hypolimnetic waters are transported into the nearshore waters; similar shifts may also occur by downwelling events or the propagation of internal Kelvin waves.Meanwhile, there is a strong current shear above and below the thermocline which eventually contributes to the generation of local vertical turbulence.The change in the vertical turbulence can be estimated by calculating the vertical diffusivity coefficient (hereafter, vertical diffusivity) during the presence of coastal upwelling events.We followed Rao and Murthy (2001b) and calculated the vertical turbulent diffusivity using a simple empirical formula as a function of Richardson number (Ri), given as: where K o = 10 À2 m 2 s À1 is an adjustable parameter and K b is the background eddy viscosity, and Ri(z) is the Richardson number profile expressed as: where u east-west (z) and u north-south (z) are the profile of mean velocity in the east-west and north-south directions, respectively, and N 2 (z) is the square of the Brunt-Väisälä frequency profile: Here, oq(z)/oz is the vertical density gradient calculated from temperature data using the UNESCO equation of state (Fofonoff and Millard, 1983), q o = 1000 kg m À3 and q(z) is the water density profile.Because velocity and temperature field data were recorded at different depths, we linearly interpolate the temperature data onto the upward looking ADCP bins in the above analysis.
Furthermore, the temporal changes in the thermal structure of nearshore waters was calculated using the vertically integrated potential energy (IPE) over the water column as described in Antenucci and Imberger (2001): IPE= R q(z,t)gzdz, where z is the water depth, and q is the water density.IPE was calculated every 30 min over the entire water column.We also calculated wind stress per unit area using the quadratic stress law s = Cd q air |W| W where W is wind speed, q air = 1.3 kg m À3 is the air density and Cd = (0.8 + 0.065 W) Â 10 À3 for W > 1 m s À1 is the drag coefficient (Wu, 1980).
Numerical Model -The occurrence of hypoxia in central Lake Erie has been well recognized and has received much attention in previous numerical modeling studies.Although box-type and 1-D models have been applied to Lake Erie to simulate DO and hypoxia since the 1980 s (Di Toro et al., 1987;Rucinski et al., 2010), the ELCOM-CAEDYM (also known as Aquatic Ecosystem Model, AEM3D) model is a three-dimensional process-based model that has been applied in Lake Erie to assess the temporal and spatial dynamics of thermal structure and hypoxia (León et al., 2011).Originally, this model was set up and calibrated for 2002, the year with the most comprehensive set of field observations (León et al., 2011(León et al., , 2012)).The ELCOM-CAEDYM model predictions of DO and temperatures were further validated during the spring-summer (Bocaniov et al., 2016;Valipour et al., 2016) and for the fallwinter-spring conditions (Oveisy et al., 2014).The model provided good predictions according to statistical measures of fit.The results of ELCOM-CAEDYM were used to investigate the sensitivity of hypoxia to variations in external nutrient loads and for the development of nutrient load -hypoxia response curves (Bocaniov et al., 2016).The ELCOM-CAEDYM model was also integrated with statistical modeling and improved the estimates of areal hypoxic extent based on the results of monitoring activities (Bocaniov and Scavia, 2016).
We have applied our existing validated high resolution threedimensional ELCOM-CAEDYM model setup to assess the development of nearshore hypoxia in Lake Erie (Valipour et al 2015a,b;2016;2019).The bathymetry of the lake was discretized with a 600 m Â 600 m uniform horizontal grid and 50 vertical layers, with spacing every 1 m for lake depths up to 45 m, and every 5 m for 45-64 m lake depths.The 600 m horizontal resolution, representing 20% of the Rossby radius of deformation (~3km; Schwab and Beletsky, 1998;Kelley et al., 2007), is deemed adequate to capture the major physical processes within the coastal boundary layer (Valipour et al., 2019).The bathymetry grids in the nearshore were carefully set up to account for very shallow depths ± 0.5 m, using the National Geophysical Data Center bathymetric data for Lake Erie (www.ngdc.noaa.gov).
The hydrodynamic model (ELCOM) solves the unsteady Reynolds-averaged Navier-Stokes equations for heat and momentum transfer across the water surface due to wind and atmospheric thermodynamics (Hodges et al., 2000) to simulate the spatial and temporal variations of the physical and transport processes.The hydrodynamic model setup and lake-wide calibration both in the offshore and nearshore waters have been described elsewhere (León et al., 2005(León et al., , 2011;;Valipour et al., 2016Valipour et al., , 2019)), demonstrating that the model is capable of resolving the predominant basin-scale offshore and nearshore physical processes, while at times having deficiencies in reproducing velocity profiles in the nearshore waters mainly due to the impact of changes in local wind forcing in the nearshore.
Meteorological data were obtained from three ECCC MET3 buoys (Fig. 1) namely WB-MET, CB-MET, and EB-MET.However, at two stations the weathervanes failed in 2013, and the wind data for stations CB-MET and EB-MET were replaced with those from nearby ECCC's weather buoys at Sta.45132 and Sta.45142,respectively.Wind data from the National Data Buoy Center (NDBC) Sta.45005 and WB-MET were similar, and thus, only the WB-MET was used.Three meteorological forcing data sets in WB, CB and EB were used to represent the spatial variability of meteorological conditions across the lake.The five major tributaries delivering nutrient loads to the lake (i.e., Detroit, Maumee, Sandusky, Cuyahoga and Grand Rivers) and the outflow (Niagara River) were used as external loading boundary conditions.
In our CAEDYM setup, the DO dynamics are driven by the following major processes: i) exchange with the atmosphere at the air-water interface; ii) photosynthetic oxygen production by phytoplankton; iii) consumption of oxygen due to respiration by phytoplankton; iv) utilization of oxygen during abiotic or biotic oxidation of reduced elements (e.g., nitrification); and, v) utilization of oxygen at the sediment/water interface (the sediment oxygen demand: SOD).
Sediment oxygen demand rates (SOD T,DO : gO 2 m À2 day À1 ) in the model are a spatial and temporal function of temperature and DO concentration at sediment-water interface, given as (Hipsey, 2008): where SOD T, DO is the modeled sediment oxygen demand rate (g O 2 m À2 day À1 ), T and DO are the water temperature (°C) and dissolved oxygen (mg L À1 ) above the sediment layer, rSOD 20 is the maximum SOD rate at 20 °C (gO 2 m À2 day À1 ), K SOD (mg O 2 L À1 ) is the constant that determines how the SOD is mediated under hypoxic conditions (K SOD = 0.5) and h is the temperature coefficient (h = 1.05).
The rSOD 20 (hereafter SOD) value is a critical parameter in our water quality model setup.However, SOD in Lake Erie's CB is poorly known (Matisoff and Neeson, 2005); therefore, a range of SOD values between 0.6 and 2.0 g m À2 per day were assessed.We compared the results of a coarse (2 km) ELCOM-CAEDYM model (León et al, unpublished) with data collected and published estimates between 2008 and 2012 (Zhou et al., 2015), and determined that SOD = 1.4 g m À2 per day best reproduces the seasonal lake-wide changes of DO concentrations.Additionally, we validated this value using 2013 nearshore data.

Wind
The most frequent seasonal wind directions were from the south-west (~225°azimuthal from true North, clockwise) with the average and maximum wind speed of 4 m s À1 and 12.7 m s À1 , respectively.There were 13 wind events that resulted in coastal upwelling events in the northern parts of CB and EB in 2013, as described below and listed in Table 2.We consider wind speeds below 5 m s À1 to be weak in Lake Erie because no forced upwelling events were observed at any monitored location below this threshold.The value is 30% lower than a previous weak wind threshold of 7 m s À1 for Lake Erie (Loewen et al., 2007).

Nearshore water temperature
The seasonal nearshore thermal structure in Lake Erie formed in the beginning of spring on approximate day of year (day) 140, as seen by the increases in solar radiation and air temperature (Fig. 2b panel i; Electronic Supplementary Material (ESM) Fig. S1).As seasonal stratification evolved (day 160), hypolimnetic waters frequently intruded to the nearshore waters via coastal upwelling events.
Coastal upwelling events -We considered a coastal upwelling event to occur when a favorable wind (e.g., from the southwest blowing toward the northeast) generated cross-shore epilimnetic currents to the offshore due to Ekman transport, and consequently pushed the thermocline front up toward the northern coast such that the temporal change in the maximum Brunt-Väisälä frequency in the water column is N(z) !0.03 (rad s À1 ) for a duration longer than the inertial period (T i = 18 h).As listed in Table 2 and consistent with episodic wind events (above), we identified 13 locally generated coastal upwelling events during the spring and summer of 2013, both in the CB and EB (Sta.RE1 and 591, respectively).The wind for these events was mainly from the south-west, had a duration of 1.3-4.7 days (i.e., ~2-10 T i ), and reached a maximum speed of 8-12 m s À1 (Table 2, Fig. 2).One exception was a wind event on day 227 in the CB which is further examined below.
Coastal Upwelling events during weak wind -On day 227, a wind event with a maximum velocity of 9.2 m s À1 from the north-west (328.2°),pushed the thermocline to the north-east portion of the CB due to Ekman transport (as observed in 15 m of waters at Sta. EL1, Fig. 3b).Later, winds subsided on day 228.8 and thus, the relaxed upwelling front returned to the northern shore of the CB.Approximately three days later (day 232.7), an upwelling front was observed at Sta. RE1 (5 km offshore) and on day 233.7 at Sta. RE2 (1 km offshore) but not at Sta. EL1, approximately 50 km away (Fig. 1).Wind events in CB (as in the WB) were relatively calm between day 230 and 235, and there were only three short duration wind events: (1) on day 231.9 with a maximum speed of 5.3 m s À1 from south-west 221°and 17 h duration, (2) on day 232.8 with a maximum speed of 5.2 m s À1 from southwest 207°and 12 h duration, and (3) on day 234.4 with a maximum speed of 6.5 m s À1 from south-west 223°and 12 h duration.The impulsive wind intervals between events (1) and (2) was 21 h, and between events (2) and (3) was 38 h.
We further explored the presence of this upwelling front at Sta. RE1 and Sta.RE2, using the wind stress components, IPE and directional currents 5 m below the surface (Fig. 4).We removed the daily and near-inertial oscillations using a low-pass filter at > 24 h (a second order Butterworth filter).As evident in the IPE timeseries, there is a local maximum at Sta. RE1 on day 232.7, and ~24 h later at Sta. RE2 on day 233.7.In either case, the currents at 5 m below the surface were negligible when IPE was at a local maximum with no substantial change in the directional velocity.
Temporal variability of DO during coastal upwelling events -Coastal upwelling events enhanced the shear currents because the current velocity above and below the thermocline were in the opposite directions (i.e., baroclinic motion).Meanwhile, the Brunt-Väisälä frequency (N) temporarily increased due to increases in the temperature gradient which resulted in decreases in the vertical diffusivity (Equation ( 1)) by at least 2 order of magnitudes from ~10 À2 m 2 s À1 to 10 À4 m 2 s À1 or lower (Table 2).Consequently, the DO supply to the bottom in the nearshore waters from the surface was restricted while the bottom layer waters in the nearshore were replaced with the offshore hypolimnetic waters resulting in lower DO and pH conditions for low redox potential biochemical processes at the sediment water interface in the nearshore.
During the 13 upwelling events observed at Sta. RE1 and Sta.591 (Table 2), near bottom DO was reduced between 15% (from 7.7 mg L À1 to 6.5 mg L À1 ) to 94% (7.5 mg L À1 to 0.4 mg L À1 ) in the CB, and between 10% (from 6.9 mg L À1 to 6.2 mg L À1 ) to 58% (8.3 mg L À1 to 3.5 mg L À1 ) in the EB.Table 2 provides detailed information on the impact of the individual coastal upwelling events on near bottom DO.Similar declines were observed for pH (Table 2).These low DO and pH conditions at the sediment-water Details of upwelling events in 2013 and their temporal impacts before, during and after the events to water temperature ( o C), DO (mg L À1 ), and pH at 1 m above the bottom at Sta. RE1 in the central basin [top] and Sta.591 in the eastern basin .N 2 (rad s À1 ) 2 and Kz (m 2 s À1 ) are maximum values in the water column.Speed and direction are the wind conditions in the offshore waters of eastern and central basins leading to the formation of upwelling events.interface created more favorable conditions for low redox potential biochemical processes (Nürnberg et al., 2019).At Sta. RE2, we also observed that DO at 2.4 m above the bottom dropped from 7.9 mg L À1 on day 231.8 to 0.1 mg L À1 on day 232.9 (98% drop), and this lasted for approximately two days while winds were relatively very weak.Multi-basin upwelling events impact on bottom DO and pH -Fig. 5 shows the temporal variability of nearshore DO and pH across the lake.It is evident that on average the nearshore DO and pH gradually decreased during the summer, with regular replenishment due to mixing or rapid changes due to upwelling events.At all stations, DO and pH followed a relatively similar trend until day 169, and also later during the season (after day 266).In the intervening time, however, distinct and site-specific upwelling events were observed at time scales on the order of days.Because the stations were located in different parts of the lake, predominant local wind direction was the key factor impacting the declines in DO and pH.For example, at Sta. RE2, DO significantly reduced on day 217 (e.g., from 6.5 mg L À1 to 1.6 mg L À1 ), however, not at Sta. 885 or Sta.365.Likewise, on day 200, DO was significantly reduced from 9.8 mg L À1 to 1.9 mg L À1 at Sta. 882, but not at Sta. RE2.At other times, impacts were observed in the WB and CB simultaneously (e.g., day 220 when DO declined at WB Sta.493 as at CB Sta.RE 2).

Sta
Relationship between DO and pH -The tight coupling of DO and pH observed in Fig. 5 is due to the influence of photosynthesis and respiration on both variables.As photosynthesis occurs, the utilization of CO 2 produces O 2 and concurrently, pH is increased as hydroxyl ions are released (Wetzel, 2001).Conversely, DO is consumed during respiration (or hypolimnetic decay of plankton and algae) and pH declines due to the production of CO 2 and H + .Within-station relationships of pH and DO (Fig. 6a-g) further demonstrate the tight coupling of these processes.Due to the high number of sensor data collected, the observed relationships are highly significant (Table 3).
Using data collected from 2008 to 2016 for GLSP (i.e., wet chemistry measured during ship-based surveys), the positive linear relationship between DO and pH is confirmed (Fig. 6h-k).The scatter between the variables increases when additional water masses are considered.For example, measurements from spring cruises (included in Fig. 6i) show higher DO concentrations because the water column at this time of year is isothermal and more freely vertically mixed.Similarly, the inclusion of additional years of monitoring (Fig. 6j-k) results in additional scatter.
Overall, the observed linear relationship between DO and pH can be expressed as ''pH = 0.2 Â DO + 6.6" and this was found to be valid both for sensor and GLSP data.The observed relationship is even stronger for DO concentrations between 2 mg L À1 and 10 mg L À1 .There was no significant relationship between DO or pH and the other monitored parameters.For example, DO and pH were not significantly related to specific conductance or alkalinity because these variables are more reflective of the characteristics of water masses from differing origins (ESM Figs.S2 to S8).The type and quantities of the dissolved constituents present have a greater impact to these variables than do the in-lake processes of respiration and photosynthesis.
The quantity of chlorophyll a (or turbidity, as it partially represents plankton) may directly impact the amount of DO present within a particular water mass (for example, through the production of O 2 during photosynthesis by phytoplankton), but the effects of other variables (such as light and temperature) are shown to have a greater effect.In this way, chlorophyll a is higher at nearshore stations, where both temperature and light are elevated.In addition, light dependent algae would be lower in hypolimnetic waters, and the relationship with DO or pH and chlorophyll or turbidity is therefore obscured.

Model
Model Validation -The accuracy of model results was previously assessed using lake-wide offshore and nearshore observations of temperature, current profiles and DO (León et al., 2005, 2012; Valipour et al., 2015a,b;2019).Here, we provide further evidence on the model accuracy to reproduce nearshore temperature and currents using observations with a focus on the CB.As shown in Fig. 7 and ESM Figures S9 to S12, the model was able to reproduce the development of upwelling events in the nearshore waters.However, there were certain upwelling events, particularly in the CB, during which the model was not as capable in reproducing the nearshore currents, such as the events on day 182 or day 235.These shortcomings can be explained by the difference between offshore wind forcing used to reproduce the nearshore processes if land-water interactions changed the local wind patterns.For these reasons, model skill and its error in reproducing temperature and currents are quantified in ESM Figs.S9-S12.
We extensively assessed the model skill in reproducing low DO in the offshore and in nearshore waters.Consistent with León et al (2011 and2012), a recent modeling effort by Rowe et al (2019) reiterated the importance of SOD to simulate lake-wide DO.To address this issue, a range of SOD values between 0.6 g m À2 and 2.0 g m À2 were examined in a multi-year (2008-2012) simulation using an ELCOM-CAEDYM setup with 2 km grid.These multi-year simulation results on the extent of offshore hypoxia at 1 m above the bottom were compared to the multi-year probabilistic and estimated observations of Zhou et al. (2015), and the error was quantified as shown in Fig. 8 and Table 4.The results demonstrate that the model is capable of simulating the seasonal and inter-annual hypoxic extents.The sensitivity tests show that SOD value of 1.4 g m À2 produced smaller errors, and is likely most representative for Lake Erie.In addition, we used the SOD value of 1.4 ± 0.4 g m À2 per day and compared the skill of the model in reproducing nearshore lake-wide DO observations in 2013 (Fig. 9).Table 5 demonstrates that the model was able to simulate nearshore DO.Specific model deficiencies in reproducing nearshore DO could result from complications in reproducing nearshore currents in the CB, where directional currents result in the transport of hypoxic waters to the nearshore.In addition, some observations were collected at depths shallower than 1 m above the bottom while the maximum vertical resolution of the model is 1 m (see Table 1).
Model results -We use the validated model to simulate the temporal and spatial variability of low DO and pH water in the nearshore.As shown in Fig. 9, the upwelling events can develop in different regions of lake, depending on the wind conditions.In turn, and depending on the offshore DO and pH levels, the upwelling events can lower the DO and pH in the nearshore due to transport of hypolimnetic waters with low concentration levels of DO and pH into the nearshore waters.Here, we use four examples shown in Fig. 10 to illustrate the spatial variability of lake-wide DO and pH, including episodic coastal upwelling events.
As shown in Figs. 9, 10a and 10b, model results reveal that hypolimnetic waters were pushed into the southern shore on day 180.5 in response to wind with a maximum speed of 6 m s À1 from the north-east (~35°) and on day 201.5 in response to wind from the south-west at 11.4 m s À1 (~219°).Wind from the north-east generated the depth-averaged nearshore currents while the offshore depth averaged currents were in the opposite direction, characteristic of topographic waves (Mortimer, 1974(Mortimer, , 2004;;Csanady, 1975).During these events, the hypolimnetic waters of the CB had a DO > 6 mg L À1 , and thus the coastal upwelling events did not significantly decrease the nearshore DO and pH levels.Later in the season, on day 210.5 and 246.5 when the hypolimnetic waters of the CB were more hypoxic, wind from the south with speed of 9.5 m s À1 (~190°) and from the south-west with speed of 6.3 m s À1 (~280°), respectively, caused coastal upwelling events which significantly decreased the nearshore DO and pH levels in the northern CB and parts of the WB (Figs. 9, 10c and 10d).
We further use the model results to better identify the spatial distribution of zones in the lake where near bottom hypoxia can be longer than two days.As shown in Fig. 11, in the CB an area as large as ~4 Â 10 3 km 2 experienced hypoxia lasting more than a month in 2013 (between day 127 and 273), while an area up to ~12 Â 10 3 km experienced hypoxia with a duration shorter than two days.

Discussion and conclusion
We present lake-wide observations on the temporal and spatial variability of nearshore hypoxic waters during coastal upwelling events with a particular focus on the north shore of Lake Erie.
Our observational results in central Lake Erie demonstrate that coastal upwelling event is not only expected after sustained wind with duration longer than 2Ti, but it can also develop due to shorter land-water breeze cycles with a duration of ~Ti and a speed of ~5 m s À1 (e.g.day 234 in Fig. 3 as discussed further below).We quantify the vertical diffusivity associated with observed coastal upwelling events in the CB and EB as this is the key parameter describing the vertical profile of DO in the nearshore during these events.We demonstrate that DO and pH are positively and linearly correlated in Lake Erie and provide a linear regression to describe this dependency.In addition, our model results show that the extent of nearshore waters affected by coastal upwelling events strongly depends on the predominant wind speed and direction.

Episodic Nearshore-offshore exchanges of hypoxic waters during coastal upwelling events
The hydrodynamic mechanisms underpinning the generation of coastal upwelling events and their potential impact on local exchange processes have been addressed elsewhere (Rao et al., 2014;Valipour et al., 2019;Jabbari et al., 2019;Rowe et al., 2019).Due to Lake Erie's unique morphology, each basin acts independently with respect to the predominant physical processes such as near-inertial waves and coastal upwelling events (Valipour et al., 2015a(Valipour et al., , 2015b(Valipour et al., , 2015c(Valipour et al., , 2019)).While there is a persistent lake-wide inter-basin transport, at times and depending on the wind and stratification conditions, episodic inter-basin exchanges of hypolimnetic waters are also expected (Boyce et al., 1980;Bartish, 1987;Jabbari et al., 2019;Rowe et al., 2019), as reported in other lakes (e.g., Lake Geneva, Umlauf and Lemmin, 2005 or Lake Simcoe, Flood et al., 2019).For the first time, we use lake-wide near bottom observations to examine the simultaneous development of coastal upwelling along the north shore of CB and EB, and compare their impacts on the temporal variability of vertical diffusivity, DO and pH in nearshore waters.
Consistent with recent studies (Jabbari et al., 2019;Rowe et al., 2019), our observations demontrate the presence of episodic hypoxia in the WB mainly driven by wind blowing from the south or south-west, bringing hypoxic waters from the CB to the WB.Our observational and model results (Figs. 9 and 10) show that these episodic hypoxic waters can partially develop in WB and rarely develop to the west of Pelee Island.When hypoxia developed in WB in 2013, its duration was<48 h (e.g., Sta 493 and Sta 885 compared to Sta 360 and Sta 365 in Fig. 9), which is consistent with recent reports on the duration of WB hypoxia (Jabbari et al., 2019).
Coastal upwelling events can simultaneously develop in the north shore of both the CB and EB if wind conditions over these basins are relatively uniform (Table 2 and Fig. 2).However, the local impacts on vertical diffusivity, DO and pH are substantially different because of spatial differences in these parameters.For example, the CB hypolimnion tends to remain hypoxic for a longer period to the parts of the EB where low DO is typically experienced only for a short period during the summer (Boyce et al., 1980).During upwelling events, the transported hypolimnetic waters into the nearshore can significantly reduce the vertical diffusivity by at least two order of magnitude from 10 À2 m 2 s À1 to values between 10 À4 and 10 À7 m 2 s À1 .These values are comparable to the reported Kz values of ~10 À5 m 2 s À1 during upwelling events in Lake Ontario (Rao and Murthy, 2001b) and 10 À4 -10 À8 m 2 s À1 in Lake Erie's WB (Loewen et al., 2007;Jabbari et al., 2019) and 10 À4 m 2 s À1 in Lake Simcoe (Cossu et al., 2017).During the development of coastal upwelling events, the reduced Kz values limit the DO supply from the epilimnion while there are potentially slight DO transfers through the thermocline.Calculation of turbulent mixing and potential DO transfer through the thermocline require highresolution measurements of velocity, temperature or DO and are beyond the scope of the current study; however, based on the findings of Bouffard et al. (2014) in Lake Erie and Chowdhury et al. (2015) in Lake Simcoe, we estimate that turbulent mixing through the thermocline can replenish ~10% of the SOD and hypolimnion oxygen demand.
The impacts of upwelling/downwelling events on DO, nutrients and other water quality parameters and their ecological consequences have been extensively studied (Haffner et al., 1984;Rao et al., 2003Rao et al., , 2014;;Rao and Schwab, 2007;Valipour et al., 2016;Nürnberg et al., 2019).Our results show that upwelling events at a lake-wide scale (Fig. 5) can significantly reduce DO and pH levels in the nearshore, but the impacts vary spatially depending on the local, offshore concentrations of DO that are transported into the nearshore during coastal upwelling events.
Development of nearshore hypoxia during weak winds.Our observations show that winds from the south-west with speed of ~5 m s À1 and duration of 12-17 h are likely responsible for the development of a coastal upwelling event in the north of CB on day 227.This conclusion is based on our observed east-west velocity components.If the upwelling front traveled from Sta. EL1, one would expect the presence of a strong westward velocity at Sta. RE2 due to a progressive westward coastal-jet.However, if coastal upwelling events are only expected by sustained winds with a duration longer than 2 Ti, we may speculate that upwelling from observed Sta.EL1 travelled to the west as a Kelvin wave.The observed upwelling front at Sta. EL1 on day 228.8 could potentially travel ~50 km from Sta EL1 and reach to Sta.RE1 on day 232.7,suggesting a propagation speed of ~0.2 m s À1 .This observed propagation speed is comparable to the analytical propagation phase speed of Kelvin waves in a two-layer system with uniform bathymetry expressed as c=(g eH e ) 0.5 (Gill, 1982), where g is the gravita- tional acceleration, e= (q 2 q 1 ) q 2

À1
, and H e = (h 1 h 2 ) (h 1 + h 2 ) À1 , h is the layer thickness, and q is the layer density and where indices 1 and 2 refer to the epilimnion and hypolimnion, respectively.On day 232, we find the analytical ratio q 2 = 999.3kg m À3 (14 °C)q 1 = 998.4kg m À3 (19 °C) of for h 1 = 15 m (uniform depth of 25 m) in the CB which gives a theoretical phase speed of c =~0.2 m s À1 .This speed is comparable to the previous estimate of 0.2 m s À1 using observations (Rao et al., 2014) and model results in the CB (Valipour et al., 2019).Due to the shallower depth of the CB, the speed is approximately half of that reported for Lake Ontario (Csanady and Scott, 1974), and the EB of Lake Erie (c =~0.4 m s À1 ; Valipour et al., 2019).
Previous observations and model results demonstrate that after winds subside, coastal upwelling fronts in the EB travel to the west as progressive Kelvin waves (Valipour et al., 2019).Likewise, our model shows the potential for the presence of a progressive Kelvin wave along the north shore of the CB, but this is not fully supported by our limited observations.Given the model performance during this period (ESM Figs.S9 to S2), we are unable to verify the presence of Kelvin wave using the model, but it suggests the presence of a Kelvin wave during this period at Sta. EL 1 (not shown).This inconsistency might be explained by (1) insufficient field observations between Sta.EL1 and Sta.RE1 which are 50 km apart and thus no wind events for ~3 days are needed for a Kelvin wave with c = 0.2 m s À1 to travel this distance or (2) model deficiency to accurately simulate the nearshore currents during this period, as discussed earlier.As shown in Figs. 3 and 4, in CB the nearshore thermocline response to small variations in wind forcing could be immediate, and thus additional field observations may be needed between Sta.EL1 and Sta.RE2 to more accurately describe the development of upwelling fronts and their propagations in response to north-west wind events.In either case, our observa-tions for the first time show that hypoxia can be developed in the nearshore waters of CB when winds are relatively calm (~5 m s À1 ) and have a duration of [12][13][14][15][16][17].
Winds in 2013 -The direction of sustained wind events and induced alongshore wind stress are critical in developing coastal upwelling events along the north shore of Lake Erie.The wind conditions in the months of June, July and August in the year 2013, used in this study, may be different in other years.For example, as shown in ESM Figs.S13 to S18, the comparison of wind roses, wind speed, and integrated wind stress (Austin and Barth, 2002) show that in 2013 wind events were mainly from the south-west while in other years the south-west and north-east frequently occur.Similarly, Rowe et al. (2019) found 2013 winds were mainly from the south-west and that this is not always the case in other years.Multi-year simulations are needed to improve understanding of the interannual variability of wind events and their impacts on the development of episodic nearshore hypoxia via coastal upwelling events.
DO and pH -The strong and positive relationship between DO and pH is due to the effects of photosynthesis and respiration to the carbonate equilibrium in lakes (Boto and John, 1981;Zang et al., 2011).The slight diel cycle of increased pH and DO due to photosynthesis during the day and their suppression due to respiration at night can be observed in the seasonal time series (Fig. 5).The diel range is greatest during the warmer growing season when both photosynthesis and respiration occur at higher rates, and when lake stratification reduces the atmospheric replenishment of lake .The linear relationship ''pH = 0.18 DO + 6.66" that we observed is remarkably similar to that observed for a marine system ''pH = 0.18 DO + 6.57" (Boto and John, 1981), suggesting that the relationship is independent of salinity and ionic conditions and may be a general relationship that be used elsewhere.Further applications to other bodies of water would be required to confirm this.Hypoxic Duration-Brief periods of hypoxia may wipe out some types of biota and have negligible effects on others in a coastal setting.For example, there was no evidence of dead fish nor strong smell immediately after an anoxic upwelling event in the CB in 2013, yet there was a windrow of dead large zooplankton; likely causalities of that event (E.T. Howell, field observations, MECPO).In contrast, similar coastal upwelling events in 2012 resulted in thousands of dead fish and a strong smell of sul-     4.

Table 4
Multi-year assessment of model performance in reproducing hypoxic areal extent using Root Mean Square Error (RMSE) between the ELCOM-CAEDYM results presented here and the probabilistic model and observations reported in Zhou et al. (2015).fur (Rao et al., 2014).These field observations suggest that more prolonged events may have the potential to extirpate sensitive benthic biota.We may conclude that the impact of anoxic upwelling events on benthic species can vary widely in their tolerance of hypoxic duration and biotic composition, and as such it may exert a selective effect on assemblages of species over time, which may require persistent biota monitoring during seasonal upwelling events.Additionally, nearshore hypoxia development has significant consequences to pH, nutrients, and apparent redox potential at the sediment-water interface (Wetzel, 2001;Golterman, 2004;Rao et al., 2014;Matisoff et al., 2016).The duration of hypoxia plays a critical role in elevating the sequence of ensuing biochemical and microbial shifts.Consistent with Loh et al. (2013), Matisoff et al. (2016) excluded the first two days of data in their study of the release of SRP from WB sediment under anaerobic conditions.Likewise, laboratory experiments demonstrated that there is a significant difference between 12 h and 4 day in shifts of bacterial communities when sediments are exposed to low DO and pH waters (Pett-Ridge and Firestone, 2005).These studies suggest that the localized nearshore impacts of hypoxia strongly depend on the duration of these events when they transfer low DO and pH waters to the nearshore.Lack of ionic exchanges from the sediment during upwelling events are also suggested in our observations at certain times when specific conductance and DO appear to be independent (ESM Figure S2).

RMSE between ELCOM-CAEDYM model and observations as reported in
Coastal upwelling events in Lake Erie, especially with respect to the redistribution of biochemically relevant constituents, have the potential to temporarily induce similar bio-chemical conditions as those observed in coastal oceans, e.g., temporal disturbances in pH, mortality, size structure of the population, reproductive behavior, and spatial distribution of aquatic communities as discussed in Chan et al (2008), Breitburg et al. (2018), and scientific reports by the Committee on Environment and Natural Resources (2010).For example, Breitburg (1992) showed that mobile fauna will recolonize disturbed areas after the development of episodic hypoxic waters in coastal waters of Chesapeake Bay, causing a temporal shift in ecological interactions.Using a 20,000 km 2 area of coastal waters along south-western Africa, among other examples, UNESCO (2018) documented that episodic nearshoreoffshore exchanges of anoxic and hypoxic waters and production of H 2 S through coastal upwelling events can kill fish or other animals in the coastal waters.Likewise, there are many reports on the development of low pH waters in the coastal oceans through coastal upwelling events, for example the intrusion of acidic waters (aragonite saturation < 1.0 and pH < 7.75) onto the continental shelf of western North America from Queen Charlotte Sound, Canada, to San Gregorio Baja California Sur, Mexico (Feely et al., 2008).
Our results (Figs. 5,10) demonstrate that the near bottom hypoxic waters can be transported across the lake every 5-10 days depending on the wind conditions.Synoptic surveys with durations of five days or less may therefore be insufficient to capture the extent of hypoxia at any given time, and are inadequate to monitor the duration of hypoxia.The deployment of buoys and

Fig. 1 .
Fig. 1.Upper panel, map of Lake Erie showing the locations of the thirteen monitoring sites in the nearshore waters of Lake Erie in 2013 -see Table 1 for more details.The solid gray contour lines are the bathymetric depths.The green circles indicate the location of ECCC GLSP sites.The triangles indicate the six offshore weather buoys namely WB-Met, CB-Met and EB-Met (western, central and eastern metrological buoys, respectively), Sta.45005 (NDBC weather buoy), Sta.45132 and Sta.45142 (Environment and Climate Change Canada weather buoys).Lower panels, the zoomed windows of deployment sites for the dashed windows areas in the upper panel identified as (a), (b) and (c).Filled contours in (a-c) are the bathymetric contour with a legend shown in panel (b).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Observed 2013 timeseries of (a) wind speed and direction at Sta.45132 every 60 min, and depth vs. time contour plots of Sta.RE1 (b) temperature (c) east-west velocity, (d) north-south velocity, (e) Brunt-Väisälä frequency N 2 and (f) the base-10 logarithmic scale of vertical diffusivity Kz, (g) the maximum base-10 logarithmic scale of Kz in the water column.(h) wind speed and direction at Sta. 45,142 every 60 min.(i-n) are the same as (b-g) at Sta. 591.The black lines on temperature contour plots in (b) and (i) show the 16 °C isotherm, characteristic of the thermocline displacement during this period.Black arrows show the upwelling events listed in Table 2. Wind direction is azimuthal from the true North (clockwise) and time is in GMT.

Fig. 3 .
Fig. 3. Observed 2013 timeseries of (a) wind speed and direction at Sta.45132 every 60 min; and (b) top, middle and bottom panels are the depth vs. time contour plots of temperature, east-west velocity and north-south velocity at Sta. EL1.(c) and (d) are the same as (b) for Sta.RE1 and RE2, respectively.The timeseries of dissolved oxygen (DO) is at 2.4 m above the bottom in Sta.RE2, and the dashed green line indicates the hypoxic line (DO = 2 mg L À1 ).The black lines on temperature contour plots in (b), (c) and (d) show the 16 °C isotherm, characteristic of the thermocline displacement during this period.Wind direction is azimuthal from the true North (clockwise) and time is in GMT.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. Observed low-pass filtered at > 24 h timeseries of (a) 60-min wind stress components at Sta. 45132; (b) integrated potential energy (IPE) at site EL1; (c) velocity components at 5 m below the surface at site EL1 (d,e) and (f,g) are the same as (b,c) for sites RE1 and RE2, respectively.Time is in GMT.

Fig. 6 .
Fig. 6.From 2013 in Lake Erie, the scatter plots of dissolved oxygen (DO) and pH at (a) Sta.885, (b) Sta.591, (c) Sta.493, (d) Sta.475, (e) Sta.365, (f) Sta.360, and (g) Sta.67.(h) the data for the scatter plots of DO and pH obtained during ECCC GLSP cruises in the nearshore (shallower than 15 m), offshore (deeper than 15 m) and lake-wide between July and September 2013.(i-k) are the same as (h), for the year of 2013, between July and September 2008-2016, and the whole years of 2008-2016, respectively.The red lines are the linear regressions and the regression coefficients are presented in Table 3. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. Comparison between modeled and observed temperature at Sta. RE 2 for day 157 through 250 in 2013.Depth vs. time contour plots of (a) observed temperature (b) ELCOM-CAEDYM model results.(right panel) Root mean square error (RMSE) at various depths between observed and modeled temperature.Time is in GMT.

Fig. 8 .
Fig. 8. Comparisons of observed vs modeled dissolved oxygen (DO) to reproduce the maximum extents of hypoxic zone in central Lake Erie between 2008 and 2012 using a range of sediment oxygen demand (SOD) from 0.6 to 2.0 g m À2 per day.Root mean square error (RMSE) between observed and modeled extents are presented in Table4.

Fig. 9 .
Fig. 9. Comparison between observed and modeled dissolved oxygen (DO) using a range of sediment oxygen demand (SOD) of 1.4 ± 0.4 g m À2 per day for day 157 to 250, 2013 at (a) Sta.885, (b) RE2 depth of 4.1 m, (c) RE2 depth of 6.1 m, (d) RE1 depth of 13.1 m, (e) WT2 depth of 6.5 m, (f) Sta.591, (g) Sta.493, (h) Sta.475, (i) Sta.365, (j) Sta.360, and (k) Sta.67.Root mean square error (RMSE) between observed and modeled extents are presented in Table 5.In each panel, the dashed green line indicates the hypoxic line (DO = 2 mg L À1 ).Time is GMT.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .
Fig. 10.2013 on day 180.5.201.5, 210.5 and 245.5, [left] modeled temperature at 1 m above the bottom and depth-averaged water circulations shown as vector plots and [right] modeled DO and pH at 1 m above the bottom (pH was estimated using pH = 0.18 DO + 6.66) and SOD of 1.4 g m À2 per day.The white contour lines are the contour depths with 10 m intervals.

Table 1
Details of field instrumentations in Lake Erie in 2013 as shown in Fig.1.The YSI parameters comprise water temperature, DO, pH, Specific conductance and Turbidity.

Table 3
Linear regression and correlation coefficients, R, calculated from DO and pH data in Lake Erie.The regression coefficients a and b are constants in the equation ''pH = a Â DO + b".Nearshore and offshore definitions relate to station (i.e., sounding) depth.

Table 5
Assessment of ELCOM-CAEDYM performance in reproducing nearshore DO using SOD = 1.4 ± 0.4 g m À2 per day.Numbers in the table show the Root Mean Square Error (RMSE) between ELCOM-CAEDYM results and nearshore observations at sites shown in Fig.1.