The Spring Runoff Event, Thermal Bar Formation, and Cross Margin Transport in Lake Superior

J. Great Lakes Res. 30 (Supplement 1):64–81 Internat. Assoc. Great Lakes Res., 2004

The Spring Runoff Event, Thermal Bar Formation, and Cross Margin Transport in Lake Superior Martin T. Auer* and Thomas L. Gatzke Department of Civil and Environmental Engineering Michigan Technological University 1400 Townsend Drive Houghton, Michigan 49931 ABSTRACT. Retention of materials delivered during the spring runoff event in the nearshore is mediated by physical phenomena such as the thermal bar. Here we investigate the timing of the runoff event vis-à-vis formation of the thermal bar and quantify its impact on the fate of materials of terrestrial origin (TSS, total suspended solids) in the nearshore waters of Lake Superior. It was estimated that the spring runoff event delivers an average of 70% of the annual TSS load for a 33-year period of record (1966–1998). The average dates for the midpoint of the spring runoff event and thermal bar formation are 9 April (range, mid-March to late April) and 2 May (range, early April to late May), respectively. This offset of 3–4 weeks provides an opportunity for materials which are discharged and retained in suspension by wave action to be transported from the nearshore before trapping by the thermal bar can occur. Model calculations indicate that 23 ± 28% (range < 1% to > 99%) of the spring runoff event TSS loading remains available for trapping at the time of thermal bar formation. In terms of mass, predicted retention varied from negligible amounts to in excess of 700,000 metric tons. The spring runoff event is predicted to have the potential to raise average nearshore TSS levels more than 4.1 ± 2.7 times above background levels and to sustain that elevation for more than 60 days. The range in the nearshore TSS response reflects interannual variability in the magnitude of the spring runoff event and its timing relative to the formation of the thermal bar. INDEX WORDS: load.

Thermal bar, spring runoff event, Lake Superior, Great Lakes, total suspended solids

INTRODUCTION The coastal margins of large lakes and oceans are among the world’s most productive waters, supporting significant sport and commercial fisheries. Coastal waters receive seasonal discharges from nutrient-rich tributaries as well as inputs from subsurface waters during wind-driven upwelling events (Chen et al. 1997). For this reason, the coastal margin is a major gateway for the transfer of materials from terrestrial environments to large lakes and oceans and is the site of intense biological, chemical, and geological processing (Henrichs et al. 2000, Klump et al. 1995). The potential for coastal ecosystems to assimilate anthropogenic inputs and to sustain a viable and healthy fishery is not well known because the flux of materials through the region and the transformations that they undergo are not fully *Corresponding

understood (Klump et al. 1995). Thus our appreciation of the role of cross-margin transport in ecosystem function and our ability to project the response of coastal environments to anthropogenic perturbation remains limited. In recognition of this need, the federal interagency Coastal Ocean Processes (CoOP) Program has set as its goal “to obtain a new level of quantitative understanding of the processes that dominate the transports, transformations and fates of biologically, chemically, and geologically important matter on the continental margins” (Brink et al. 1992). The Laurentian Great Lakes include some of the most densely populated and heavily utilized of the coastal waters of the U.S. and, as closed basins dominated by their coastal nature, represent sites where the influence of biogeochemical processes are magnified beyond that of many marine environments (Klump et al. 1995). The KITES Project (Keweenaw Interdisci-

author. E-mail: [email protected]

64

Thermal Bar Retention of Runoff to Lake Superior plinary Transport Experiment in Superior), sponsored by NSF as part of the CoOP Program, explored physical processes mediating the cross-margin transport of materials discharged to (tributaries), introduced to (upwellings) and generated within (biogeochemical transformations) the coastal waters of Lake Superior. The research presented here focuses on the role of thermal fronts in mediating the cross margin transport of terrigenous materials delivered during the spring runoff event in that lake. Thermal fronts are characterized by strong horizontal gradients in temperature and density and are commonly observed in large lakes and coastal oceans. Ullman et al. (1998) provide an excellent description of the phenomenon and identify the essentially ubiquitous presence of fronts (both temporally and geographically) in the Laurentian Great Lakes. Tikhomirov (1963) was among the first to study the thermal bar, describing a front which forms in spring as temperatures rise more rapidly in shallow nearshore waters than in deeper offshore regions. Convergent circulation where warm (> 4°C) nearshore and cold (< 4°C) offshore water masses meet establishes a frontal boundary, the thermal bar, corresponding to the 4°C isotherm (Ullman et al. 1998). Continued warming of nearshore waters forces the frontal boundary lakeward and the bar persists as a migrating feature of the lake’s thermal regime until the onset of vertical stratification. Thermal fronts act as boundaries between water bodies, limiting the exchange of heat and mass (Huang 1972, Blanton 1986, Lathrop et al. 1990). The spring thermal bar, originally described for the Great Lakes by Rodgers (1965), has been the subject of scientific investigation for many years and is recognized as a system-wide feature of the thermal regime of all the lakes (Rodgers 1965, Bolgrien and Brooks 1992). The ability of the thermal bar to inhibit mixing and maintain separation of water masses in the Great Lakes has been clearly demonstrated (Rodgers 1966, Noble and Anderson 1968, Hubbard and Spain 1973, Spain et al. 1976). Subsequent studies in the Great Lakes and elsewhere have confirmed that waters inshore of the thermal bar become enriched in materials of terrigenous origin as tributary discharges associated with the spring runoff event become trapped shoreward of the front (Hubbard and Spain 1973; Spain et al. 1976; Moll et al. 1993a, 1993b; Likhoshway et al. 1996). Nearshoreoffshore gradients in biological materials (bacteria, Menon et al. 1971; phytoplankton, Nalewajko 1967; Stoermer 1968; Moll et al. 1993a, 1993b; Likhoshway et al. 1996; and zooplankton, Evans et al. 1980)

65

