1. Introduction

The Great Lakes of North America span 1200 kilometres west to east and the area surrounding them is home to more than one quarter of Canada’s population and approximately one tenth of the United States of America’s population ^[1,2]. The system of lakes and waterways provides a means of transportation and yields water for consumptive uses. It is a source of electrical power and offers a platform for recreational activities such as swimming, fishing, and boating. Some of the world’s largest concentrations of industrial activity take place on the shores of the Great Lakes. It is equally important for farming, supporting approximately 25% of the Canadian agricultural output. Spanning 244,000 km², the Great Lakes Basin is the largest system of surface freshwater in the world supplying the surrounding area with 23,000 km³ of drinking water ^[1,2]. These characteristics make the basin a valuable water resource and one that needs to be properly maintained in order for it to continue to offer these significant benefits.

The study of contaminated sediments can be a proxy indicator of water quality that has implications for assessment, management, remediation and restoration efforts in areas of elevated contamination. More specifically, the measurability of sediment chemistry makes it possible to create guidelines and threshold levels within ecosystems and can lead to proper regulation of the waterways to ensure the health of the ecosystem ^[3]. This article will investigate lead which is one of the more well-known metallic contaminants in the Great Lakes Basin.

---

1.1. Causes of contamination

Contaminated sediments have adverse effects on the aquatic food web and environment. Many pollution stressing factors have and continue to be placed on each of the Great Lakes. Industrial and agricultural expansion as a result of the demand from large populations living on the edges of the lake and river systems has led to increases in soil runoff, discharge from cities, and leachate from disposal sites ^[1]. Small particles in the water are very efficient at picking up certain hazardous chemicals and holding them on their surface. These chemicals and metals include but are not limited to dioxins, phenols, phosphorus, lead, mercury, and pesticides. As these chemicals are bound to the particles, they make their way to the bottom of the lakes/rivers and settle there. Additionally, sediment can be re-suspended through the process of harbour dredging, shipping and navigating, wind and wave action caused by storms, and biotic disturbances which can result in the associated contaminants being mobile and/or more bioavailable in the water column ^[1].

Lead is a nonessential trace element that is toxic to biological life at elevated concentrations. Adverse consequences include an increase in mortality, decrease in benthic invertebrate abundance and diversity, and abnormal development in species living in highly contaminated areas ^[4,5]. The primary historical source of lead to the Great Lakes prior to the mid-1970s was leaded gasoline; lead was used as a gasoline additive beginning in the 1920s. After 1975, there were dramatic reductions in the use of leaded gasoline in concert with the use of unleaded fuels and catalytic converters to reduce combustion engine emissions ^[6]. In addition, the potential adverse human health effects of lead contamination received greater regulatory scrutiny in the 1980s. The dates of the Lake Ontario lake-wide sediment surveys for which extensive lead data are available (1968 and 1998) bracket the period during which lead loadings were dramatically reduced; therefore the kriging analysis and resulting data interpretations provide important knowledge of the impacts of management decisions in restoring and conserving large lake systems that are heavily impacted by human activities.

---

1.2. Threshold effect level and probable effect level

The Canadian federal government specifies Threshold Effect Level (TEL) and Probable Effect Level (PEL) guidelines for sediment contamination. The TEL refers to the concentration below which adverse biological effects are expected to occur rarely, while the PEL defines the level above which adverse effects are expected to occur frequently ^[5,7,8,9,10]. The TEL for lead is 35.0 µg/g while the PEL is 91.3 µg/g. The TEL and PEL have been used to study sediment contamination throughout the Great Lakes region ^{[5,9,10,11,12,13,14,15,16,17,18]}.

---

1.3. Kriging and bathymetry

Kriging interpolation methods were initially developed for mining applications ^[8,11]. They utilize statistical models that incorporate autocorrelation among a group of measured points to create prediction surfaces. Specifically, weights are assigned to measurement points on the basis of distance; in which spatial autocorrelation is quantified in order to weight the spatial arrangement of measured sampling locations ^[19]. By accounting for statistical distance with a variogram model, as opposed to Euclidean distance utilized in deterministic interpolation, customization of the estimation method to a specific analysis is possible ^[11]. If the pattern of spatial continuity of the data can be described visually using a variogram model, it is difficult to improve on the estimates that can be derived in the kriging process. Furthermore, kriging accounts for both the clustering of nearby samples and for their distance to the point to be estimated ^[20]. Given the statistical properties of this method, measures of certainty or accuracy of the predictions can be produced using a cross-validation process. It is arguable that kriging is the optimal interpolation method on the basis of its functionality and its ability to assess error statistically, when generating prediction surfaces ^[11].

