- Email: [email protected]
Numerical modeling of hydrodynamics and tracer dispersion during ice-free period in Lake Winnipeg
Journal of Great Lakes Research 38 (2012) 147–157
Contents lists available at SciVerse ScienceDirect
Journal of Great Lakes Research j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j g l r
Numerical modeling of hydrodynamics and tracer dispersion during ice-free period in Lake Winnipeg Jun Zhao a, Yerubandi R. Rao a,⁎, Leonard I. Wassenaar b a b
Aquatic Ecosystem Management Research Division (AEMRD), Water Science and Technology Directorate, Environment Canada, 867 Lakeshore Rd., Burlington, ON, Canada L7R 4A6 Environment Canada, 11 Innovation Blvd., Saskatoon, SK, Canada S7N 3H5
a r t i c l e
i n f o
Article history: Received 8 April 2010 Accepted 3 December 2010 Available online 21 March 2011 Communicated by David Schwab Index words: Lake Winnipeg 3D hydrodynamic model Circulation Passive tracer Flushing time Dispersion
a b s t r a c t The three-dimensional hydrodynamic Estuary Lake Coastal Ocean Model (ELCOM) was used to study the water circulation and thermal structure of Lake Winnipeg, Canada. To assess the model performance, we simulated circulation and temperature distribution of the lake in 2007 and compared the model results with the field observations of water levels, temperature, and currents. The model showed considerable success in reproducing the thermal structure and circulation. The model predicted isothermal conditions in the South Basin and Narrows for the whole period and produced weak thermal stratification during the early summer in the North Basin. The modeled currents were used to examine the transport, dispersion of passive tracers, and local flushing time in the lake. Simulations using passive tracers qualitatively agreed well with the field measurements of hydrogen isotopes in Lake Winnipeg taken during the study period. The mean transport is towards the north in the lake with two cyclonic gyres in the North Basin and one gyre in the South Basin. The model results show that the local flushing time of Red River waters in the South Basin is about 45 days. Crown Copyright © 2011 Published by Elsevier B.V. on behalf of International Association for Great Lakes Research. All rights reserved.
Introduction Lake Winnipeg is the 10th largest freshwater lake in the world and the 5th largest in Canada. It has two distinct basins, the North Basin (100 km wide) and the South Basin (40 km wide), which are separated by the Narrows, a 2.5 km wide channel. The dual-basin lake is elongated in shape, extends 436 km from north to south, and is relatively shallow with a mean depth of 9 m and 12 m in south and north basins, respectively. With abundant wind energy and shallow water depths, little or no significant thermal stratification occurs in the lake (Kenney, 1979; Hamblin, 1976; Brunskill et al., 1980). The lake is subject to wind driven turbidity that stirs up sediments from the bottom (Brunskill et al., 1980). The Red, Winnipeg, and Saskatchewan rivers are the major rivers flow into Lake Winnipeg, contributing more than 70% of the water received. The Nelson River drains Lake Winnipeg and runs northward into Hudson Bay. Significant changes in water transparency, biological species composition, primary productivity, and sediment chemistry suggest that the lake is on a trajectory of progressive eutrophication that may affect ecosystem sustainability (Lake Winnipeg Stewardship Board, 2006). Stressors on the lake require continued scientific attention to conditions in the lake that are essential to identify and quantify chemical inputs, organisms, and processes and to determine the
⁎ Corresponding author. Tel.: +1 905 336 4785. E-mail address: [email protected] (Y.R. Rao).
effects on the lake ecosystem. One of the primary goals of the Federal Lake Winnipeg Action Plan is to generate basic understanding of the ecology of Lake Winnipeg in order to ensure that resource management decisions lead to restoring Lake Winnipeg and assuring its long term sustainability. The Red and Winnipeg rivers are the primary sources of flow and nutrients to the south basin of Lake Winnipeg, whereas the Saskatchewan River discharges into the north basin. There is little information on how materials and nutrients introduced primarily from the Red and Winnipeg rivers would deposit and mix within the lake in the South Basin. In addition to basin scale agricultural nutrient inputs, the city of Winnipeg also discharges treated sewage water to the lake through the Red River. Brunskill et al. (1980) collected and analyzed several baseline physical, chemical, and biological samples from 50 lake-wide stations during 6 open water cruises. He examined both spatial and temporal variabilities in the physical, chemical, and biological components of the lake. Since 1969, there have been few smaller studies carried out by the Department of Fisheries and Oceans (DFO) and the Province of Manitoba. Einarrson and Lowe (1968) and Hamblin (1976) attempted to model the seiches and surges in the lake. They observed that the water levels may increase up to 2 m at the downwind shore locations during storms. Constrictions at the Narrows cause distinct winddriven circulations in the north and south basins. Inter-basin water exchange is critical in transporting algae or nutrients from one basin to another. The water mass characteristics indicate that physical processes influence the water exchange flows (Brunskill et al., 1980). Water retention times of the south and north basins will be affected
0380-1330/$ – see front matter. Crown Copyright © 2011 Published by Elsevier B.V. on behalf of International Association for Great Lakes Research. All rights reserved. doi:10.1016/j.jglr.2011.02.005
148
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
Fig. 1. Map of instrument locations and selected bathymetric features within the numerical model domain of Lake Winnipeg. River inflows and outflow are identified with arrows.
evolution of thermal stratification and seasonal and transient circulations in the entire lake. Numerical hydrodynamic models on the other hand are well suited to testing the nature of physical processes in complex coastal environments and are used extensively to explain processes that disperse and transport pollutants in the coastal waters and predict their effects on the ecosystem. Models are fundamental to the management of lakes because of the limitation of site measurements over temporal and spatial scales. As part of this project, both field data collection and numerical model development were carried out in Lake Winnipeg. During 2007, we measured currents, water temperature, wind, solar radiation, water levels, waves, and some water quality parameters at fixed mooring stations in the South and North basins and in the Narrows. The Estuary Lake Coastal Ocean Model (ELCOM) has been used in several lakes (Hodges et al., 2000). Recently, this model has also been applied in the Great Lakes and its embayments (Leon et al., 2005; Rao et al., 2009). ELCOM has an advantage because it can be coupled with the water quality simulator Computational Aquatic Ecosystem Dynamics Model (CAEDYM). For these reasons, as a first step, ELCOM was chosen to simulate the circulation and thermal structure patterns in Lake Winnipeg. The primary objective of this paper is to examine the 3D water circulation, temperature structure, and water levels in Lake Winnipeg based on the model results. Field data including water quality parameters and water isotopes collected in 2007 offered the opportunity to carry out a detailed verification of the hydrodynamic model for the lake. Dispersal patterns of nutrients, pollutants, and of marine organisms that are passive or have limited locomotion are highly dependent on the hydrodynamics of the lake. Because of this, as a secondary objective, model currents were used to investigate the transport and dispersion of passive tracers and estimate local flushing time in Lake Winnipeg. The simulated circulation and estimates of mixing and flushing are critical in assessing the fate and transport of water quality constituents in the lake. Materials and methods Field observations
by these exchange dynamics. Kenney (1979) showed that significant water exchange takes place between the two basins through the Narrows. The surface currents driven by surface wind stress and water level oscillations are about 20–90 cm/s in the Narrows. With the large surface area and shallow water depth, the water column in Lake Winnipeg is generally isothermal with some minor stratification established during weak winds in summer; however, most of that information is based on spatially few spot measurements. Environment Canada deployed two temperature moorings in 2006 to examine the development of thermal stratification and its impact in the lake (Rao, unpublished). These studies and ship based measurements carried out during this period suggest that thermal stratification develops in deeper sections and may play a role in the development of hypoxia in the North Basin. However, these studies were limited to very few observations and insufficient to describe the
From May to October, 2007, we obtained a time series measurements of vertical temperature structure, currents, light transmission, and bottom layer oxygen at several moorings. The air temperature, wind speed and direction, relative humidity, and solar radiation data were obtained from a meteorological buoy which was deployed at station 501 in the South Basin and a land station which was established on George Island and was later relocated to Gimli airport. Fig. 1 shows the positions of the mooring stations and land based station at Gimli. Water temperature data were obtained from 6 thermistor strings located at stations of 500, 502, 503, 504, 505, and 506. The Onset Tidbit type thermistors were deployed on moorings at 1–3 m intervals in the vertical. Temperature data were recorded at 10 min interval and it is accurate to the order of 0.2 °C. Because of the missing of the upper part (spar buoy) on retrieval, temperature
Table 1 Mooring positions during 2007 in Lake Winnipeg. Station no.
Type
Instrument/depth
500 501 502 503 504 505 506 Gimli airport
ADCP and temperature chain Meteorology ADCP and temperature chain Temperature chain Temperature chain ADCP and temperature chain Temperature chain Meteorology
Temperature/ADCP (1, 2,4, 6, 8, 10, 11,11 m)
Tidbits in bold were not recovered due to the damage of spar buoys.
Temperature (1, 3, 4, 6,8, 10, 11, 12, 13,14, 15, 16, 17,17.6 m) Temperature (1, 3, 4, 6,8, 10, 11, 12, 13,14, 15, 16 m) Temperature (1, 3, 4, 4,6, 8, 9, 10, 11, 12,13, 14, 15, 16 m) Temperature/ADCP (1, 3, 4, 6, 8, 10, 11, 12, 13,14, 15, 16, 17, 18,16.6, 18.1 m) Temperature (1, 2, 4, 6, 8, 10, 13, 14, 14.4 m)
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
149
Fig. 2. Time series of surface meteorological observations at station 501 in 2007 (a) relative humidity, (b) solar radiation, (c) air temperature, (d) air pressure, (e) wind speed, and (f) wind direction.
loggers at 1, 3, 4, and 6 m were lost at stations 502, 503, and 505. Three current meter moorings (Acoustic Doppler Current Profilers) were installed at 505 in the North Basin, 500 in the South Basin, and 502 in the Narrows, respectively. The velocity profiles were obtained at 1 m bins. The long-term accuracy of velocity profiles obtained from broadband ADCPs are on the order of ±0.2%. An Acoustic Doppler Velocimeter (ADV) mooring was deployed at station 506 in the North
Basin of Lake Winnipeg. The detailed information about the mooring positions is shown in Table 1. In addition, three Environment Canada weather buoys of Weather Buoy System (WBS), 45144 in the North Basin, 45145 in the Narrows, and 45140 in the South Basin, also provided meteorological data. In 2007, three scientific surveys of Lake Winnipeg took place from May 25 to June 15, July 19 to August 10, and September 12 to October
Fig. 3. Time series of daily river inflows and outflow in Lake Winnipeg in 2007. (a) Red River, (b) Winnipeg River, (c) Saskatchewan River, (d) Dauphin River, and (e) Nelson River.
150
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
5. A suite of samples including δD, δ18O were collected from 56, 51, and 53 stations around the lake in the spring, summer, and fall surveys, respectively. In addition, monthly samples were collected by Environment Canada from the Red, Winnipeg, and Saskatchewan rivers to isotopically characterize the riverine inputs. All samples were analyzed for their δD and δ18O values; however, only δD is used for this study.
Numerical model Numerical simulations of water circulation and thermal structure in Lake Winnipeg were conducted with the Estuary model and Lake Coastal Ocean Model (ELCOM). ELCOM is a nonlinear, fully threedimensional, free surface, z-level model (Hodges et al., 2000). It is designed for practical numerical simulations of the hydrodynamics and the thermodynamics of inland and coastal waters. The model solves the unsteady Reynolds-averaged Navier–Stokes equations and scalar transport equations for incompressible flow using Boussinesq approximation and neglecting the non-hydrostatic pressure terms (Hodges and Dallimore, 2006). Free-surface evolution is modeled using the semi-implicit method of Casulli and Cattani (1994). ELCOM uses Arakawa C-grid; velocities are defined on cell faces with the freesurface height and scalar concentrations on cell centers. The fundamental numerical scheme of ELCOM is the TRIM approach of Casulli and Cheng (1992) with some modifications for accuracy. This provides computational efficiency and allows sharper gradients to be maintained with coarse grid resolution. ELCOM models the vertical Reynolds stress terms in momentum and transport equations with a 3D mixed-layer approach derived from the mixing energy budgets (Hodges et al., 2000). In the horizontal direction, turbulence is treated using a constant eddy viscosity coefficient. Density is calculated from salinity and temperature through an algebraic equation, while the salinity and temperature are calculated through differential equations of diffusive transport. The wind momentum is approximated as a uniform distribution over the mixed layer (Imberger and Patterson, 1981). Wind stress is calculated from wind velocity using a bulk formulation: τw = ρa CD juw juw
ð1Þ
where ρa is the density of air taken as 1.2 kg/m3, CD is drag coefficient, and uw is the wind velocity measured 10 m above the water surface. The computations were performed on a horizontal grid of 2 × 2 km (Fig. 1). For the purpose of estimating exchange flows between the north and south basins, the Narrows constriction has to be resolved in the model. A horizontal grid resolution of 2 km was chosen to represent the Narrows (~ width, 2.5 km) with the availability of computational resources. Although this grid resolution may be coarse to resolve the characteristics of exchange flows, it is sufficient for the present application because our primary goal is not the precise quantification of the exchange flow but the general circulation of the lake. The bathymetry was derived from Canadian Hydrographic Service Charts with a mean lake level of 217.44 m. There were 21 levels in the vertical dimension with 1 m resolution for properly simulating temperature distribution. A time step of 300 s was found to be adequate for obtaining numerical stability. Simulations were started from rest, with horizontal free surface and isopycnals. Bottom boundaries were modeled as turbulent benthic boundary layer and free-slip boundary conditions were used for side land boundaries.
Results and discussion Fig. 2 shows the results of a time series of surface meteorological variables recorded at station 501 and WBS buoy 45144 in the North Basin. Because of the instrument failure, the measured solar radiation data were only available from August 9 to September 28 at meteorological station 501. Solar radiation data were not typically recorded at weather stations, but they may be predicted from other meteorological measurements (Davies and Mckay, 1982; Won, 1977). The meteorological records at nearby Gimli airport weather station were used to compute the solar radiation to fill the data gaps. The average wind speed was about 5.2 m/s in the North Basin and 4.5 m/s in the South Basin. The predominant wind direction was towards the north. The major wind event (maximum speed is 14.3 m/s) occurred from day 213 to 214. Heat and momentum fluxes were estimated from these local meteorological variables to provide the surface boundary conditions of the numerical model. There are four major river inflows in Lake Winnipeg, the Red River and the Winnipeg River in the South Basin, the Saskatchewan River and to a lesser extent the Dauphin River in the North Basin. The outflow is the Nelson River from
Fig. 4. Time series of observed and modeled water level deviations from the mean lake level at (a) station 506 in the North Basin and (b) station 502 in the Narrows.
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
the North Basin. The ELCOM model uses the archived hourly hydrometric data provided by Water Survey of Canada (Fig. 3). As Lake Winnipeg is more than 400 km in meridional direction, winds from the North Basin were different from the winds over the South Basin. We provide wind stress based on the wind data collected at 3 buoys, 45144, 45140, and 45145, of WBS and the meteorological station at Gimli (see Fig. 1). ELCOM was developed for spatially uniform forcing, but limited spatial variability can be incorporated by dividing the model domain into three basins, and in each basin, the forcing is set to be uniform and equal to the observations at the nearby station. Wind stress was modeled as a momentum source distributed vertically over the surface wind-mixed layer (Hodges et al., 2000). In order to better represent the smaller streams and unmeasured flows to the lake, we increased the inflows by 10%–20% based on the water balance estimates (Zhang and Rao, 2012). The temperature profiles measured in the first week at stations 505 and 500 were linearly interpolated to the grid points for providing initial temperature conditions of the North Basin and the South Basin, respectively. The salinity was set as constant, 0.2 ppt. Model validation The ELCOM model has been successfully applied in lakes, estuaries and marine coastal areas (Hodges et al., 2000; Laval et al., 2003; Leon et al., 2005, Rao et al., 2009). Most of these studies involved comparison of field
151
and model results for temperature and salinity; however, only limited comparisons of currents and water levels have been carried out (Rao et al., 2009). This is the first time the ELCOM model was used to simulate the large-scale circulation and thermal structure in Lake Winnipeg. Because of the complicated nature of the topography and the bathymetry of the lake, adequate evaluation of ELCOM with field data was essential for gaining confidence of hydrodynamics for water quality simulations. There were 6 temperature moorings and 3 ADCP measurements in the lake (Fig. 1, Table 1). Water level data were only collected at station 506 in the North Basin. However, water level gauges operated by Water Survey of Canada provided data at some coastal stations. Although all data have been used in the validation of the model, only representative visual comparisons have been presented in this section. We first assess ELCOM performance in reproducing the water levels and the depth mean currents in the lake. Water level changes in the lake result from wind-induced setup, inflows and outflows, precipitation, and evaporation. Fig. 4(a) shows the comparison of the model computed and the observed water level deviations from mean lake level at station 506 in North Basin. As the high frequency motions were not well-resolved in the model, the low pass filtered values (three-point moving average) of observed water levels were shown in the figure. The water level data of Water Survey of Canada at Pine Dock station at the Narrows were used to compare with the model results at station 502 (Fig. 4(b)). We measured the model performance in terms of the γ2 value, which is defined as the variance of the
Fig. 5. Time series of northward (a) and eastward (b) wind velocity and time series of northward and eastward components of observed and model calculated currents at depths of 1 m (c and d), 5 m (e and f), and 9 m (g and h) at station 505 in the North Basin.
152
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
model errors (differences between the observations and model results) normalized by the observed variance. 2
γ = VarðO−M Þ = VarðOÞ
ð2Þ
where O and M are the observed and model calculated variables, respectively, and Var is the variance (Thompson and Sheng, 1997). The γ2 value is a measure of the variance of the model error upon the variance of the observations. The smaller the γ2 is, the better the model results fit the observations. The simulated water levels at 2 sites (station 506 in the North Basin and 502 at Narrows) agreed reasonably well with the observations, with γ2 values of about 0.356 at station 506 and 0.397 at station 502. ELCOM simulates water level very well during calm period but underestimates the storm surges and surface seiches during storm events, which could be because of coarse horizontal resolution and uniform basin winds in the model. ELCOM uses implicit scheme for time-stepping; therefore, slightly larger (300 s) time step was used in the model integration. In order to assess if shorter time steps will provide better estimates of short term variations in the water levels during the storm conditions, we conducted sensitivity tests with 60 s and 120 s time steps. However, the predicted water levels are not affected. We also calculated the γ2 value of the depth-averaged currents by using the model simulations and the ADCP measurements at station 505 in the North Basin and station 500 in the South Basin (Fig. 4 of Rao and Zhao, 2010). The γ2 value of N–S component is about 0.74 for the model results at station 505 and 0.83 at station 500. The model predicted the E–W component
less well than the N–S component; the γ2 values are 0.78 at station 505 and 0.93 at station 500. These results indicate that the model is able to capture critical aspects of the depth-averaged (barotropic) dynamics reasonably well. As a further verification, we examined the model's ability to reproduce the time series of vertical structure of observed currents at all ADCP stations. Figs. 5 and 6 show the time series of wind velocity and the comparison of time series of observed and modeled currents at different water depths at stations of 505 and 500, respectively. In general, the low-pass filtered currents appear to be largely driven by prevailing winds. The model calculated surface currents at 505 were in good agreement with field observations (Fig. 5). However, below the surface layer, γ2 values progressively increased, indicating poorer fit. The γ2 values were about 0.68–1.50 for the eastward components and 0.94–3.03 for the north components at station 505 in the North Basin. In the Narrows, the model underestimated the currents (not shown). The discrepancies between the observed and modeled currents could be due to resolving the complexity of topography in the Narrows. The model currents were sensitive to the unresolved small-scale topography with the present resolution. The values of γ2 ranged from 1.09 to 1.70 for the eastward components and 1.08 to 1.17 for the northward components at station 502 in the Narrows. In the South Basin, γ2 values ranged from 0.95 to 1.24 for the eastward components and 0.94 to 1.52 for northward components at station 500 (Fig. 6). We also examined the model results in simulating the vertical structure of temperature. The time-depth distributions of simulated temperatures were compared with observed temperatures at stations
Fig. 6. Time series of northward (a) and eastward (b) wind velocity and time series of northward and eastward components of observed and model calculated currents at depths of 1 m (c and d), 5 m (e and f), and 7 m (g and h) at station 500 in the South Basin.
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
153
Fig. 7. Model calculated time-depth temperature distribution compared with observations at 3 thermistors stations, 500 (a, c), 502 (b, e), and 505 (c, f).
500, 502, and 505 (Fig. 7). At stations 505 and 502, the temperature loggers at 1, 3, 4, and 6 m were lost because of the missing of spar buoys; only thermistors below 6 m depth returned good quality data, whereas at station 500 all thermistors provided hourly values of temperature at several depths. These observed temperatures were linearly interpolated to model levels for comparing with model results. The ELCOM model captured the dominant features of evolving thermal structure of the water column well. The shallow South Basin warmed quickly and remained isothermal for most of the study period (Fig. 7(a) and (d)). The water temperature was isothermal in the Narrows (Fig. 7(b) and (e)). Stratification was only observed in the deeper North Basin in late June but it vanished rapidly by the end of July due to the high wind that occurred during this period (Fig. 7(c) and (f)). The comparisons of the model results with the field observed temperatures demonstrated that the model can accurately depict the observed thermal structure at all stations. To quantify the temperature misfit between the observation and the model, we also calculated the γ2 values for temperature. Table 2 shows the statistical evaluation (γ2 values) of model performance at different water depths of stations 500 and 506. As before, the model predicted temperatures were in good agreement with observations. The γ2 values were low and in the range of 0.069–0.077 at station 506 in the North Basin and 0.048– 0.066 at station 500 in the South Basin.
Table 2 Statistical evaluation (γ2 values) of model performance for water temperature. Station no.
Depth (m)
γ2 value
500
1.0 4.0 8.0 1.0 4.0 8.0
0.077 0.072 0.069 0.066 0.053 0.048
506
Temperature and circulation patterns Water circulation is responsible for the transport and distribution of contaminants in Lake Winnipeg. In order to examine how the contaminants could be distributed throughout the lake, mean summer circulation of the model runs for the period (6 June to 6 September, 2007) were obtained (Fig. 8). The large scale patterns of the summer mean surface circulation produced by the model followed the general wind pattern over the area. In general, the surface currents flowed towards the northeast in the North Basin (Fig. 8(a)), with several small-scale circulation features due to complex topography near the long point and Georges Island. The currents were mainly northward in the Narrows. A weak cyclonic gyre was found in the South Basin, with a southward current along the Winnipeg Beach and a northward current along the Victoria Beach. The mean surface temperature during this period was ~ 17 °C in the North Basin and ~20 °C in the South Basin. Fig. 8(b) shows the modeled depth-averaged circulation in Lake Winnipeg. The depth-averaged circulation in the North Basin was characterized by spatially variable and nearly northward to northeastward currents in the center of the lake and two cyclonic gyres near the Long point. There was an eastward to south-eastward return flow along the north shore to the mouth of the Nelson River. These gyres and in particular the cyclonic gyre in the South Basin may play a role in transporting and mixing of material and nutrients within the lake. The depth-averaged temperature in the North Basin was about 14 °C in the centered deep region and about 17 °C near the coast. In the South Basin and the Narrows, the depth-averaged temperature was about 19 °C. From the present simulations, examples of snapshots of surface temperatures are presented on 12 and 30 July (Fig. 9(a) and (b)). During these days, winds were moderate and less cloudy; therefore, good quality surface temperatures from Moderate-Resolution Imaging Spectroradiometer (MODIS) are available (Liu et al., 2007). The model produced surface temperature distributions that are qualitatively in agreement with satellite data (Fig. 9(c) and (d)). The possible
154
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
Fig. 8. Time mean circulation and temperature at (a) surface and (b) depth-averaged values. Current vectors are plotted at every second grid.
cause of the discrepancy between model solutions and observations is that the satellite observations of skin temperature can be higher than the surface layer temperature as simulated by the model. Modeling water isotope distribution in the lake The stable isotopes of water (δ2D, δ18O) are conservative tracers of the water molecule, and are of longstanding importance in hydrology (Clark and Fritz, 1997). The three main contributing watersheds to Lake Winnipeg differ in their precipitation and water balance history and processes; hence, stable isotopes provide unique conservative tracers of watershed origin. We used the natural abundance values of the incoming rivers and the patterns in the lake to further validate the tracer component of the ELCOM model. The observed δD values across the lake in spring of 2007 (May 25 to June 15) were interpolated to the model grids as the initial concentrations (C0). The δD values for model inflows were based on the measurements conducted at the Red, Winnipeg, Saskatchewan, and Dauphin rivers in 2007 and used as a passive tracer in the numerical model. The model was integrated from day 160 to 248. The governing equations for the concentrations of the hydrogen isotope tracer were the same as for temperature and salinity in the model, except that the model pressure and flow fields were not affected by the tracer concentration. To validate the performance of the model in simulating passive tracers in Lake Winnipeg, we compare the simulated results with observed summer mean (July and August) δD values in the surface water (Fig. 10). The observed δD values in July/August of 2007 (Fig. 10(a)) were characterized by low δD values in the North Basin and high values in the South Basin. The lowest δD values were at the mouth of the Saskatchewan River of about 120 because of the low δD waters from the Saskatchewan River. The low δD values subsequently spread both eastward and southward after turning at the Long Point. In the eastern coastal area of the North Basin, the δD values were
about 85 per mil. Because of the higher δD water inputs from the Red and Winnipeg rivers, the δD values in the South Basin remained as low as 70 per mil. The model results reproduced most of the features of the observed δD spatial distributions (Fig. 10(b)). The simulated δD values exhibit lowest values of about 120 per mil at the Saskatchewan River mouth. However, the model produced slightly low values of δD along the Long Point and slightly underestimated δD values in the South Basin and in the Narrows. Dispersion of passive tracers and local flushing time Given the confidence of the model to reasonably well simulate surface currents, temperature, and the isotopic patterns of Lake Winnipeg, we used the model to examine the advection and dispersion of point-source passive tracers (eg., Red River) and estimate the local flushing time in the South Basin. The idealized passive conservative tracer was released on day 121 over a specific sub-area of the South Basin near the mouth of the Red River. The plume of the Red River is distinct and the mixing zone is typically well established in the south-eastern part of the lake, which is similar to the size of the sub-area in the model (Rao and Zhao, 2010). In this run, the initial (normalized) concentration of the tracer is set to 1.0 over the sub-area and 0.0 outside. Fig. 11 shows the evolution of near-surface concentrations of the passive tracer (Fig. 11(a)). Since there was no additional source (or sink) of tracers added to the model after the initialization, the nearsurface concentration variations shown in the figure were due only to the advection and dispersion by the model currents. By day 151 (30 days after the initialization), some of the passive tracer dispersed northward onto the central region of the South Basin, and the counter-clockwise current flowing along the western shoreline had brought in some tracer free water. The concentration of passive tracer near the mouth of the Red River was diluted by the zero tracer concentration of the river input (Fig. 11(b)). Fig. 11(c) shows that a
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
155
Fig. 9. Model calculated surface temperature on (a) July 12 and (b) July 30, 2007, and Lake surface temperature derived form MODIS data on (c) July 12 and (d) July 30, 2007.
large amount of passive tracer was advected to the Narrows. Due to the weak counter-clockwise gyre, the tracer concentration in the middle of the South Basin was higher than the coastal area by day 181. After 90 days of simulation, the surface concentration of passive tracer
was reduced to about 0.2 in the South Basin as some of the passive tracer was dispersed to the North Basin (Fig. 11(d)). To quantify the effect of the hydrodynamics on the aquatic ecosystem in Lake Winnipeg, we calculated the flushing time from
156
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
Fig. 10. July and August mean (a) observed and (b) model calculated deuterium distribution in Lake Winnipeg.
the volume averaged concentrations (VAC) defined as the volume integrated concentration of tracer concentrations over a specific subarea normalized by the total volume of the sub-area. The local flushing time was defined as the e-folding time (Te) for the temporal decay of the volume averaged concentration (Jouon et al., 2006; Sheng et al., 2009). Fig. 12 shows time series of VAC during the 90 day period between day 121 and 210. The VAC of the passive tracer decreased exponentially with time as the passive tracers are carried and dispersed mainly by the 3D model currents. The time series of VAC can be approximated by the equation, as follows: C = C0 expð−t = Te Þ
ð3Þ
where C is the VAC at time t, C0 is the initial value of VAC and equal to 1.0 at day 121, and Te is the e-folding time. Based on the time series of VAC during the 90 day period shown in Fig. 12, the estimated efolding time Te was about 45 days for the passive tracer released in the southern part of the South Basin. Brunskill et al. (1980) previously estimated a water residence time of about 0.5 years for the South Basin, which was calculated from the volume of the South Basin divided by the difference between the annual river inflow and the annual evaporation. By considering the inter-basin fluxes, Kenney (1979) estimated a water residence time of 0.17 years for the South Basin based on a simplified model of differential mass transport of water in the Narrows. The local flushing time of our study 45 days (~0.12 years) for the southern part of South Basin was similar to Kenney's estimation for the whole South Basin. Summary and conclusions A three dimensional, z-level hydrodynamic model, known as ELCOM, has been used to simulate water levels, temperature, currents,
and the transport of stable isotopes in Lake Winnipeg. This study is a first attempt of using a 3D model for developing high-resolution circulation in Lake Winnipeg. The model was validated against observations of thermal structure, surface currents, and tracer concentrations in 2007 in the lake. The simulated water levels agreed reasonably well with water level observations during typical wind conditions in Lake Winnipeg during the study period, with the normalized variances of model errors around 0.36. However, the model under-predicted water levels during storm events. Although surface currents are predicted well, the sub-surface currents in Lake Winnipeg were not produced well by the numerical model, with the normalized variances of model errors greater than 1. The discrepancies between the observed and modeled currents could be due to coarse horizontal resolution of the model, which may not have resolved the complex topography of the lake. Furthermore, providing higher resolution winds may improve the prediction of storm surges and sub-surface circulation in the lake. The performance of the ELCOM model in predicting the development of temperature structure is very encouraging, with the normalized variances approximately 0.06. The model also reproduced the evolution and development of thermal stratification very well as well as the distribution of natural water isotopic tracers. A summary of mean summer circulation patterns in the lake can be derived from the simulations. The mean circulation at the surface of the lake followed the mean winds over the region. The mean depthaveraged circulation resulted in two cyclonic gyres in the North Basin and one gyre in the South Basin. These gyres will probably play a role in transporting and mixing of water quality constituents within the lake. The simulated δD values agreed reasonably well with the observed deuterium data collected in the lake during 2007. The numerical model was also used to study the transport and dispersion of idealized passive tracers in the South Basin. Based on the evolution of the concentration of passive tracers, we calculated that the local
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
157
Fig. 12. Time series of volume averaged concentrations (VAC) for passive tracers.
constructive comments that helped to improve the original manuscript.
References
Fig. 11. Surface concentrations of passive tracers released in the lower South Basin at days (a) 121, (b) 151, (c) 181, and (d) 211.
flushing time is about 45 days for the southern part of the South Basin. Our modeling in Lake Winnipeg provided us an understanding of the evolution of lake thermal structure and circulation in the lake. Some preliminary understanding of mixing and transport of material from Red River is also obtained. As a next step, we have to include the water quality components (CAEDYM) into the model so that ecosystem response can be examined under different nutrient conditions. Acknowledgments We thank the Center for Water Research, University of Western Australia, for providing access to ELCOM model. We thank the National Water Research Institute's Research Support Branch and Lake Winnipeg Research Consortium for providing assistance in the fieldwork. We also thank two anonymous reviewers for providing
Brunskill, G.J., Elliott, S.E.M., Campbell, P., 1980. Morphometry, hydrology, and watershed data pertinent to the limnology of Lake Winnipeg. Canadian Manuscript Report of Fisheries and Aquatic Sciences, No. 1556. Casulli, V., Cattani, E., 1994. Stability, accuracy and efficiency of a semi-implicit method for three dimensional shallow water flow. Comput. Math. Application. 22, 99–112. Casulli, V., Cheng, R.T., 1992. Semi-implicit finite difference methods for threedimensional shallow water flow. Int. J. Numer. Methods Fluids 15, 629–648. Clark, I., Fritz, P., 1997. Environmental Isotopes in Hydrogeology. Lewis Publishers, NewYork, U.S.A. Davies, J.A., McKay, D.C., 1982. Estimating solar irradiance and components. Sol. Energy 29, 55–64. Einarrson, E., Lowe, A.B., 1968. Seiches and set-up on Lake Winnipeg. Limnol. Oceanogr. 13, 257–271. Hamblin, P.F., 1976. Seiches, circulation, and storm surges of an ice-free Lake Winnipeg. J. Fish. Res. Board Can. 33, 2377–2391. Hodges, B.R., Dallimore, C., 2006. Estuary, Lake and Coastal Ocean Model: ELCOM, User's Guide. Center for Water Research, University of Western Australia Technical Publication, p. 62. Hodges, B.R., Imberger, J., Saggio, A., Winters, K., 2000. Modeling basin scale internal waves in a stratified lake. Limnol. Oceanogr. 45, 1603–1620. Imberger, J., Patterson, J.C., 1981. A dynamic reservoir simulation model—DYRESM 5. In: Fisher, H. (Ed.), Transport Models for Inland and Coastal Waters. Academic Press, pp. 310–361. Jouon, A., Douillet, P., Ouillon, S., Fraunie, P., 2006. Calculation of hydrodynamic time parameters in a semi-open coastal zone using a 3D hydrodynamic model. Cont. Shelf Res. 26, 1395–1415. Kenney, B.C., 1979. Lake surface fluctuations and the mass flow through the Narrows of Lake Winnipeg. J. Geophys. Res. 84 (B3), 1225–1235. Lake Winnipeg Stewardship Board, 2006. Nutrient Literature Review. http://www.lakewinnipeg.org/web/downloads/Nutrient_Literature_Review_LakeWinnipeg_Final _May_2006.pdf. Laval, B., Imberger, J., Hodges, B.R., Stocker, R., 2003. Modelling circulation in lakes: spatial and temporal variations. Limnol. Oceanogr. 48, 983–994. Leon, L., Imberger, J., Smith, R.E.H., Hecky, R.E., Lam, D.C.L., Schertzer, W.M., 2005. Modelling as a tool for nutrient management in Lake Erie: a hydrodynamics study. J. Great Lakes Res. 31 (sup2), 309–318. Liu, J., Hirose, T., Kapfer, M., Bennett, J., McCullough, G., Hochheim, K., Stainton, M., 2007. Operational Water Quality Monitoring over Lake Winnipeg Using Satellite Remote Sensing Data. Conference Proceeding. CRSS/ASPRS 2007 Specialty Conference, Our Common Borders—Safety, Security, and the Environment Through Remote Sensing. Rao, Y.R., Marvin, C.H., Zhao, J., 2009. Application of a numerical model for circulation, temperature, and pollutant distribution in Hamilton Harbour. J. Great Lakes Res. 35, 61–73. Rao, Y.R., Zhao, J., 2010. Numerical simulation of the influence of a Red River flood on circulation and contaminant dispersion in Lake Winnipeg. Nat. Hazardsdoi:10.1007/s11069-010-9534-5 Sheng, J., Zhao, J., Zhai, L., 2009. Examination of circulation, dispersion, and connectivity in Lunenburg Bay of Nova Scotia using a nested-grid circulation model. J. Mar. Sys. 77, 350–365. Thompson, K.R., Sheng, J., 1997. Subtidal circulation on the Scotian Shelf: assessing the hindcasting skill of a barotropic model. J. Geophys. Res. 102, 24987–25003. Won, T., 1977. The simulation of hourly global radiation from hourly reported meteorological parameters—Canadian Prairie area. Paper presented at the 3rd Conference, Canadian Solar Energy Society, Edmonton, Alberta, 23 p. Zhang, W., Rao, Y.R., 2012. Application of a eutrophication model for assessing water quality in Lake Winnipeg. J. Great Lakes Res. 38 (Suppl. 3), 158–173.
.
Contents lists available at SciVerse ScienceDirect
Journal of Great Lakes Research j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j g l r
Numerical modeling of hydrodynamics and tracer dispersion during ice-free period in Lake Winnipeg Jun Zhao a, Yerubandi R. Rao a,⁎, Leonard I. Wassenaar b a b
Aquatic Ecosystem Management Research Division (AEMRD), Water Science and Technology Directorate, Environment Canada, 867 Lakeshore Rd., Burlington, ON, Canada L7R 4A6 Environment Canada, 11 Innovation Blvd., Saskatoon, SK, Canada S7N 3H5
a r t i c l e
i n f o
Article history: Received 8 April 2010 Accepted 3 December 2010 Available online 21 March 2011 Communicated by David Schwab Index words: Lake Winnipeg 3D hydrodynamic model Circulation Passive tracer Flushing time Dispersion
a b s t r a c t The three-dimensional hydrodynamic Estuary Lake Coastal Ocean Model (ELCOM) was used to study the water circulation and thermal structure of Lake Winnipeg, Canada. To assess the model performance, we simulated circulation and temperature distribution of the lake in 2007 and compared the model results with the field observations of water levels, temperature, and currents. The model showed considerable success in reproducing the thermal structure and circulation. The model predicted isothermal conditions in the South Basin and Narrows for the whole period and produced weak thermal stratification during the early summer in the North Basin. The modeled currents were used to examine the transport, dispersion of passive tracers, and local flushing time in the lake. Simulations using passive tracers qualitatively agreed well with the field measurements of hydrogen isotopes in Lake Winnipeg taken during the study period. The mean transport is towards the north in the lake with two cyclonic gyres in the North Basin and one gyre in the South Basin. The model results show that the local flushing time of Red River waters in the South Basin is about 45 days. Crown Copyright © 2011 Published by Elsevier B.V. on behalf of International Association for Great Lakes Research. All rights reserved.
Introduction Lake Winnipeg is the 10th largest freshwater lake in the world and the 5th largest in Canada. It has two distinct basins, the North Basin (100 km wide) and the South Basin (40 km wide), which are separated by the Narrows, a 2.5 km wide channel. The dual-basin lake is elongated in shape, extends 436 km from north to south, and is relatively shallow with a mean depth of 9 m and 12 m in south and north basins, respectively. With abundant wind energy and shallow water depths, little or no significant thermal stratification occurs in the lake (Kenney, 1979; Hamblin, 1976; Brunskill et al., 1980). The lake is subject to wind driven turbidity that stirs up sediments from the bottom (Brunskill et al., 1980). The Red, Winnipeg, and Saskatchewan rivers are the major rivers flow into Lake Winnipeg, contributing more than 70% of the water received. The Nelson River drains Lake Winnipeg and runs northward into Hudson Bay. Significant changes in water transparency, biological species composition, primary productivity, and sediment chemistry suggest that the lake is on a trajectory of progressive eutrophication that may affect ecosystem sustainability (Lake Winnipeg Stewardship Board, 2006). Stressors on the lake require continued scientific attention to conditions in the lake that are essential to identify and quantify chemical inputs, organisms, and processes and to determine the
⁎ Corresponding author. Tel.: +1 905 336 4785. E-mail address: [email protected] (Y.R. Rao).
effects on the lake ecosystem. One of the primary goals of the Federal Lake Winnipeg Action Plan is to generate basic understanding of the ecology of Lake Winnipeg in order to ensure that resource management decisions lead to restoring Lake Winnipeg and assuring its long term sustainability. The Red and Winnipeg rivers are the primary sources of flow and nutrients to the south basin of Lake Winnipeg, whereas the Saskatchewan River discharges into the north basin. There is little information on how materials and nutrients introduced primarily from the Red and Winnipeg rivers would deposit and mix within the lake in the South Basin. In addition to basin scale agricultural nutrient inputs, the city of Winnipeg also discharges treated sewage water to the lake through the Red River. Brunskill et al. (1980) collected and analyzed several baseline physical, chemical, and biological samples from 50 lake-wide stations during 6 open water cruises. He examined both spatial and temporal variabilities in the physical, chemical, and biological components of the lake. Since 1969, there have been few smaller studies carried out by the Department of Fisheries and Oceans (DFO) and the Province of Manitoba. Einarrson and Lowe (1968) and Hamblin (1976) attempted to model the seiches and surges in the lake. They observed that the water levels may increase up to 2 m at the downwind shore locations during storms. Constrictions at the Narrows cause distinct winddriven circulations in the north and south basins. Inter-basin water exchange is critical in transporting algae or nutrients from one basin to another. The water mass characteristics indicate that physical processes influence the water exchange flows (Brunskill et al., 1980). Water retention times of the south and north basins will be affected
0380-1330/$ – see front matter. Crown Copyright © 2011 Published by Elsevier B.V. on behalf of International Association for Great Lakes Research. All rights reserved. doi:10.1016/j.jglr.2011.02.005
148
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
Fig. 1. Map of instrument locations and selected bathymetric features within the numerical model domain of Lake Winnipeg. River inflows and outflow are identified with arrows.
evolution of thermal stratification and seasonal and transient circulations in the entire lake. Numerical hydrodynamic models on the other hand are well suited to testing the nature of physical processes in complex coastal environments and are used extensively to explain processes that disperse and transport pollutants in the coastal waters and predict their effects on the ecosystem. Models are fundamental to the management of lakes because of the limitation of site measurements over temporal and spatial scales. As part of this project, both field data collection and numerical model development were carried out in Lake Winnipeg. During 2007, we measured currents, water temperature, wind, solar radiation, water levels, waves, and some water quality parameters at fixed mooring stations in the South and North basins and in the Narrows. The Estuary Lake Coastal Ocean Model (ELCOM) has been used in several lakes (Hodges et al., 2000). Recently, this model has also been applied in the Great Lakes and its embayments (Leon et al., 2005; Rao et al., 2009). ELCOM has an advantage because it can be coupled with the water quality simulator Computational Aquatic Ecosystem Dynamics Model (CAEDYM). For these reasons, as a first step, ELCOM was chosen to simulate the circulation and thermal structure patterns in Lake Winnipeg. The primary objective of this paper is to examine the 3D water circulation, temperature structure, and water levels in Lake Winnipeg based on the model results. Field data including water quality parameters and water isotopes collected in 2007 offered the opportunity to carry out a detailed verification of the hydrodynamic model for the lake. Dispersal patterns of nutrients, pollutants, and of marine organisms that are passive or have limited locomotion are highly dependent on the hydrodynamics of the lake. Because of this, as a secondary objective, model currents were used to investigate the transport and dispersion of passive tracers and estimate local flushing time in Lake Winnipeg. The simulated circulation and estimates of mixing and flushing are critical in assessing the fate and transport of water quality constituents in the lake. Materials and methods Field observations
by these exchange dynamics. Kenney (1979) showed that significant water exchange takes place between the two basins through the Narrows. The surface currents driven by surface wind stress and water level oscillations are about 20–90 cm/s in the Narrows. With the large surface area and shallow water depth, the water column in Lake Winnipeg is generally isothermal with some minor stratification established during weak winds in summer; however, most of that information is based on spatially few spot measurements. Environment Canada deployed two temperature moorings in 2006 to examine the development of thermal stratification and its impact in the lake (Rao, unpublished). These studies and ship based measurements carried out during this period suggest that thermal stratification develops in deeper sections and may play a role in the development of hypoxia in the North Basin. However, these studies were limited to very few observations and insufficient to describe the
From May to October, 2007, we obtained a time series measurements of vertical temperature structure, currents, light transmission, and bottom layer oxygen at several moorings. The air temperature, wind speed and direction, relative humidity, and solar radiation data were obtained from a meteorological buoy which was deployed at station 501 in the South Basin and a land station which was established on George Island and was later relocated to Gimli airport. Fig. 1 shows the positions of the mooring stations and land based station at Gimli. Water temperature data were obtained from 6 thermistor strings located at stations of 500, 502, 503, 504, 505, and 506. The Onset Tidbit type thermistors were deployed on moorings at 1–3 m intervals in the vertical. Temperature data were recorded at 10 min interval and it is accurate to the order of 0.2 °C. Because of the missing of the upper part (spar buoy) on retrieval, temperature
Table 1 Mooring positions during 2007 in Lake Winnipeg. Station no.
Type
Instrument/depth
500 501 502 503 504 505 506 Gimli airport
ADCP and temperature chain Meteorology ADCP and temperature chain Temperature chain Temperature chain ADCP and temperature chain Temperature chain Meteorology
Temperature/ADCP (1, 2,4, 6, 8, 10, 11,11 m)
Tidbits in bold were not recovered due to the damage of spar buoys.
Temperature (1, 3, 4, 6,8, 10, 11, 12, 13,14, 15, 16, 17,17.6 m) Temperature (1, 3, 4, 6,8, 10, 11, 12, 13,14, 15, 16 m) Temperature (1, 3, 4, 4,6, 8, 9, 10, 11, 12,13, 14, 15, 16 m) Temperature/ADCP (1, 3, 4, 6, 8, 10, 11, 12, 13,14, 15, 16, 17, 18,16.6, 18.1 m) Temperature (1, 2, 4, 6, 8, 10, 13, 14, 14.4 m)
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
149
Fig. 2. Time series of surface meteorological observations at station 501 in 2007 (a) relative humidity, (b) solar radiation, (c) air temperature, (d) air pressure, (e) wind speed, and (f) wind direction.
loggers at 1, 3, 4, and 6 m were lost at stations 502, 503, and 505. Three current meter moorings (Acoustic Doppler Current Profilers) were installed at 505 in the North Basin, 500 in the South Basin, and 502 in the Narrows, respectively. The velocity profiles were obtained at 1 m bins. The long-term accuracy of velocity profiles obtained from broadband ADCPs are on the order of ±0.2%. An Acoustic Doppler Velocimeter (ADV) mooring was deployed at station 506 in the North
Basin of Lake Winnipeg. The detailed information about the mooring positions is shown in Table 1. In addition, three Environment Canada weather buoys of Weather Buoy System (WBS), 45144 in the North Basin, 45145 in the Narrows, and 45140 in the South Basin, also provided meteorological data. In 2007, three scientific surveys of Lake Winnipeg took place from May 25 to June 15, July 19 to August 10, and September 12 to October
Fig. 3. Time series of daily river inflows and outflow in Lake Winnipeg in 2007. (a) Red River, (b) Winnipeg River, (c) Saskatchewan River, (d) Dauphin River, and (e) Nelson River.
150
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
5. A suite of samples including δD, δ18O were collected from 56, 51, and 53 stations around the lake in the spring, summer, and fall surveys, respectively. In addition, monthly samples were collected by Environment Canada from the Red, Winnipeg, and Saskatchewan rivers to isotopically characterize the riverine inputs. All samples were analyzed for their δD and δ18O values; however, only δD is used for this study.
Numerical model Numerical simulations of water circulation and thermal structure in Lake Winnipeg were conducted with the Estuary model and Lake Coastal Ocean Model (ELCOM). ELCOM is a nonlinear, fully threedimensional, free surface, z-level model (Hodges et al., 2000). It is designed for practical numerical simulations of the hydrodynamics and the thermodynamics of inland and coastal waters. The model solves the unsteady Reynolds-averaged Navier–Stokes equations and scalar transport equations for incompressible flow using Boussinesq approximation and neglecting the non-hydrostatic pressure terms (Hodges and Dallimore, 2006). Free-surface evolution is modeled using the semi-implicit method of Casulli and Cattani (1994). ELCOM uses Arakawa C-grid; velocities are defined on cell faces with the freesurface height and scalar concentrations on cell centers. The fundamental numerical scheme of ELCOM is the TRIM approach of Casulli and Cheng (1992) with some modifications for accuracy. This provides computational efficiency and allows sharper gradients to be maintained with coarse grid resolution. ELCOM models the vertical Reynolds stress terms in momentum and transport equations with a 3D mixed-layer approach derived from the mixing energy budgets (Hodges et al., 2000). In the horizontal direction, turbulence is treated using a constant eddy viscosity coefficient. Density is calculated from salinity and temperature through an algebraic equation, while the salinity and temperature are calculated through differential equations of diffusive transport. The wind momentum is approximated as a uniform distribution over the mixed layer (Imberger and Patterson, 1981). Wind stress is calculated from wind velocity using a bulk formulation: τw = ρa CD juw juw
ð1Þ
where ρa is the density of air taken as 1.2 kg/m3, CD is drag coefficient, and uw is the wind velocity measured 10 m above the water surface. The computations were performed on a horizontal grid of 2 × 2 km (Fig. 1). For the purpose of estimating exchange flows between the north and south basins, the Narrows constriction has to be resolved in the model. A horizontal grid resolution of 2 km was chosen to represent the Narrows (~ width, 2.5 km) with the availability of computational resources. Although this grid resolution may be coarse to resolve the characteristics of exchange flows, it is sufficient for the present application because our primary goal is not the precise quantification of the exchange flow but the general circulation of the lake. The bathymetry was derived from Canadian Hydrographic Service Charts with a mean lake level of 217.44 m. There were 21 levels in the vertical dimension with 1 m resolution for properly simulating temperature distribution. A time step of 300 s was found to be adequate for obtaining numerical stability. Simulations were started from rest, with horizontal free surface and isopycnals. Bottom boundaries were modeled as turbulent benthic boundary layer and free-slip boundary conditions were used for side land boundaries.
Results and discussion Fig. 2 shows the results of a time series of surface meteorological variables recorded at station 501 and WBS buoy 45144 in the North Basin. Because of the instrument failure, the measured solar radiation data were only available from August 9 to September 28 at meteorological station 501. Solar radiation data were not typically recorded at weather stations, but they may be predicted from other meteorological measurements (Davies and Mckay, 1982; Won, 1977). The meteorological records at nearby Gimli airport weather station were used to compute the solar radiation to fill the data gaps. The average wind speed was about 5.2 m/s in the North Basin and 4.5 m/s in the South Basin. The predominant wind direction was towards the north. The major wind event (maximum speed is 14.3 m/s) occurred from day 213 to 214. Heat and momentum fluxes were estimated from these local meteorological variables to provide the surface boundary conditions of the numerical model. There are four major river inflows in Lake Winnipeg, the Red River and the Winnipeg River in the South Basin, the Saskatchewan River and to a lesser extent the Dauphin River in the North Basin. The outflow is the Nelson River from
Fig. 4. Time series of observed and modeled water level deviations from the mean lake level at (a) station 506 in the North Basin and (b) station 502 in the Narrows.
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
the North Basin. The ELCOM model uses the archived hourly hydrometric data provided by Water Survey of Canada (Fig. 3). As Lake Winnipeg is more than 400 km in meridional direction, winds from the North Basin were different from the winds over the South Basin. We provide wind stress based on the wind data collected at 3 buoys, 45144, 45140, and 45145, of WBS and the meteorological station at Gimli (see Fig. 1). ELCOM was developed for spatially uniform forcing, but limited spatial variability can be incorporated by dividing the model domain into three basins, and in each basin, the forcing is set to be uniform and equal to the observations at the nearby station. Wind stress was modeled as a momentum source distributed vertically over the surface wind-mixed layer (Hodges et al., 2000). In order to better represent the smaller streams and unmeasured flows to the lake, we increased the inflows by 10%–20% based on the water balance estimates (Zhang and Rao, 2012). The temperature profiles measured in the first week at stations 505 and 500 were linearly interpolated to the grid points for providing initial temperature conditions of the North Basin and the South Basin, respectively. The salinity was set as constant, 0.2 ppt. Model validation The ELCOM model has been successfully applied in lakes, estuaries and marine coastal areas (Hodges et al., 2000; Laval et al., 2003; Leon et al., 2005, Rao et al., 2009). Most of these studies involved comparison of field
151
and model results for temperature and salinity; however, only limited comparisons of currents and water levels have been carried out (Rao et al., 2009). This is the first time the ELCOM model was used to simulate the large-scale circulation and thermal structure in Lake Winnipeg. Because of the complicated nature of the topography and the bathymetry of the lake, adequate evaluation of ELCOM with field data was essential for gaining confidence of hydrodynamics for water quality simulations. There were 6 temperature moorings and 3 ADCP measurements in the lake (Fig. 1, Table 1). Water level data were only collected at station 506 in the North Basin. However, water level gauges operated by Water Survey of Canada provided data at some coastal stations. Although all data have been used in the validation of the model, only representative visual comparisons have been presented in this section. We first assess ELCOM performance in reproducing the water levels and the depth mean currents in the lake. Water level changes in the lake result from wind-induced setup, inflows and outflows, precipitation, and evaporation. Fig. 4(a) shows the comparison of the model computed and the observed water level deviations from mean lake level at station 506 in North Basin. As the high frequency motions were not well-resolved in the model, the low pass filtered values (three-point moving average) of observed water levels were shown in the figure. The water level data of Water Survey of Canada at Pine Dock station at the Narrows were used to compare with the model results at station 502 (Fig. 4(b)). We measured the model performance in terms of the γ2 value, which is defined as the variance of the
Fig. 5. Time series of northward (a) and eastward (b) wind velocity and time series of northward and eastward components of observed and model calculated currents at depths of 1 m (c and d), 5 m (e and f), and 9 m (g and h) at station 505 in the North Basin.
152
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
model errors (differences between the observations and model results) normalized by the observed variance. 2
γ = VarðO−M Þ = VarðOÞ
ð2Þ
where O and M are the observed and model calculated variables, respectively, and Var is the variance (Thompson and Sheng, 1997). The γ2 value is a measure of the variance of the model error upon the variance of the observations. The smaller the γ2 is, the better the model results fit the observations. The simulated water levels at 2 sites (station 506 in the North Basin and 502 at Narrows) agreed reasonably well with the observations, with γ2 values of about 0.356 at station 506 and 0.397 at station 502. ELCOM simulates water level very well during calm period but underestimates the storm surges and surface seiches during storm events, which could be because of coarse horizontal resolution and uniform basin winds in the model. ELCOM uses implicit scheme for time-stepping; therefore, slightly larger (300 s) time step was used in the model integration. In order to assess if shorter time steps will provide better estimates of short term variations in the water levels during the storm conditions, we conducted sensitivity tests with 60 s and 120 s time steps. However, the predicted water levels are not affected. We also calculated the γ2 value of the depth-averaged currents by using the model simulations and the ADCP measurements at station 505 in the North Basin and station 500 in the South Basin (Fig. 4 of Rao and Zhao, 2010). The γ2 value of N–S component is about 0.74 for the model results at station 505 and 0.83 at station 500. The model predicted the E–W component
less well than the N–S component; the γ2 values are 0.78 at station 505 and 0.93 at station 500. These results indicate that the model is able to capture critical aspects of the depth-averaged (barotropic) dynamics reasonably well. As a further verification, we examined the model's ability to reproduce the time series of vertical structure of observed currents at all ADCP stations. Figs. 5 and 6 show the time series of wind velocity and the comparison of time series of observed and modeled currents at different water depths at stations of 505 and 500, respectively. In general, the low-pass filtered currents appear to be largely driven by prevailing winds. The model calculated surface currents at 505 were in good agreement with field observations (Fig. 5). However, below the surface layer, γ2 values progressively increased, indicating poorer fit. The γ2 values were about 0.68–1.50 for the eastward components and 0.94–3.03 for the north components at station 505 in the North Basin. In the Narrows, the model underestimated the currents (not shown). The discrepancies between the observed and modeled currents could be due to resolving the complexity of topography in the Narrows. The model currents were sensitive to the unresolved small-scale topography with the present resolution. The values of γ2 ranged from 1.09 to 1.70 for the eastward components and 1.08 to 1.17 for the northward components at station 502 in the Narrows. In the South Basin, γ2 values ranged from 0.95 to 1.24 for the eastward components and 0.94 to 1.52 for northward components at station 500 (Fig. 6). We also examined the model results in simulating the vertical structure of temperature. The time-depth distributions of simulated temperatures were compared with observed temperatures at stations
Fig. 6. Time series of northward (a) and eastward (b) wind velocity and time series of northward and eastward components of observed and model calculated currents at depths of 1 m (c and d), 5 m (e and f), and 7 m (g and h) at station 500 in the South Basin.
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
153
Fig. 7. Model calculated time-depth temperature distribution compared with observations at 3 thermistors stations, 500 (a, c), 502 (b, e), and 505 (c, f).
500, 502, and 505 (Fig. 7). At stations 505 and 502, the temperature loggers at 1, 3, 4, and 6 m were lost because of the missing of spar buoys; only thermistors below 6 m depth returned good quality data, whereas at station 500 all thermistors provided hourly values of temperature at several depths. These observed temperatures were linearly interpolated to model levels for comparing with model results. The ELCOM model captured the dominant features of evolving thermal structure of the water column well. The shallow South Basin warmed quickly and remained isothermal for most of the study period (Fig. 7(a) and (d)). The water temperature was isothermal in the Narrows (Fig. 7(b) and (e)). Stratification was only observed in the deeper North Basin in late June but it vanished rapidly by the end of July due to the high wind that occurred during this period (Fig. 7(c) and (f)). The comparisons of the model results with the field observed temperatures demonstrated that the model can accurately depict the observed thermal structure at all stations. To quantify the temperature misfit between the observation and the model, we also calculated the γ2 values for temperature. Table 2 shows the statistical evaluation (γ2 values) of model performance at different water depths of stations 500 and 506. As before, the model predicted temperatures were in good agreement with observations. The γ2 values were low and in the range of 0.069–0.077 at station 506 in the North Basin and 0.048– 0.066 at station 500 in the South Basin.
Table 2 Statistical evaluation (γ2 values) of model performance for water temperature. Station no.
Depth (m)
γ2 value
500
1.0 4.0 8.0 1.0 4.0 8.0
0.077 0.072 0.069 0.066 0.053 0.048
506
Temperature and circulation patterns Water circulation is responsible for the transport and distribution of contaminants in Lake Winnipeg. In order to examine how the contaminants could be distributed throughout the lake, mean summer circulation of the model runs for the period (6 June to 6 September, 2007) were obtained (Fig. 8). The large scale patterns of the summer mean surface circulation produced by the model followed the general wind pattern over the area. In general, the surface currents flowed towards the northeast in the North Basin (Fig. 8(a)), with several small-scale circulation features due to complex topography near the long point and Georges Island. The currents were mainly northward in the Narrows. A weak cyclonic gyre was found in the South Basin, with a southward current along the Winnipeg Beach and a northward current along the Victoria Beach. The mean surface temperature during this period was ~ 17 °C in the North Basin and ~20 °C in the South Basin. Fig. 8(b) shows the modeled depth-averaged circulation in Lake Winnipeg. The depth-averaged circulation in the North Basin was characterized by spatially variable and nearly northward to northeastward currents in the center of the lake and two cyclonic gyres near the Long point. There was an eastward to south-eastward return flow along the north shore to the mouth of the Nelson River. These gyres and in particular the cyclonic gyre in the South Basin may play a role in transporting and mixing of material and nutrients within the lake. The depth-averaged temperature in the North Basin was about 14 °C in the centered deep region and about 17 °C near the coast. In the South Basin and the Narrows, the depth-averaged temperature was about 19 °C. From the present simulations, examples of snapshots of surface temperatures are presented on 12 and 30 July (Fig. 9(a) and (b)). During these days, winds were moderate and less cloudy; therefore, good quality surface temperatures from Moderate-Resolution Imaging Spectroradiometer (MODIS) are available (Liu et al., 2007). The model produced surface temperature distributions that are qualitatively in agreement with satellite data (Fig. 9(c) and (d)). The possible
154
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
Fig. 8. Time mean circulation and temperature at (a) surface and (b) depth-averaged values. Current vectors are plotted at every second grid.
cause of the discrepancy between model solutions and observations is that the satellite observations of skin temperature can be higher than the surface layer temperature as simulated by the model. Modeling water isotope distribution in the lake The stable isotopes of water (δ2D, δ18O) are conservative tracers of the water molecule, and are of longstanding importance in hydrology (Clark and Fritz, 1997). The three main contributing watersheds to Lake Winnipeg differ in their precipitation and water balance history and processes; hence, stable isotopes provide unique conservative tracers of watershed origin. We used the natural abundance values of the incoming rivers and the patterns in the lake to further validate the tracer component of the ELCOM model. The observed δD values across the lake in spring of 2007 (May 25 to June 15) were interpolated to the model grids as the initial concentrations (C0). The δD values for model inflows were based on the measurements conducted at the Red, Winnipeg, Saskatchewan, and Dauphin rivers in 2007 and used as a passive tracer in the numerical model. The model was integrated from day 160 to 248. The governing equations for the concentrations of the hydrogen isotope tracer were the same as for temperature and salinity in the model, except that the model pressure and flow fields were not affected by the tracer concentration. To validate the performance of the model in simulating passive tracers in Lake Winnipeg, we compare the simulated results with observed summer mean (July and August) δD values in the surface water (Fig. 10). The observed δD values in July/August of 2007 (Fig. 10(a)) were characterized by low δD values in the North Basin and high values in the South Basin. The lowest δD values were at the mouth of the Saskatchewan River of about 120 because of the low δD waters from the Saskatchewan River. The low δD values subsequently spread both eastward and southward after turning at the Long Point. In the eastern coastal area of the North Basin, the δD values were
about 85 per mil. Because of the higher δD water inputs from the Red and Winnipeg rivers, the δD values in the South Basin remained as low as 70 per mil. The model results reproduced most of the features of the observed δD spatial distributions (Fig. 10(b)). The simulated δD values exhibit lowest values of about 120 per mil at the Saskatchewan River mouth. However, the model produced slightly low values of δD along the Long Point and slightly underestimated δD values in the South Basin and in the Narrows. Dispersion of passive tracers and local flushing time Given the confidence of the model to reasonably well simulate surface currents, temperature, and the isotopic patterns of Lake Winnipeg, we used the model to examine the advection and dispersion of point-source passive tracers (eg., Red River) and estimate the local flushing time in the South Basin. The idealized passive conservative tracer was released on day 121 over a specific sub-area of the South Basin near the mouth of the Red River. The plume of the Red River is distinct and the mixing zone is typically well established in the south-eastern part of the lake, which is similar to the size of the sub-area in the model (Rao and Zhao, 2010). In this run, the initial (normalized) concentration of the tracer is set to 1.0 over the sub-area and 0.0 outside. Fig. 11 shows the evolution of near-surface concentrations of the passive tracer (Fig. 11(a)). Since there was no additional source (or sink) of tracers added to the model after the initialization, the nearsurface concentration variations shown in the figure were due only to the advection and dispersion by the model currents. By day 151 (30 days after the initialization), some of the passive tracer dispersed northward onto the central region of the South Basin, and the counter-clockwise current flowing along the western shoreline had brought in some tracer free water. The concentration of passive tracer near the mouth of the Red River was diluted by the zero tracer concentration of the river input (Fig. 11(b)). Fig. 11(c) shows that a
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
155
Fig. 9. Model calculated surface temperature on (a) July 12 and (b) July 30, 2007, and Lake surface temperature derived form MODIS data on (c) July 12 and (d) July 30, 2007.
large amount of passive tracer was advected to the Narrows. Due to the weak counter-clockwise gyre, the tracer concentration in the middle of the South Basin was higher than the coastal area by day 181. After 90 days of simulation, the surface concentration of passive tracer
was reduced to about 0.2 in the South Basin as some of the passive tracer was dispersed to the North Basin (Fig. 11(d)). To quantify the effect of the hydrodynamics on the aquatic ecosystem in Lake Winnipeg, we calculated the flushing time from
156
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
Fig. 10. July and August mean (a) observed and (b) model calculated deuterium distribution in Lake Winnipeg.
the volume averaged concentrations (VAC) defined as the volume integrated concentration of tracer concentrations over a specific subarea normalized by the total volume of the sub-area. The local flushing time was defined as the e-folding time (Te) for the temporal decay of the volume averaged concentration (Jouon et al., 2006; Sheng et al., 2009). Fig. 12 shows time series of VAC during the 90 day period between day 121 and 210. The VAC of the passive tracer decreased exponentially with time as the passive tracers are carried and dispersed mainly by the 3D model currents. The time series of VAC can be approximated by the equation, as follows: C = C0 expð−t = Te Þ
ð3Þ
where C is the VAC at time t, C0 is the initial value of VAC and equal to 1.0 at day 121, and Te is the e-folding time. Based on the time series of VAC during the 90 day period shown in Fig. 12, the estimated efolding time Te was about 45 days for the passive tracer released in the southern part of the South Basin. Brunskill et al. (1980) previously estimated a water residence time of about 0.5 years for the South Basin, which was calculated from the volume of the South Basin divided by the difference between the annual river inflow and the annual evaporation. By considering the inter-basin fluxes, Kenney (1979) estimated a water residence time of 0.17 years for the South Basin based on a simplified model of differential mass transport of water in the Narrows. The local flushing time of our study 45 days (~0.12 years) for the southern part of South Basin was similar to Kenney's estimation for the whole South Basin. Summary and conclusions A three dimensional, z-level hydrodynamic model, known as ELCOM, has been used to simulate water levels, temperature, currents,
and the transport of stable isotopes in Lake Winnipeg. This study is a first attempt of using a 3D model for developing high-resolution circulation in Lake Winnipeg. The model was validated against observations of thermal structure, surface currents, and tracer concentrations in 2007 in the lake. The simulated water levels agreed reasonably well with water level observations during typical wind conditions in Lake Winnipeg during the study period, with the normalized variances of model errors around 0.36. However, the model under-predicted water levels during storm events. Although surface currents are predicted well, the sub-surface currents in Lake Winnipeg were not produced well by the numerical model, with the normalized variances of model errors greater than 1. The discrepancies between the observed and modeled currents could be due to coarse horizontal resolution of the model, which may not have resolved the complex topography of the lake. Furthermore, providing higher resolution winds may improve the prediction of storm surges and sub-surface circulation in the lake. The performance of the ELCOM model in predicting the development of temperature structure is very encouraging, with the normalized variances approximately 0.06. The model also reproduced the evolution and development of thermal stratification very well as well as the distribution of natural water isotopic tracers. A summary of mean summer circulation patterns in the lake can be derived from the simulations. The mean circulation at the surface of the lake followed the mean winds over the region. The mean depthaveraged circulation resulted in two cyclonic gyres in the North Basin and one gyre in the South Basin. These gyres will probably play a role in transporting and mixing of water quality constituents within the lake. The simulated δD values agreed reasonably well with the observed deuterium data collected in the lake during 2007. The numerical model was also used to study the transport and dispersion of idealized passive tracers in the South Basin. Based on the evolution of the concentration of passive tracers, we calculated that the local
J. Zhao et al. / Journal of Great Lakes Research 38 (2012) 147–157
157
Fig. 12. Time series of volume averaged concentrations (VAC) for passive tracers.
constructive comments that helped to improve the original manuscript.
References
Fig. 11. Surface concentrations of passive tracers released in the lower South Basin at days (a) 121, (b) 151, (c) 181, and (d) 211.
flushing time is about 45 days for the southern part of the South Basin. Our modeling in Lake Winnipeg provided us an understanding of the evolution of lake thermal structure and circulation in the lake. Some preliminary understanding of mixing and transport of material from Red River is also obtained. As a next step, we have to include the water quality components (CAEDYM) into the model so that ecosystem response can be examined under different nutrient conditions. Acknowledgments We thank the Center for Water Research, University of Western Australia, for providing access to ELCOM model. We thank the National Water Research Institute's Research Support Branch and Lake Winnipeg Research Consortium for providing assistance in the fieldwork. We also thank two anonymous reviewers for providing
Brunskill, G.J., Elliott, S.E.M., Campbell, P., 1980. Morphometry, hydrology, and watershed data pertinent to the limnology of Lake Winnipeg. Canadian Manuscript Report of Fisheries and Aquatic Sciences, No. 1556. Casulli, V., Cattani, E., 1994. Stability, accuracy and efficiency of a semi-implicit method for three dimensional shallow water flow. Comput. Math. Application. 22, 99–112. Casulli, V., Cheng, R.T., 1992. Semi-implicit finite difference methods for threedimensional shallow water flow. Int. J. Numer. Methods Fluids 15, 629–648. Clark, I., Fritz, P., 1997. Environmental Isotopes in Hydrogeology. Lewis Publishers, NewYork, U.S.A. Davies, J.A., McKay, D.C., 1982. Estimating solar irradiance and components. Sol. Energy 29, 55–64. Einarrson, E., Lowe, A.B., 1968. Seiches and set-up on Lake Winnipeg. Limnol. Oceanogr. 13, 257–271. Hamblin, P.F., 1976. Seiches, circulation, and storm surges of an ice-free Lake Winnipeg. J. Fish. Res. Board Can. 33, 2377–2391. Hodges, B.R., Dallimore, C., 2006. Estuary, Lake and Coastal Ocean Model: ELCOM, User's Guide. Center for Water Research, University of Western Australia Technical Publication, p. 62. Hodges, B.R., Imberger, J., Saggio, A., Winters, K., 2000. Modeling basin scale internal waves in a stratified lake. Limnol. Oceanogr. 45, 1603–1620. Imberger, J., Patterson, J.C., 1981. A dynamic reservoir simulation model—DYRESM 5. In: Fisher, H. (Ed.), Transport Models for Inland and Coastal Waters. Academic Press, pp. 310–361. Jouon, A., Douillet, P., Ouillon, S., Fraunie, P., 2006. Calculation of hydrodynamic time parameters in a semi-open coastal zone using a 3D hydrodynamic model. Cont. Shelf Res. 26, 1395–1415. Kenney, B.C., 1979. Lake surface fluctuations and the mass flow through the Narrows of Lake Winnipeg. J. Geophys. Res. 84 (B3), 1225–1235. Lake Winnipeg Stewardship Board, 2006. Nutrient Literature Review. http://www.lakewinnipeg.org/web/downloads/Nutrient_Literature_Review_LakeWinnipeg_Final _May_2006.pdf. Laval, B., Imberger, J., Hodges, B.R., Stocker, R., 2003. Modelling circulation in lakes: spatial and temporal variations. Limnol. Oceanogr. 48, 983–994. Leon, L., Imberger, J., Smith, R.E.H., Hecky, R.E., Lam, D.C.L., Schertzer, W.M., 2005. Modelling as a tool for nutrient management in Lake Erie: a hydrodynamics study. J. Great Lakes Res. 31 (sup2), 309–318. Liu, J., Hirose, T., Kapfer, M., Bennett, J., McCullough, G., Hochheim, K., Stainton, M., 2007. Operational Water Quality Monitoring over Lake Winnipeg Using Satellite Remote Sensing Data. Conference Proceeding. CRSS/ASPRS 2007 Specialty Conference, Our Common Borders—Safety, Security, and the Environment Through Remote Sensing. Rao, Y.R., Marvin, C.H., Zhao, J., 2009. Application of a numerical model for circulation, temperature, and pollutant distribution in Hamilton Harbour. J. Great Lakes Res. 35, 61–73. Rao, Y.R., Zhao, J., 2010. Numerical simulation of the influence of a Red River flood on circulation and contaminant dispersion in Lake Winnipeg. Nat. Hazardsdoi:10.1007/s11069-010-9534-5 Sheng, J., Zhao, J., Zhai, L., 2009. Examination of circulation, dispersion, and connectivity in Lunenburg Bay of Nova Scotia using a nested-grid circulation model. J. Mar. Sys. 77, 350–365. Thompson, K.R., Sheng, J., 1997. Subtidal circulation on the Scotian Shelf: assessing the hindcasting skill of a barotropic model. J. Geophys. Res. 102, 24987–25003. Won, T., 1977. The simulation of hourly global radiation from hourly reported meteorological parameters—Canadian Prairie area. Paper presented at the 3rd Conference, Canadian Solar Energy Society, Edmonton, Alberta, 23 p. Zhang, W., Rao, Y.R., 2012. Application of a eutrophication model for assessing water quality in Lake Winnipeg. J. Great Lakes Res. 38 (Suppl. 3), 158–173.
.