have been observed as well and are often attributed to the trapping of nutrients in nearshore waters. OBJECTIVES Although recognized as a prominent feature of Great Lakes limnology for almost 40 years, the role of the thermal bar in establishing and maintaining gradients across the coastal margin has not been evaluated in a systematic and quantitative manner. The potential for the thermal bar to trap materials of terrigenous origin is coupled to the timing of the spring runoff event, vis-à-vis formation of the front. Where the bar is in place over a significant portion of the runoff period, there is potential for capture, nearshore waters may become enriched, and offshore waters would receive nutrient re-supply through migration of the front and/or resuspension events. In cases where the spring runoff event is essentially complete at the time of thermal bar formation, the potential for trapping is much less and the delivery of terrigenous materials to offshore waters maximized. Thus, interannual variations in the timing of the spring runoff event and of thermal bar formation may have important effects on nutrient levels and ecosystem function in nearshore and offshore environments. Here, we characterize the timing and duration of the spring runoff event and the timing of thermal bar formation at a site on Lake Superior for a 33-year period of record (1966–1998). Tributary total suspended solids (TSS) loadings are calculated and serve as input to a mathematical model of the nearshore region. Retention of the spring runoff event solids loading by the thermal bar is then estimated for each of the 33 years, documenting interannual variation. This paper provides the first published quantification of the role of the thermal bar in mediating cross-margin transport in the Great Lakes. STUDY SITE The thermal regime of Lake Superior is unique among the Great Lakes, having the longest period of convective spring mixing, the shortest duration of summer stratification and the lowest maximum surface temperature (Bennett 1978). The thermal bar forms latest in Lake Superior (relative to the other Great Lakes) and temperature gradients between shallow nearshore waters and deeper offshore waters are generally maintained throughout the summer (Ullman et al. 1998). This research focuses on a 250

66

Auer and Gatzke

FIG. 1. Study area tributaries. The three named streams contribute approximately 75% of the annual hydrologic input. km stretch of coastal margin on Lake Superior’s southern shore adjacent to the Keweenaw Peninsula of Michigan and includes three major and ten minor tributaries and their watersheds (Fig. 1). These tributaries provide ~8% of the flow to the lake (mean annual discharge 1966–1989: Ontonagon River, 41 m3⋅s–1, Bad River, 30 m3⋅s–1, Sturgeon River, 26 m3 ⋅ s–1, and minor tributaries, 40 m3⋅s–1). The drainage basins represent ~9% of the Lake Superior watershed (84% of that between Odanah, Wisconsin, and Copper Harbor, Michigan) and are representative of regional land use (largely forest) and soil types (Robertson 1997). Streams in the region draining watersheds with predominantly clay soils (the Bad and Ontonagon rivers) are highly turbid, while those draining sand and gravel soils (the Sturgeon River and minor tributaries) are less so. All streams are highly colored (humic materials), although this is masked in some cases by turbidity. The lacustrine portion of the study site is taken to include the nearshore waters of Lake Superior to a distance of three kilometers from shore. The lakeward boundary is coincident with the position of the thermal bar

during its early development. Lake water depth in the study site increases from west to east, averaging 13 m near the Bad and Ontonagon rivers and 86 m at Copper Harbor. Variations in water depth across the study site lead to differences in the timing and position of the thermal bar (cf. Spain et al. 1976). METHODS Retention of runoff event solids in the nearshore at the time of thermal bar formation over the period 1966-1998 was quantified by: (1) developing a time series of TSS loads, (2) determining the date of thermal bar formation, and (3) calculating solids fate and transport in the nearshore. Methods supporting the first two steps are treated here and the third in a subsequent section devoted to modeling. TSS Loads Development of a daily time series of TSS loads requires estimates of tributary flows and TSS concentrations. Instantaneous and mean daily flow data