For a kriging spatial interpolation model to provide accurate predictions, the Mean Prediction Error (MPE) should be close to 0, the Average Standard Error (ASE) should be as small as possible (below 20), and the Standardized Root-Mean-Squared Prediction Error (SRMSPE) should be close to1 ^[10,11]. If the SRMSPE is greater than 1, there is an underestimation of the variability of the predictions and if the SRMSPE is less than 1, overestimation of the variability is the result ^[19]. Should the errors exceed the accepted values it is necessary to log-transform the data in order to achieve more suitable prediction outcomes.

Bathymetry and sedimentology data offer interesting opportunities for the assessment of contaminant distribution patterns. Limited information can be derived from sediment sample distribution maps using traditional dot or proportional circle representations ^[5,9,12]. Careful design, implementation, and evaluation of the geovisualization of lake bathymetry together with interpolated contaminant surfaces may help in developing improved representations of contamination patterns, their comprehension by stakeholders, and assist in identifying areas where higher pollution concentrations exist ^[5,9]. The use of three-dimensional (3D) geovisualization of interpolated two-dimensional (2D) surfaces in the Great Lakes has only recently been explored ^[5,9]. The main reason for this is that bathymetry data for all of the lakes did not become generally available until recently ^[21]. In this article, 3D refers to the overlay of contamination surfaces that are draped over the bathymetry data ^[5,9].

---

2. Study area

Lake Ontario is the second smallest lake in the Great Lakes Basin (Figure 1). It has an average depth of 86 metres and a water retention time of approximately six years ^[1]. The surface area is approximately 18960 km^{2 [22]}. The western and northwestern Canadian shores are developed with major urban industrial centres such as Hamilton, Toronto and Oshawa. Conversely, the southern U.S. shores are much less urbanized ^[1].

Field research was conducted in 1968 and again in 1998 under the auspices of the Environment Canada Great Lakes Sediment Assessment Program. The core samples that were acquired consisted of the top three centimetres of the surficial sediment which were obtained using the mini-core box sampling method ^[23,24]. Each sample was analyzed for organic contaminants and metals including lead. A total of 249 samples were obtained in 1968, with 68 samples collected in 1998. The characteristics of the datasets are outlined in Table 1.

---

3. Materials and Method

Marvin et al. ^[24] mapped sediment mercury levels throughout the Great Lakes. It is a legacy contaminant originating from a range of sources including atmospheric deposition, industrial discharges, and runoff from urban/industrial watersheds and therefore serves as a proxy for the occurrence and distribution of a range of pollutants. Traditionally, two-dimensional analyses such as dot maps or proportion circle representation have been utilized to examine pollutant levels ^[24,25]. The generation of area estimates allows for expanded spatial analysis of contaminant distribution patterns. The ordinary kriging technique was used to generate lake-wide pollutant distribution surfaces. This method was chosen over other techniques such as Inverse Distance Weighting (IDW) ^[26] or other kriging methods such as Indicator Kriging ^[27] as it has proven useful in similar lake and river analyses ^{[8,10,11,15,16,17,28,29]}. The 1968 and 1998 model statistics are presented in Table 2. It should be noted that log-normalization of the data was necessary as the error statistics for the non-logged models were not satisfactory. The kriged results are presented in Figure 2 (1968) and Figure 3 (1998).

For 1968 and 1998, the statistics indicate slight under predictions of the estimation surfaces due to the SRMSPE values of 1.079 and 1.111 respectively. The models (spherical and exponential) with the best results were however not identical which may be due to the higher number of sample points in 1968 versus 1998. For both years a minimum of one and a maximum of five nearest neighbours were used with 25000/50000 metre and 50000/100000 metre search ranges respectively. The search direction was 90° for each year. All values were determined through experimentation and examination of the model statistics. The kriging map results are contained within a bounding box made up of the outermost points in each direction. This results in some interpolated area values beyond the sediment samples however the trend beyond these points is already determined based on the calculations within the box and is limited by the lake shoreline.