Thermal Bar Retention of Runoff to Lake Superior for the Ontonagon, Bad, and Sturgeon rivers were obtained for the period 1966–1998 from the USGS on-line database (http://waterdata.usgs.gov/nwis). Flows at USGS gauging stations were adjusted to include the entire drainage area by multiplying reported flows by the ratio of the gauged to ungauged watershed area. Daily TSS concentration estimates were derived from concentration vs. flow (C/Q) plots. Samples were collected by autosampler (4 times daily) on the Ontonagon and Sturgeon rivers and by grab collection (7 times) on the Bad River over the period April–October 1998 and in the spring of 1999. TSS analysis followed the procedure outlined by APHA (1998). Additional TSS data for the Bad River were obtained from the National Stream Quality Accounting Network (NASQAN, Robertson 1997). The resulting TSS database consisted of approximately 50 data points for each major tributary and 7 data points for each minor tributary. Land use for the watersheds making up the study site are predominantly forested (90–97%; Robertson 1997) and thus the C/Q relationships developed here are considered appropriate for application over the 30+ year period of record. All C/Q plots exhibited positive slopes (runoff-generated loads) and were flow-stratified yielding separate regressions for base and high flow regimes (Fig. 2). Time series of TSS loads for the Ontonagon, Bad, and Sturgeon rivers were calculated as the product of measured daily flows and estimated TSS concentration using the U.S. Army Corps of Engineers’ program FLUX (Walker 1996). It is noted that the Sturgeon River discharges to Portage Lake and then splits, flowing to the study site via the North Entry to the Keweenaw Waterway and/or to Keweenaw Bay via the South Entry to the Keweenaw Waterway. An exact quantification of the loading contribution to the study site from the Sturgeon River is problematic because of the dynamics of flow partitioning (cf. Churchill et al. 2004a) and losses of solids within the waterway. Here, it is assumed that the Sturgeon River TSS load behaves conservatively in the waterway and that all of the flow exits at the North Entry. This conservative assumption provides an upper bound for the estimate of the Sturgeon River contribution to TSS levels at the study site. The ten minor tributaries, accounting for ~29% of the system watershed area and flow, are ungauged. Development of a paired C/Q database for these streams, of the size required to support loading calculations, is beyond the scope of this study. Instead, a unit area load approach was utilized. Mean base flow TSS concentrations for the minor tributaries

67

FIG. 2. Concentration vs. flow (C/Q) relationships for total suspended solids (TSS) in the Bad, Ontonagon, and Sturgeon rivers.

were compared with those of the three gauged rivers. It was determined that the base flow TSS concentrations most closely compare with those of the Sturgeon River. The daily time-series TSS load for the Sturgeon River was thus multiplied by the ratio of the minor tributary drainage area to the Sturgeon River drainage area to yield individual daily time-series TSS loads for the minor streams.

68

Auer and Gatzke

Thermal Bar Formation The value of remote sensing as a means for tracking the timing and structure of the thermal bar has been recognized for some time (Hubbard and Spain 1973, Strong 1974, Mortimer 1988) and the approach is presently enjoying widespread application in the Great Lakes (cf. Schwab et al. 1992, Budd et al. 1998). Satellite images (NOAA TIROS-N, AVHRR) characterizing Lake Superior surface temperatures over the period 1995–1998 were used to identify the date of thermal bar formation. The thermal bar was considered to be in place when persistent warming (> 4°C) of nearshore surface waters was observed in satellite images to a distance of 3 km from shore. That distance is consistent with ground truth observations of the spatial extent of early thermal bar formation in Lake Superior during early development (Hubbard and Spain 1973). Daily temperature data from the Ontonagon, Michigan, drinking water supply intake were used to refine the precision of this estimate and to extend the analysis to years for which satellite images are not available. The intake is 0.9 km offshore and 0.5 km west of the Ontonagon River mouth. Examination of turbidity data from the water supply confirms that the intake largely tracks lake, as opposed to river conditions. RESULTS AND DISCUSSION TSS Loads Daily load time-series for TSS (Fig. 3) were calculated for each of the thirteen tributaries for the period 1966-1998. The Bad River was found to deliver the greatest percentage (34%) of the annual average load, followed by the Ontonagon (30%) and Sturgeon (13%) rivers. The 10 minor tributaries collectively accounted for 23%, with no single stream contributing more than 7% of the annual total. The Bad and Ontonagon rivers are unique among those sampled in that their watersheds are largely characterized by easily mobilized clay soils, yielding higher TSS concentrations than those tributaries draining predominantly sand and gravel substrates. Thus, the Bad and Ontonagon rivers carry TSS loads (averaging 64%) out of proportion to their contribution to the study site drainage area (55%). Annual TSS loads varied by more than a factor of ten over the 33-year period of record (Fig. 4) ranging from 0.12 × 106 M.T. in 1987 to 1.54 × 106 M.T. in 1996. The load for 1998 (0.46 × 106 M.T.), the year of C/Q sampling, was 64% of the 33-year mean (0.72 × 106 M.T.). Interannual differences in loads

FIG. 3. TSS loading time series for 1996 for the three major tributaries (Bad, Ontonagon, and Sturgeon rivers) and one minor tributary (Iron River). are closely tied, but not linearly related, to differences in flow. For example, the 1998 flow was half that of 1996, but the TSS load decreased by a factor of 3 (4.6 × 105 M.T. to 1.5 × 106 M.T.). TSS loads are influenced not only by the magnitude, but also by the structure of the flow regime. Delivery of loads within a brief (days), but intense, event tends to deliver more TSS than that occurring over an extended (weeks), but less intense, runoff period. Years with large spring runoff events (1997, 1998) delivered greater loads than those with comparably smaller spring peaks but similar annual flows (1994, 1995). In the Lake Superior region, peak flows generally occur during spring snowmelt, not during rain

Thermal Bar Retention of Runoff to Lake Superior

FIG. 4.

69

Annual TSS load summed for all major and minor tributaries.

events, and yearly snowfall has a substantial influence on spring flows. Years with heavy snowfall and rapid melting periods generate the largest and most intense TSS loads. The spring runoff event (see Fig. 3) is the most prominent feature of the TSS loading regime, with an average of 71 ± 17% of the annual load (range, 32–98%) contributed during that period. Considerable interannual variation was observed in the magnitude, timing, and duration of this phenomenon. In order to quantify this variability, event timing was characterized according to two criteria: (1) the event was considered to be initiated when Ontonagon River TSS loads rose above baseflow levels and (2) the event extended until baseflow loads were sustained for a period of 3 days. A representative date for the event was established as the day in which one-half of the total event load had been delivered, termed here the center of mass. The mean date for the center of mass was 10 April (± 12 days; range, 15 March–29 April) and the average event duration was 47 ± 15 days, (range, 17–77 days). Characteristics of the timing and duration of the spring runoff

event are summarized for the 33-year period of record in Figure 5. The daily time series of TSS loads calculated here are employed subsequently as input for a Lake Superior nearshore TSS model. This provides an excellent opportunity to explore the characteristics of the trapping phenomenon within the context of natural variation in regional hydrology. Thermal Bar Formation The date of thermal bar formation was established for each of the 33 years in the period of record. Dates were established from satellite images for the interval 1995–1998. The progressive warming of nearshore waters, from the time of cold isothermal conditions through full formation of the thermal bar, is illustrated for 1998 in Figure 6. Based on these images, the thermal bar is judged to have been formed between 9 and 18 April, with a midpoint of 14 April. Proceeding in a like manner, the date of thermal bar formation was identified for the other

70

Auer and Gatzke TABLE 1. Timing of thermal bar formation as identified through satellite images and temperature data from the Village of Ontonagon water intake.

Year 1995 1996 1997 1998

FIG. 5. Dates of the spring runoff event and thermal bar formation, 1966–1998.

years available from the satellite image data base (1995–1998; Table 1). Establishment of dates of thermal bar formation for years outside the satellite image data base required an alternate approach. An inspection of the water temperature records for the Ontonagon, Michigan, drinking water intake revealed that the midpoint date of thermal bar formation, identified from satellite images, coincided with an intake temperature of 5°C. This criterion was applied to identify the date of thermal bar formation for each year over the interval 1966–1998 (Fig. 5). Considerable variability was noted in the date of formation of the thermal bar, ranging from 6 April (1969) to 28 May (1972) with an average date of 2 May (± 13 days). To some degree, this variability is due to meteorological conditions. For example, one of the earliest dates of thermal bar formation (13 April) occurred in 1998 (Table 1, Fig. 5), an El Niño

Period of TB Date of TB Formation—Identified Formation –Identified by Inspection of by Inspection of Thermal Images Water Supply Data 26 Apr–7 May 29 Apr 12 May–4 Jun 23 May 10 May–27 May 17 May 9 Apr–18 Apr 13 Apr

year. Other effects on the thermal regime of large lakes (upwellings, Lake Superior, Budd 2004; duration of the ice-free season, Great Slave Lake, Schertzer et al. 2004) were reported for that year as well. There is also variability in the date of thermal bar formation within the study site, e.g., the bar tends to form earlier in the shallow nearshore waters off Ontonagon than it does in the deeper waters characteristic of that portion of the study site to the northeast. Thus estimates of the date of thermal bar formation based on observations at Ontonagon tend to provide an upper bound with respect to the potential for trapping. Inspection of the results presented in Figure 5 provides information regarding the relative timing of the spring runoff event (center of mass) and formation of the thermal bar. On average, thermal bar formation lags the center of mass in the spring runoff event by 3 weeks, leaving a substantial margin in time for terrigenous materials to be removed from the nearshore by mass transport and settling prior to thermal bar formation. There is considerable variation in the lag period, however, with bar formation occurring within a week of the center of mass on six occasions and lagging the event by more than a month eight times over the period of record. Thus it remains for the modeling exercise described below to quantify the mass and fraction of the spring runoff event solids loading available for trapping and to establish the level of natural variability in that amount. MODELING The significance of the thermal bar as a barrier to cross-margin transport is determined in part by the levels of TSS present at the time of bar formation. As shown in Figure 5, the spring runoff event most often occurs well before thermal bar formation. While this suggests in a qualitative way that there is

Thermal Bar Retention of Runoff to Lake Superior

71

FIG. 6. The progression of thermal bar (TB) development 1998 as evidenced in satellite images.

opportunity for material to be lost from the nearshore prior to formation of any barrier, it remains to quantify that loss. This is accomplished here through application of a mathematical model which accommodates those physical processes which affect the fate of TSS subsequent to their discharge via tributaries. The goal of the modeling ex-

ercise is to incorporate TSS loads and physical transport mechanisms in a simulation of nearshore TSS concentrations at the time of TB formation. The model utilizes the TSS load delivered during the spring runoff event to generate a concentration time series for each of 53 model cells and then determines the spring runoff event TSS mass and fraction re-

72

Auer and Gatzke

maining in the system at the time of formation of the thermal bar. The modeling analysis is applied to the 33-year period of record, representing a range of intervals between the runoff event center of mass and thermal bar formation, i.e., time for TSS to be lost from the system. This admittedly simple, coarsescale model does not permit fine scale resolution of TSS levels near river plumes and does not incorporate inter-annual and intra-annual variability in mass transport processes, however, we believe that it does facilitate a valuable first cut estimate of the significance of the relative timing of the spring runoff event and thermal bar formation. Model Framework and Solution Technique The model consists of 53 completely-mixed cells, each with a 5 km longshore length and 3 km offshore width, accommodating the 250 km study site from the mouth of the Bad River (near Odanah, Wisconsin) northeast to Copper Harbor, Michigan. Model cells have an average depth of 13 m from the Bad River to 5 Mile Point, Michigan, and an average depth of 86 m from that location to Copper Harbor. The model is developed upon the principles of mass balance and accommodates loads from the 13 study tributaries, inputs from resuspension, losses to sedimentation and exchange through mass transport (advection and diffusion): V⋅

dCi = Wi + Wr + Q ⋅ (Ci −1 − Ci ) + Elong ′ ⋅ dt

(1)

(Ci +1 − Ci ) + (Ci −1 − Ci )  + Eoff ′ ⋅ (Cbc − Ci ) − Sv ⋅ As ⋅ Ci

where: V is cell volume (m 3), C is the TSS concentration (g⋅m –3 ), W i is the tributary TSS load (g⋅d–1), Wr is the resuspension TSS load (g⋅d–1), Q is the advective flow (mean current velocity) into and out of each cell (m3⋅d–1), E9 is the bulk diffusion coefficient in the offshore (E9off) and longshore (E’long) directions (m3⋅d–1), Cbc is the background TSS concentration (g⋅m–3), Sv is the settling velocity (m⋅d–1), As is the cell surface area (m2), and i refers to the cell number. Equation 1 was applied to each of the 53 linked model cells and solved using a first order Euler integrator (Chapra and Canale 1988). Model runs extended from the start of the spring runoff event until formation of the thermal bar. In cases where the spring runoff event extended beyond the date of thermal bar formation, the model was run through the conclusion of the event. Tributary loads were input to the model only for

the period of simulation to permit isolation of effects directly associated with the event. Sensitivity analyses were performed to identify model inputs for use in the calibration process. The model was then calibrated, confirmed, and applied in determining TSS concentrations in the Lake Superior nearshore during and after the spring runoff event. Model output consisted of an annual daily time series of TSS concentrations for each of the 53 model cells for each of the 33 years in the period of record. Concentration data were used to calculate the TSS mass loaded during the event and the TSS mass present at the time of thermal bar formation (including TSS mass loaded subsequently if the event extended beyond the date of thermal bar formation). From this, the fraction of the SRE TSS load available for trapping by the thermal bar is determined. Model Inputs Of the 53 model cells, 13 receive a tributary TSS load, calculated as described above. Advection, flow into and out of each cell, was calculated as the product of the cell’s cross-sectional area (Ax, m2) and the current velocity (Cv, m⋅d-1). Estimates of current velocities ranged from 0.1–10 cm⋅s-1 (Churchill et al. 2004b), consistent with those reported for the study region by Viekman and Wimbush (1993). Mean values for the offshore and longshore diffusion coefficients (9.5 × 105 and 4.0 × 107 m2⋅d–1, respectively) were determined from BASS (Bathymetric Acoustic Survey System) deployments by Churchill et al. (2004a). The offshore diffusion coefficient was reduced by a factor of 100 when the thermal bar was in place, consistent with the range of reduction in thermocline vertical transport coefficients observed during stratified and unstratified conditions (400, Schnoor 1996; 100, Doerr et al. 1996; 20, Thomann and Mueller 1987). A boundary condition of 0.9 mgTSS⋅L–1 was applied in mass transport calculations, based on the average background TSS concentration recorded at the Village of Ontonagon drinking water intake during the spring runoff event. Settling velocities, averaging 2.4 ± 2.2 m⋅d–1 were determined using sediment traps (Urban et al. 2004). Sensitivity Analysis A sensitivity analysis was performed to examine the impact of coefficient selection on model performance and to identify coefficients to be used in the model calibration process. Values for many of the model inputs (cell geometry, initial and boundary

Thermal Bar Retention of Runoff to Lake Superior conditions, tributary loads) are well described, while others (settling velocity, current velocity, and offshore and longshore diffusion) are determined through field measurement and carry a degree of natural variability or uncertainty in their estimation. The sensitivity analysis consisted of the independent adjustment of these four variables to define their relative importance with respect to system response, i.e., the timing and magnitude of TSS peaks in the Lake Superior nearshore. The sensitivity analysis was performed for 1998, a year with a well-defined spring runoff event. Values for current velocity, settling velocity, and diffusion were adjusted within the range of their natural variability or experimental uncertainty and the response examined at two locations within the nearshore framework: the cell receiving the Ontonagon River load (response dominated by loads) and a cell 50 km down current (dominated by fate and transport mechanisms). The magnitude and temporal structure of the TSS response to the SRE load were used as indicators of model sensitivity. Longshore diffusion and current velocity have the greatest impact on the temporal structure of TSS concentrations, advancing and retarding (respectively) longshore movement of the TSS peak by 2 days for the range of coefficient values examined here (Fig. 7). Settling velocity and offshore diffusion are most important with respect to the magnitude of the TSS peak (Fig. 7). Settling velocity and current velocity were selected for adjustment in the calibration process, as their values are reasonably well supported here through field measurements (Urban et al. 2004 and Churchill et al. 2004a, respectively). Average measured values are applied for offshore and longshore diffusion. Model Calibration and Confirmation The model was calibrated for the Ontonagon River loading cell, the only location for which estimates of TSS concentrations are available. TSS calibration data were derived from turbidity measurements collected at the Village of Ontonagon drinking water intake. Turbidity values (NTU) for the period 1989–1998 were converted to TSS concentrations using a calibration curve developed during the 2000 spring runoff event (TSS = 0.8373⋅ NTU; R2 = 0.86). Model calibration was first performed for baseline conditions, i.e. periods of minimal tributary input, adjusting the value of the zero-order resuspension term (final value, 0.9 g⋅m–2 ⋅ d -1) until background TSS concentrations (~0.8g⋅

73

m –3 , Jeong 2002; 0.9 g⋅m –3 , water intake) were matched. With calibration to baseline conditions established, attention was then focused on modeling of the loading events. Again, calibration was performed on the Ontonagon loading cell and data from the Ontonagon water supply intake were used for comparison. A best fit was sought for both the magnitude of the peak TSS concentration and the timing of attenuation of that peak to the background concentration. TSS loading time-series for multiple years (1992, 1995, 1998) were used in the calibration process in order to fit the model under a variety of loading conditions, and therefore establish validity over the entire study period. Current velocity and settling velocity were used as calibration parameters, with offshore and longshore diffusion coefficients set at their mean values, 1.0 × 105 cm2⋅s–1 and 4.2 × 106 cm2⋅s–1, respectively. A satisfactory model fit (Fig. 8) was obtained using values of 5 cm⋅s–1 for current velocity and 1 m⋅d–1 for settling velocity. In general, current velocity controls the time to background and settling velocity determines concentration magnitude. Thus adjustment of these coefficients has independent, non-compensating effects and imparts some confidence in coefficient selection. The model was then confirmed by comparing data from the Ontonagon water intake to the predicted TSS concentrations for the other 7 years in the 10year calibration-confirmation data base (Fig. 8). The model performed well in predicting the timing of peaks in TSS, especially the occurrence of the ascending and descending limbs of runoff events. Model output failed to match observed spikes in TSS (Fig. 8), but the general magnitude of predicted concentrations was otherwise quite satisfactory. Although uncertainty in model inputs contributes to this situation, failure to match the magnitude of TSS peaks is thought to be largely an artifact of the model framework. Model-calculations assume complete mixing, i.e., that the TSS load is evenly distributed over the 15 × 106 m2 area of the model cell. Calibration-confirmation data are drawn from a single point within that area, the drinking water intake. The model’s ability to match TSS peaks is dependent on the extent to which that single point is representative of the larger area. When the river plume lies in close proximity to the water intake, the completely mixing assumption is violated, the field data do not appropriately represent the model cell and peaks in TSS are under-predicted. A better fit to those peak runoff conditions could be obtained with a more

74

Auer and Gatzke

FIG. 7. Model sensitivity to current velocity, settling velocity, longshore diffusion, and offshore diffusion for the cell receiving the Ontonagon River discharge and a cell 50 km downcurrent.

Thermal Bar Retention of Runoff to Lake Superior

FIG. 8. Model fit to TSS data from the Village of Ontonagon water intake for the 10-year period, 1989-1998. The years 1992, 1995, and 1998 were used in calibration; other years represent model confirmation.

75

76

Auer and Gatzke

FIG. 9. Percent of the spring runoff event TSS load available for trapping at the time of thermal bar formation: (a) distribution over the 33-year period of record and (b) as a function of the date of thermal bar formation relative to the SRE center of mass. A negative value for the thermal bar lag indicates that the front formed prior to occurrence of the center of mass.

fine-scale hydrodynamic framework, however, model inputs supporting such an effort are not presently available. Model Application The confirmed model was applied in estimating the percentage and mass of the spring runoff event TSS load remaining in the nearshore and available for trapping at the time of thermal bar formation. An average of 23.3 ± 28.5% of the event load was found to be available (Fig. 9a). As evidenced by the high standard deviation, considerable interannual variation is observed in the percent available, i.e., a minimum of < 1% in 1990 and a maximum of > 99% in

1969. These estimates represent a conservative lower bound for TSS availability because TSS lost to settling prior to thermal bar formation may subsequently be reintroduced through resuspension. Model runs were made with the settling term set equal to zero to simulate the fate of dissolved substances. Here, an average of 30.7 ± 27.4% was retained, with a minimum of 4% and a maximum of 100%. Thus, it may be concluded that, on average, one-quarter to one-third of the mass of terrigenous materials loaded during the spring runoff event are available for trapping by the thermal bar. Substantial interannual variation in the availability of TSS for trapping is predicted due to differences in the timing of the spring runoff event relative to the timing of thermal bar formation. The years of highest percentage TSS retention are typically those where the thermal bar forms in close temporal proximity to the runoff event center of mass (Fig. 9b). Where the bar forms before the event center of mass (negative values on the x-axis of Fig. 9b), in excess of 75% of the load is retained. A 1–2 week lag in thermal bar formation reduces retention to 25–50% and the amount available for trapping falls to < 25% for a lag of three weeks or more. The relationship between TSS retention and lag time is a function of the breadth of the peak of the runoff event and the magnitude of the sink terms (assumed to be constant among years). Model-calculated TSS concentrations return to within 5% of background levels within 3.7 ± 1.8 weeks following cessation of loads. On average, thermal bar formation lags the event center of mass by approximately 3 weeks and thus a relatively small fraction (< 25%) of the event load would typically remain available for trapping. Additional sensitivity analyses were performed to assess the impact of uncertainty in model coefficients on estimates of TSS retention. A ±50% change was applied to each of four model coefficients (current velocity, settling velocity, and longshore and offshore diffusion) and the model was run for years of low (1998) and high (1993) TSS retention. The calculation is most sensitive to settling velocity and offshore diffusion with a variation in the retention estimate of < 5%. The calculation is insensitive to uncertainty in estimates of longshore diffusion and advective transport (current velocity). These simulations suggest that uncertainties in model coefficients do not influence the reliability of estimates of the fraction of the spring runoff event load available for trapping by the thermal bar. A similar calculation, quantification of the mass available for trapping, additionally accommodates

Thermal Bar Retention of Runoff to Lake Superior

FIG. 10. Mass delivered during the spring runoff event TSS load available for trapping at the time of thermal bar formation: (a) distribution over the 33year period of record and (b) as a function of the date of thermal bar formation relative to the SRE center of mass. A negative value for the thermal bar lag indicates that the front formed prior to occurrence of the center of mass. interannual variation in the amount of TSS delivered in the spring runoff event. This effect can be significant. For example, 1987 ranked 19th in %TSS available for trapping, but 30th in mass TSS available, a reflection of the total amount of solids delivered by that year’s spring runoff event. The range in mass remaining (Fig. 10a) over the 33-year period of record is remarkable, ranging from negligible amounts in several years (1967–68, 1987–90) to in excess of 700,000 metric tons in others (1969, 1979). As with percentage remaining, variability in mass remaining is closely linked to the timing of the thermal bar (Fig. 10b) but there is more scatter for mass remaining due to the effect of interannual variations in the magnitude of the runoff event load. The potential significance of the runoff event and

77

thermal bar trapping in stimulating biochemical processes (e.g., bacterioplankton and phytoplankton growth) was explored by examining the degree to which terrestrial inputs elevated TSS concentrations in the nearshore over baseline levels. TSS concentrations were calculated on a daily basis and averaged over the 53-cell network for the entire spring runoff period. Daily average concentrations were then divided by the baseline concentration to yield elevation over baseline as a multiplier, X, of the baseline concentration. Additionally, daily average concentrations (less the baseline concentration) were summed over the duration of the event and expressed as concentration days (CD, mg⋅L–1⋅d), an integrative representation which reflects the impact of both magnitude and duration of runoff. Interannual differences in the predicted impact on the nearshore are quite dramatic (Fig. 11). In 1979 nearshore TSS concentrations were elevated by a factor of 30 and remained elevated for almost 2 months, yielding a CD value of 805 mg⋅L–1⋅d. The response was significantly more muted in 1968, where nearshore TSS levels were elevated for a shorter time and the elevation over background was only a factor of 1–2, yielding a CD value of 104 mg⋅L–1⋅d. Conditions in 1985 and 1992 were intermediate to these extremes and bracket the mean response for the 33-year period of record where CD averaged 300 ± 194 mg⋅L–1⋅d. This level and duration of elevation in nearshore TSS levels is striking. While a direct transfer of these estimates to parameters of biochemical importance (labile organic carbon, phosphorus) would be speculative, such an increase would vastly exceed that related to any other nutrient cycling or supply phenomenon experienced in the lake. An increase of this magnitude would be expected to play an important role in driving bacterioplankton and phytoplankton population dynamics. SUMMARY AND SYSTEM IMPACT Observations of nearshore-offshore gradients in temperature and a variety of water quality parameters have regularly been reported in the Great Lakes literature. The thermal bar has been implicated as a factor governing the development of these gradients, potentially leading to the stimulation of biological activity where materials of terrestrial origin (nutrients) become trapped shoreward of the bar. The degree to which trapping is a significant factor in elevating nearshore concentrations depends on the availability of those materials at the time of thermal

78

Auer and Gatzke

FIG. 11. Total suspended solids concentration in the Lake Superior nearshore during the spring runoff event, expressed for 4 years as elevation over baseline (i.e. X, times the baseline concentration). The overall effect, integration under the curve, is expressed as concentration-days (mg⋅L–1⋅d).

bar formation. Hydrologic records, satellite images and water supply intake data were utilized to establish the timing of thermal bar formation relative to the spring runoff event for a 33-year period of record in the waters of Lake Superior near Michigan’s Keweenaw Peninsula. The center of mass (midpoint) of the runoff event occurs, on average, on 10 April (± 12 days; range, 15 March–29 April) and the event has an average duration of 47 ± 15 days (range, 1777 days). The average date for thermal bar formation is 2 May (± 13 days; range, 6 April–28 May), lagging the midpoint of the runoff event by approximately 3 weeks. The potential exists for significant quantities of terrestrially-derived materials to be transported from the nearshore in the period between the midpoint of

the runoff event and the time of thermal bar formation, impacting the significance of the trapping phenomenon. A mass balance model for total suspended solids in the nearshore region was used to quantify this effect. On average, 23.3 ± 28.5% of the spring runoff event load remained available for trapping. The variability in that amount (less than 15% retained half of the time, more than 50% retained one quarter of the time; Fig. 12a) reflects interannual differences in the relative timing of the runoff event and thermal bar formation. Estimation of the mass remaining additionally accommodates interannual variation in the magnitude of the runoff event. The range here is particularly striking, with negligible amounts in several years and in excess of 700,000 metric tons in others (< 50,000 M.T. more than 50%

Thermal Bar Retention of Runoff to Lake Superior

FIG. 12. Frequency distribution for (a) percentage of spring runoff event load available, (b) mass of spring runoff event load available, and (c) nearshore TSS elevation above background (as concentration days) for the 33-year period of record.

79

of the time, > 250,000 M.T. one fifth of the time; Fig. 12b). It is clear from these results that significant interannual variation in nearshore levels of terrestrially-derived materials may be expected as a result of differences in the magnitude of the runoff event and the timing of thermal bar formation. The spring runoff event is seen to be one of the dominant phenomena impacting the interaction between Lake Superior and its watershed. Based on the USGS flow record and independently-developed relationships between flow and tributary TSS concentrations, it was determined that the spring runoff event delivers ~70% of the annual tributary TSS load to the study system for the 33-year period of record. Model calculations indicate that the runoff event has the potential to increase nearshore TSS concentrations by as much as a factor of 30. The average product of elevation above background and duration of elevated concentrations (concentration days, CD) was 300 ± 194 mg⋅L-1⋅d. For an average spring runoff event duration of 47 days and a postevent residence time of 26 days, this corresponds to a daily average TSS elevation over background of a factor of 4.1 ± 2.7. As with the percent and mass retained, there is considerable variability in the nearshore system concentration response, with a factor of 3 or less predicted 40% of the time and a factor of 7 or more expected 25% of the time (Fig. 12c). This variability reflects interannual differences in the magnitude and duration of the runoff event. It is concluded from this analysis that the spring runoff event has the potential to increase levels of terrestrially-derived substances, including biologically important materials, in nearshore waters by as much as a factor of 30. The magnitude of these predicted increases is such that they could play an important role in mediating food web dynamics. There also exist significant interannual differences in both the magnitude of the runoff event and the nature of the mass transport environment (thermal bar presence) which may influence nearshore-offshore concentration gradients. As a result of these high levels of natural environmental variability, considerable year-to-year differences in the response of bacterioplankton and phytoplankton communities to the spring runoff event may be expected. The significance of the offset in the timing of the spring runoff event and formation of the thermal bar documented here and the level of uncertainty in inter-annual and intra-annual variability in mass transport processes suggest that development of a more comprehensive hydrodynamic modeling of TSS in the coastal waters of Lake Superior is merited.

80

Auer and Gatzke

ACKNOWLEDGMENTS The authors would like to thank Dr. Jim Churchill and Dr. Noel Urban for sharing the results of field studies of the Lake Superior nearshore in support of the modeling exercise. Dr. Dale Robertson kindly provided TSS data from the NASQAN Program which substantially increased the quality of loading estimates. We are indebted to the staff of the Ontonagon, Michigan, water supply system for sharing historical records of water quality at their intake. The authors are especially grateful to Dr. Judy Budd for providing satellite images used in establishing dates of thermal bar formation. The authors offer their appreciation to Dr. William M. Schertzer, Dr. Joseph V. DePinto, and an anonymous review for their comments and constructive criticisms. This paper is a contribution of the Keweenaw Interdisciplinary Transport Experiment in Superior (KITES Project) funded by the National Science Foundation under Grant No. OCE-9712872. REFERENCES APHA (American Public Health Association). 1998. Standard Methods for the Examination of Water and Wastewater, 20 th Edition. American Public Health Association, American Water Works Association, and Water Environment Federation, Washington, D.C. Bennett, E.B. 1978. Characteristics of the thermal regime of Lake Superior. J. Great Lakes Res. 4:310–319. Blanton, J.O. 1986. Coastal frontal zones as barriers to offshore fluxes of contaminants. Rapp. P.-V. Cons. Int. Explor. Mer. 186:18–30. Bolgrien, D.W., and Brooks, A.S. 1992. Analysis of thermal features of Lake Michigan from AVHRR satellite images. J. Great Lakes Res. 18:259–266. Brink, K.H., Bane, J.M., Church, T.M., Fairall, C.W., Geernaert, G.L., Hammond, D.E., Henrichs, S.M., Martens, C.S., Nittrouer, C.A., Rogers, D.P., Roman, M.R., Roughgarden, J.D., Smith, R.L., Wright, L.D., and Yoder, J.A. 1992. Coastal Ocean Processes: A Science Prospectus. Woods Hole Oceanographic Institution Technical Report, WHOI-92-18, Woods Hole, MA. Budd, J.W. 2004. Large-scale transport phenomena in the Keweenaw region of Lake Superior: the Ontonagon plume and the Keweenaw eddy. J. Great Lakes Res. 30 (Suppl. 1):467–480. ——— , Kerfoot, W.C., and Maclean, A.L. 1998. Documenting complex surface temperature patterns from Advanced Very High Resolution Radiometer (AVHRR) imagery of Saginaw Bay, Lake Huron. J. Great Lakes Res. 24:582–594. Chapra, S.C., and Canale, R.P. 1988. Numerical Methods

for Engineers, 2nd Edition. New York, NY: McGrawHill Book Company. Chen, C., Wiesenburg, D., and Xie, L. 1997. Influences of river discharges on biological production over the inner shelf: a coupled biological and physical model of the Louisiana-Texas shelf. J. Marine Res. 55:293–320. Churchill, J.H., Kerfoot, W.C., and Auer, M.T. 2004a. Exchange of water between the Keweenaw Waterway and Lake Superior: characteristics and forcing mechanisms. J. Great Lakes Res. 30 (Suppl. 1):55–63. ——— , Williams, A.J. and Ralph, E.A. 2004b. Sediment resuspension and transport over the shelf and slope off of Lake Superior’s Keweenaw Peninsula. In Review. Journal of Geophysical Research—Oceans. Doerr, S.M., Canale, R.P., and Effler, S.W. 1996. Development and testing of a total phosphorus model for Onondaga Lake. Lake and Reservoir Management 12(1):141–150. Evans, M.S., Hawkins, B.E., and Sell, D.W. 1980. Seasonal features of zooplankton assemblages in the nearshore area of southeastern Lake Michigan. J. Great Lakes Res. 6:275–289. Henrichs, S., Bond, N., Garvine, R., Kineke, G., and Lohrenz, S. 2000. Coastal Ocean Processes (CoOP): Transport and Transformation Processes over Continental Shelves with Substantial Freshwater Inflows. University of Maryland Technical Report TS-237, Cambridge, MD. Huang, J.C.K. 1972. The thermal bar. Geophysical Fluid Dynamics 3:1–25. Hubbard, D.W., and Spain, J.D. 1973. The structure of the early spring thermal bar in Lake Superior. In Proc. 16th Conf. Great Lakes Res., pp. 735–742. Internat. Assoc. Great Lakes Res. Jeong, J. 2002. Distributions and dynamics of copper and particulates in the Keweenaw Peninsula region of Lake Superior. Ph.D. thesis, Michigan Technological University. Klump, J.V., Bedford, K.W., Donelan, M.A., Eadie, B.J., Fahnenstiel, G.L., and Roman, M.R. 1995. Coastal Ocean Processes (CoOP): Cross-Margin Transport in the Great Lakes. University of Maryland Technical Report TS-148, Cambridge, MD. Lathrop, R.G. Jr., Vande Castle, J.R., and Lillesand, T.M. 1990. Monitoring river plume transport and mesoscale circulation in Green Bay, Lake Michigan, through satellite remote sensing. J. Great Lakes Res. 16: 471–484. Likhoshway, Y.V., Kuzmina, A.Y., Potyemkina, T.G., Potyemki, V.L., and Shimaraev, M. 1996. The distribution of diatoms near a thermal bar in Lake Baikal. J. Great Lakes Res. 22:5–14. Menon, A.S., Dutka, B.J., and Jurkovic, A.A. 1971. Preliminary bacteriological investigation of the Lake Ontario thermal bar. In Proc. 14th Conf. Great Lakes Res., pp. 59–68. Internat. Assoc. Great Lakes Res. Moll, R.A., Johengen, T., Bratkovich, A., Saylor, J.,

Thermal Bar Retention of Runoff to Lake Superior Meadows, G., Meadows, L., and Pernie, G. 1993a. Vernal thermal fronts in large lakes: A case study from Lake Michigan. Verh. Internat. Verein. Limnol. 25:65–68. ——— , Bratkovich, A., Chang, W.Y.B., and Peimin, P. 1993b. Physical, chemical, and biological conditions associated with the early stages of the Lake Michigan thermal front. Estuaries 16:92–103. Mortimer, C.H. 1988. Discoveries and testable hypotheses arising from Coastal Zone Color Scanner imagery of southern Lake Michigan. Limnol. Oceanogr. 33: 203–226. Nalewajko, C. 1967. Phytoplankton distribution in Lake Ontario. In Proc. 10 th Conf. Great Lakes Res., pp. 63–69. Internat. Assoc. Great Lakes Res. Noble, V.E., and Anderson, R.F. 1968. Temperature and current in the Grand Haven, Michigan, vicinity during thermal bar conditions. In Proc. 11th Conf. Great Lakes Res., pp. 470–479. Internat. Assoc. Great Lakes Res. Robertson, D.M. 1997. Regionalized loads of sediment and phosphorus to Lakes Michigan and Superior—high flow and long-term average. J. Great Lakes Res. 23:416–439. Rodgers, G.K. 1965. The thermal bar in the Laurentian Great Lakes. In Proc. 8th Conf. Great Lakes Res., pp. 358–363. Internat. Assoc. Great Lakes Res. ——— . 1966. The thermal bar in Lake Ontario. In Proc. 9 th Conf. Great Lakes Res., pp. 369–374. Internat. Assoc. Great Lakes Res. Schertzer, W.M., Rouse, W.R., Blanken, P.D., and Walker, A.E. 2004. Over-lake meteorology and estimated bulk heat exchange of Great Slave Lake in 1998 and 1999. Journal of Hydrometeorology 4(4):649–659. Schnoor, J.L. 1996. Environmental Modeling: Fate and Transport of Pollutants in Water, Air and Soil. New York: Wiley Interscience. Schwab, D.J., Leshkevich, G.A., and Muhr, G.C. 1992. Satellite measurement of surface water temperature in

81

the Great Lakes: Great Lakes Coastwatch. J. Great Lakes Res. 18:247–258. Spain, J.D., Wernert, G.M., and Hubbard, D.W. 1976. The structure of the thermal bar in Lake Superior. J. Great Lakes Res., 2:296–306. Stoermer, E.F. 1968. Nearshore phytoplankton populations in the Grand Haven, Michigan vicinity during thermal bar conditions. In Proc. 11th Conf. Great Lakes Res., pp. 137–150. Internat. Assoc. Great Lakes Res. Strong, A.E. 1974. Great Lakes temperature maps by satellite (IFYGL). In Proc. 17 th Conf. Great Lakes Res., pp. 321–333. Internat. Assoc. Great Lakes Res. Thomann, R.V., and Mueller, J.A. 1987. Principles of Surface Water Quality Modeling and Control. New York: Harper and Row. Tikhomirov, A.I. 1963. The thermal bar of Lake Ladoga. Bulletin (Izvestiya) All-Union Geographical Society 95:134–142. American Geophysical Union Translation, Soviet Hydrology: Selected Papers Number 2. Ullman, D., Brown, J., Cornillon, P., and Mavor, T. 1998. Surface temperature fronts in the Great Lakes. J. Great Lakes Res. 24:753–775. Urban, N.R., Lu, X., Chai, Y., and Apul, D.S. 2004. Sediment trap studies in Lake Superior: insights into resuspension, cross-margin transport, and carbon cycling. J. Great Lakes Res. 30 (Suppl. 1):147–161. Viekman, B.E., and Wimbush, M. 1993. Observations of the vertical structure of the Keweenaw Current, Lake Superior. J. Great Lakes Res. 19:470–479. Walker, W.W., Jr. 1996. Simplified Procedures for Eutrophication Assessment and Prediction: FLUX, PROFILE, and BATHTUB Models. U.S. Army Corps of Engineers. Submitted: 9 February 2004 Accepted: 28 July 2004 Editorial handling: Joseph V. DePinto