The lake-bottom bathymetry data were obtained from the National Oceanic and Atmospheric Administration (NOAA) website where they have been compiled for all of the Great Lakes ^[21]. The data come in ASCII format and these were converted to a raster grid using the ArcGIS software. Specifically, the data were downloaded in ARC ASCII file format. The data were then imported into ArcGIS. The conversion tool ASCII to Raster was used to convert the data to the bathymetry surface (Figure 4). Lighter tones indicate deeper portions of the lake which may seem counter-intuitive but was the best visualization option. A viewpoint from the west was chosen as this allows for examination of most of the lake-bottom. Due to the steep drop-off of the southern shoreline, a conventional view from a southern viewpoint does not allow for satisfactory 3D geovisualization opportunities.

---

4. Results and Discussion

The sample points and kriged maps were divided into three main categories: below the TEL, TEL to less than the PEL, and greater than or equal to the PEL. Each of these was subsequently divided into three concentration/prediction intervals. Spatial variability was observed in both the 1968 and 1998 datasets. Area analysis of the kriged surfaces was performed to quantify the results in terms of prediction interval area (in km² and percentage) for each class. The results are presented in Tables 3 and 4. The 1968 kriged map encompassed an area of 18961.13 km². Of this, 11100.56 km² or 59.17% of the lake had contamination levels above the PEL. In the highest class (>=128.8 µg/g), an extraordinary 6288.39 km² of 33.52% of the lake was in the highest contamination category. The 1998 results encompassed an area of 17611.65 km². Substantially lower contaminant concentrations were observed with 510.18 km² or 2.90% of the lake in the highest PEL category. The less than PEL areas covered almost 80% of the lake.

Both kriged surfaces were clipped to their common extents (i.e., where they overlapped) and a comparative assessment of changes over time was made (Tables 5 and 6). The resulting analysis area was slightly reduced to 17441.57 km². As stated above, contamination levels in 1998 dropped substantially, especially in areas that were the most heavily polluted and above the Probable Effect Level (6189.07 km² or 62.00% of the common analysis area in 1968 versus 4660.23 km² or 26.72% in 1998). TEL to <PEL areas were more prominent in 1998 than 1968 (making up 53.30% versus 32.71% respectively). In addition, areas below the Threshold Effect Level increased from 922.09km² (5.29%) in 1968 to 3484.22 km² (19.98%) in 1998. The >=PEL categories were mostly smaller in 1998 as compared to 1968. Contamination area reductions were observed in the two highest >=PEL intervals. In the <TEL and TEL to <PEL categories, increases in area were observed indicating lower overall pollution levels in 1998.

Bathymetry features (Figure 5) and annualized lake currents (Figure 6) together with an understanding of key sources of contaminants can assist with explaining the observed patterns. For 1968, a 3D perspective was first generated with the sediment sample point data layered on the bathymetry (Figure 7). The major source of water for Lake Ontario is the Niagara River and Welland Canal system. Noticeable from the point distribution map is that the majority of the >=PEL points are found in the deeper parts of the lake and in the direction of annualized circulation patterns which flow to the east along the southern lake shoreline. It is not however until the kriged surface is layered on the bathymetry (Figure 8) that the distribution patterns become much more interpretable. Three areas of higher contaminant concentrations can be clearly identified in the Niagara, Mississauga, and Rochester Basins. These are interrupted by two underwater sills (the Whitby-Olcott in the western part of the lake and the Scotch-Bonnet in the east). In addition, the Duck-Galloo Sill in the northeastern part of the lake acts as a barrier to sediment transport out through the St. Lawrence River. The shallower northern shoreline mostly has concentrations below the TEL. This is due to a lack of major industrial activity in this area which can act as sources for lead contamination and the fact that the deeper-water depositional areas lie much farther offshore. Resuspension and redistribution of sediments through storm activity and the influence of prevailing current regimes can also affect the observed patterns. The TEL and PEL isolines assist in identifying the main distribution patterns. The bathymetry data greatly assist in the examination and interpretation of the kriged data.

For 1998, the point data were also draped on the bathymetry layer (Figure 9). The smaller sample size is clearly evident. The general pattern however remains similar in that most of the >=PEL sample points are again located in the deeper lake basins or in proximity to known source areas such as the harbor in Hamilton, Ontario. The kriging generated surface (Figure 10) allows for a more direct comparison with the 1968 results despite the smaller input sample dataset. The >=PEL area has diminished in size which given the phase-out of lead as an additive in gasoline during the 1980s in the USA ^[6] and by 1990 in Canada (made under the Canadian Environmental Protection Act that came into force on April 26, 1990) ^[32] is a likely consequence. Industrial activity that was also the source of some lead inputs has also been subject to more stringent environmental controls. Forsythe et al. ^[5] also observed reduced lead contamination levels in Lake St. Clair (located between Lakes Huron and Erie) when comparing 1970 and 1974 datasets to newer 2001 data.

The <TEL area along the northern shore of Lake Ontario in 1998 has become larger and former above PEL areas along the southern shore of the lake are now below the PEL. The three deep lake basins are still above the PEL however the areas in the highest category (>=128.8 µg/g) have been reduced to one area in the Mississauga Basin to the east of the Whitby-Olcott Sill. Concentrations on this sill are now below the TEL whereas in 1968 they were between the TEL and PEL. A similar pattern can also be observed in the area of the Scotch-Bonnet Sill.

Due the distribution of the sediment sample survey points in 1968 and 1998, the interpolated surfaces are not directly comparable (unless they are clipped to common extents) as they do not encompass the same portions of the lake-bottom. It can however be concluded that lead contamination levels were higher in larger portions of the lake in 1968 as compared to 1998. Shallow and underwater ridge areas tend to have lower levels of contamination than deeper lake basins or contaminant inflow areas as fine-grained particulates, i.e., silts and clays, are ultimately deposited in these areas and become a component of the sedimentary record. The 1968 dataset likely provides a more detailed estimation surface as there were more points available for interpolation procedures. Kriging does however allow for the generation of lake-wide contamination surfaces irrespective of the number of points that are available.

The dates of the historical comprehensive sediment survey of Lake Ontario (1968) and the more recent assessment (1998) bracket the periods during which the primary driver of loadings of lead contamination to the lower Great Lakes, i.e., use of leaded gasoline, was substantially reduced. Previous work by Marvin et al. ^[33] using sediment cores showed that accumulation of lead in Lake Ontario bottom sediments peaked in the late 1970s. This reduction in lead loadings is clearly evident by the comparison of both the magnitude of contamination and the associated spatial distributions in 1968 (Figure 8) and 1998 (Figure 10). The distribution of lead contamination in 1968 implicated both atmospheric deposition and runoff from urbanized/industrialized watersheds, including those of the Niagara River, Genesee River and Oswego Rivers, as sources and vectors for loadings of lead from vehicular emissions to Lake Ontario (Figure 8). The spatial distribution of lead contamination in 1998 de-emphasized the major tributaries as vectors for lead loadings, and showed the majority of the residual historical lead contamination to be confined to the deepest areas of the three major depositional basins (Figure 10). The 1998 lead distribution also indicates that sediment monitoring efforts should specifically target these deep water depositional areas in order to provide the most accurate representation of lead contamination in Lake Ontario.

---

5. Conclusion

Bathymetry data allowed for enhanced interpretation and an improved understanding of observed lead pollution patterns. They were especially useful when combined with knowledge of the annualized circulation patterns for the lake. The general trend for contamination in Lake Ontario was studied and found to be one of lower overall lead pollution levels. The areas that are still above the PEL were observed to be much smaller in 1998 compared to 1968. Lead is still being input into the system however at much reduced levels when compared to historical levels.

The results of the comprehensive kriging analysis of Lake Ontario sediment lead contamination could have implications for future survey design and contaminant monitoring. There are currently a number of general approaches used in lake-wide survey design, including site selection based on a grid pattern, site selection based on historical stations in order to ensure as accurate a temporal comparison as possible, a randomized design, or a combination of these concepts. The use of kriging can enable assessment of areas for sampling that represent highest priority and assist with sample site selection. It can also support the determination of the required degrees of temporal and spatial resolution in order to track progress toward satisfactory restoration and/or protection of large aquatic systems.

---

Conflict of interest

The authors declare no conflict of interest.

---