AM  Vol.7 No.15 , September 2016
Distribution of Sediment Measurements in Lake Michigan as a Case Study: Implications for Estimating Sediment and Water Interactions in Eutrophication and Bioaccumulation Models
ABSTRACT
Lake Michigan, the sixth largest freshwater lake in the world by surface area, was utilized as a water body for assessment. Field data collected at sampling sites throughout the lake in an intensive monitoring effort were utilized for evaluation of the distribution of sediment measurements. An assessment of sediment nutrient and carbon measurements within Lake Michigan was completed to recognize strata resulting from the hydrodynamics of the system. Nonparametric comparison tests revealed that significant differences exist between measurements of sediment nutrients and organic carbon in the lake using strata based upon water column depth (all results demon-strated a p < 0.05, α = 0.05). Cross-validation analysis was applied to the field-collected samples, revealing that large errors occur when estimating sediment flux of carbon or nutrients at a given location in the lake without considering stratification of the distributions of these measurements. Errors in estimating sediment concentrations of nutrients and carbon specific to a location in the lake demonstrated a statistically significant increase when stratification of sediment measurements wasn’t employed among sites. For example, distributions of errors in estimating all nutrients and organic carbon concentrations, whereby distance squared inverse interpolation methods were applied, demonstrated a statistically significant increase in absence of stratification (all p < 0.001, α = 0.05). These results have implications for characterization, monitoring, and modeling sediment and water interaction as related to eutrophication, as well as to contaminant exposure and bioaccumulation for chemicals within Lake Michigan and large water bodies where stratification of the sediment based upon physics of the system exists.

1. Introduction

A tight coupling exists between sediments and the overlying water column in lakes and coastal ecosystems via sediment resuspension. Within the Great Lakes, resuspension of large inventories of nutrients and contaminants deposited over the past decades presently represents substantial fluxes to the water column that can match or exceed external inputs [1] - [5] . Understanding how inventories of nutrients and contaminants are distributed in the sediment within a water body is an important preliminary step to describing water column and sediment interaction for the formulation of eutrophication and bioaccumulation models that rely on accurate estimates of nutrient and/or contaminant fluxes from sediment sources.

Concentrations of nutrients and carbon in sediments are resultant of the bathymetry and hydrodynamics of a given physical system which have influenced rates of distribution of particles in the system over time. Nutrients (phosphorus, silica, and nitrogen) that are present in the sediment are important for algae growth in both marine and freshwater ecosystems when they are reintroduced into the water column. Eutrophication models are commonly developed for water bodies to understand the dynamics between algae populations, the lower food web, and available nutrient concentrations [6] - [16] . In such models, the sediment plays a very important role as both a source and a sink for nutrients from the water column, with nutrients moving between the sediment and the water column by resuspension, settling, and diffusion. Accurately estimating nutrient fluxes from sediments that are available for algae growth within such models is critical to understand and quantify eutrophication processes.

A precise estimation of nutrient fluxes utilized by algae from sediment sources is, likewise, critically important to understand contaminant dynamics involving algae and thus contaminant dynamics within the lower food web of a water body. Algae populations, the base of the food chain, serve as a food source for predators and represent a large biomass in the food web that is capable of making a significant contribution to cycling of nutrients and carbon.

Particulate organic carbon (e.g., fine-grained sediments, plankton) acts as a carrier for toxic chemicals dependent upon a given chemical’s hydrophobicity. Exposure to contaminants such as PCBs that reside in and cycle between lake sediments and the overlying water column has been evaluated using a carbon-based modeled representation of the toxic contaminant within the water body [10] [17] [18] . These models operate by examining the circulation and distribution of the contaminant within the water body as the contaminant interacts with carbon within the system and is transported between environmental compartments (water, sediment, atmosphere). Such models have been applied to define the exposure concentrations to a specific fish species within the water body [10] [17] [18] , and describe changes in the species’ vitality, such as the ability to survive or reproduce [19] . Alternatively, understanding contaminant cycling between sediment and the overlying water column can be used within a bioaccumulation modeling context designed to estimate toxic chemical concentration within an organism that is accumulated from contaminated prey within the organism’s food web [20] - [23] . Bioaccumulation estimates are often linked to a human health risk [24] - [26] . Thus, the ability to accurately represent the concentration of a contaminant within the sediment of a given water body is critically important in modeling a contaminant’s distribution in the water column, bioaccumulation, biomagnification throughout a food web, effects on an organism’s vitality (ecological risk) and human health risks associated with the contaminant.

In the present study, Lake Michigan is used as a case study. Lake Michigan is the sixth largest freshwater lake in the world by surface area with a hydraulic residence time of approximately 62 years [27] . However, it has been demonstrated using radiotracer studies that the internal removal of particle reactive constituents via sedimentation is very rapid. Radiotracer studies using 239 Pu (half life is approximately 25,000 years) and 137 Cs (half-life is approximately 30.2 years) documented that a high percentage (>95%) of these tracers were removed from the water column and transferred to the sediment within a few years [2] [28] . The concentrations remaining in the water column were primarily the result of an annual cycle of sediment resuspension and redeposition that released constituents from the sediment back into the water [1] [3] . Within the present study, we 1) evaluate the distribution of sediment nutrient and carbon concentrations within Lake Michigan to recognize strata that are relevant to the interaction with the hydrodynamics of the system; 2) provide an estimate of the error that can occur when utilizing sediment measurements within a modeling framework without considering distributions of these measurements; and 3) discuss implications that these results have for modeling sediment-water interaction as they relate to eutrophication, contaminant exposure, and bioaccumulation of chemicals within Lake Michigan and other large water bodies where stratification of sediment is likely to exist.

2. Materials and Procedures

2.1. Sediment Samples

The Lake Michigan Mass Balance Study (LMMB) was conducted by the United States Environmental Protection Agency during the 1994 and 1995 time period. As part of this study, a thorough sampling of Lake Michigan sediment was conducted over a series of cruises commencing in July 1994 and ending in May 1996 [29] . Sediment samples included measurements in mg/g of phosphorus (both total phosphorus and sodium hydroxide extractable phosphorus, PO4NaOH), silica (SiO2), organic nitrogen, and organic carbon. A detailed description of sampling techniques and sample analysis can be found in the Lake Michigan Mass Balance Project Methods Compendium [29] . All data collected as part of the Lake Michigan Mass Balance project were subjected to water quality assurance procedures [29] . Data at any sample station that failed the water quality assurance procedures was not included in the present analysis. In the present study, we utilize data collected at stations in Lake Michigan that were widely dispersed throughout the lake both with respect to latitude and longitude (Figure 1), and station depth (m). This data set comprises the most extensive collection of sediment samples in the lake in the past two decades, and provides a thorough coverage of the lake bed appropriate for application of geospatial methods. The data were collected using box cores,

Figure 1. Sediment measurements of nutrients and carbon were collected at sediment sampling stations in Lake Michigan that were widely dispersed throughout the lake both with respect to latitude and longitude [29] .

gravity cores, and Ponars, and the data utilized in the present study represent the surficial sediment (the top 1 cm of the sediment).

2.2. Estimation of Strata

The distribution of the lake sediment measurements was investigated based upon depth of the water column overlaying the sediment bed at each sampling location. The water depth of interaction with the lake bottom can be estimated by Airy theory (proposed by George Biddell Airy in 1841), and also commonly referred to as linear wave theory [30] . This is a core theory of ocean surface waves used in ocean and coastal engineering. Both the shape and speed of a wave are resultant of the displacement of water particles. In deep water, the shape of the wave is a sine wave. As the wave approaches shore the wave motion is affected by bottom friction at a water depth equal to one-half the wavelength. Thus, the potential for wave and sediment interaction exists: sediments beneath a water column shallower than a critical depth will be influenced by wave-induced current and shear stress, while sediments beneath a water column deeper than the critical depth will not. In applying linear wave theory for deep water gravity waves in Lake Michigan, stratification of the sediment was hypothesized to occur based upon the physics of the system. Wind energies result in dissipation of wave energy that impacts resuspension and settling of nutrients and carbon at depths corresponding to the magnitude of the wave energy. The water depth of interaction with the lake bottom was calculated using the equation [30] - [34] :

(1)

where D is the water column depth at which a wave of a given wavelength begins to interact with the sediment (m), g is acceleration due to gravity (m/s2), and T is wave period (s). The National Oceanic and Atmospheric Administration (NOAA) of the U.S. Department of Commerce maintains a total of two buoys within Lake Michigan that are used to record various measurements of the physical system. Records of NOAA buoy data for Lake Michigan at buoys 45,002 and 45,007 were assessed to determine the distribution of water depth of interaction with the lake bottom using Equation (1) over a period of inception to the year 2002 for each buoy, which thus allowed for a continuous 20 year interval of observations at each buoy. For buoy 45,002, there were a total of 120,743 observations included dating from 1979; whereas, for buoy 45,007 a total of 123,704 observations were recorded dating from 1981. These data represent long term data sets for Lake Michigan deep-water waves. Using these distributions, strata including depositional, transitional, and non-depositional were estimated based upon frequency of water column depth at which a wave of a given wavelength begins interaction with the sediment.

2.3. Statistical Analysis

Once boundaries for each stratum (including depositional, transitional, and non-depo- sitional strata) were estimated for Lake Michigan, the distribution of sediment measurements collected in the field for nutrients including total phosphorus, sodium hydroxide extractable phosphorus, silica, and organic nitrogen (hereafter referred to simply as sediment nutrients) as well as organic carbon were investigated with respect to the strata. Measurements for each nutrient and organic carbon, collected at the 116 field sampling locations, were separated into corresponding strata. Independent statistical measures including descriptive statistics, nonparametric comparisons tests and several novel geospatial statistical methods based upon examination of the continuation, or similarity with respect to spatial distribution, were utilized to examine sediment samples. A set of descriptive statistics including calculation of the range; 1st, 2nd, and 3rd quartiles; outliers; and extreme values [35] were employed to compare the distribution of values between strata for each variable (nutrient or organic carbon) and to compare each strata to the distribution of all measurements for a given variable (the distribution that would occur in the absence of stratification).

A series of nonparametric comparison tests, Wilcoxin sum of ranks tests, were applied to examine if the distribution of values within each strata was significantly different from the distribution of the entire data set; whereby, the entire data set refers to the total unstratified data for each nutrient or organic carbon measurement. For each nutrient or organic carbon sediment data collection, a nonparametric ANOVA, the Kruskal Wallis test, was used to examine if the central tendency of the data collected within each stratum was the same across all strata.

Estimation of strata using Linear Wave Theory and based upon water column depth was further explored by application of novel geospatial statistical methods designed to examine spatial continuity of the sediment data sets [36] . Geospatial statistics tools incorporate spatially dependent features including data continuity and data distribution trends to derive better statistical inferences. Two spatial statistical summaries were derived by Xia and Miller [36] . One of them is:

(2)

where is the average perpendicular distance from each point in the h-scatterplot to the line y = x, is the total number of paired samples separated by a lag of h, (x0, y0) are the coordinates of a point within the h-scatterplot which represents the values of the two samples separated by a lag h, and the summation takes place over all possible pairs of samples with a separation of h. The other is:

(3)

where is a normalized index between 0 and 1, and denotes the bigger coordinate component of (x0, y0). The measurement represents an averaged ratio of the actual variation of paired samples over their maximum possible variation for all paired samples with separation h. The smaller the ratio η is, the more continuous the data set under concern. For any given lag h, = 0 indicates all possible points on the associated h-scatterplot are on the y = x line, implying all paired samples with the same values; = 1 indicates that every pair of samples are maximally different and 0 < < 1 measures the degree of the closeness (or diffuseness) to the diagonal line y = x from all points on the h-scatterplot. The sediment data sets for sediment nutrients and organic carbon were analyzed using these geostatistical measures. For a complete derivation and further description of Equations ((2) and (3)) please refer to Xia and Miller [36] .

2.4. Cross-Validation Analysis Using Estimated Strata

We investigated quantification of the error that would occur in employment of interpolation methods to the Lake Michigan sediment nutrients and organic carbon field- collected data using a stratification based upon depth of the overlying water column versus using no stratification. A version of the Princeton Ocean Model (POM) by Blumberg and Mellor [37] , and later configured for Lake Michigan geometry [38] was utilized for its high resolution bathymetric grid. The POM was developed as a high resolution model for Lake Michigan that operates on a grid system with 44,042 water column cells and 19 sigma layers with the 19th sigma layer of grid cells consisting of 2318 5 km × 5 km cells that lie at the surface of the lake bed. Interpolation methods were used to estimate the distribution of surficial sediment concentrations of the sediment nutrients and organic carbon that could be used in sediment and water exchange with the 19th sigma layer of POM, and could be used for a more complex sediment diagenesis model. In employing interpolation methods to the data grouped into strata for a given nutrient or carbon, a separate interpolation was created for each strata. Cross- validation was applied to estimate the error that would occur if interpolation of a given nutrient or organic carbon data set was completed using stratification versus using no stratification. Two different interpolation methods that are commonly applied for geospatial analysis were chosen for investigation, including inverse distance weighting (i.e. inverse distance squared) and natural neighbor. In application of cross-validation, the data was divided into two groups. One group which we will refer to as the “training sample” was used to develop the model, and a second group called the “prediction set” was used to evaluate the predictive ability of the model [39] . In our study, a training sample equivalent to N-1 was used, whereby N was the total number of measurements collected for a given variable. The prediction set included the single value excluded from the training sample, Pi. A prediction error, Ei, was calculated as the absolute value of the difference between the prediction set, Pi, and the interpolated value generated from the training sample corresponding to the location of the single value in the prediction set, Ti:

. (4)

This process was completed over all samples for a given variable of interest so that a total of N prediction errors were generated (one prediction error corresponding to each combination of training set and prediction set). In considering data splitting, the training set was maximized in order to construct the most robust interpolations using each respective interpolation method. In this way, utilizing cross-validation not only allows an understanding of error attributed to different interpolation methods, but more importantly an estimation of the error that would occur in failing to utilize a stratification of sediment data when employing any given interpolation method [39] .

A series of nonparametric comparison tests, Wilcoxin sum of ranks tests, were applied to examine if the distribution of error associated with the interpolation applied using stratification based upon the depth of the overlying water column was significantly different (smaller than) from the distribution of error associated with using interpolation methods with the data set left unstratified. Nonparametric comparison tests were used to make comparisons for stratified and unstratified data sets for the sediment nutrients and organic carbon used with both inverse distance squared and natural neighbor interpolation methods.

3. Results

Data for recorded wavelength of surface waves in Lake Michigan were analyzed from two U.S. Department of Commerce, NOAA buoys, buoy 45,002 located in the northern basin and buoy 45,007 found in the southern basin of the lake. A distribution of 244,447 observations was examined to estimate the frequency of overlying water column depth at which a wave of a given wavelength began interaction with the sediment, (see Equation (1)). From these observations, three strata were defined, a non-depositional stratum representing the overlying water column depth range of 0 to 40 m, a transitional stratum with a depth of the overlying water column ranging from greater than 40 m to 100 m, and a depositional stratum containing overlying water column depths greater than 100 m (Figure 2). The non-depositional zone represents the overlying water column depths at which approximately the 98th percentile of recorded wavelengths at both buoys indicated an interaction with the underlying sediment begins to occur (for buoy 45007 it represents the 98.3089th percentile and for buoy 45002 it represents the 98.0545th percentile). The depositional zone characterizes the overlying water column depths at which no interaction with the underlying sediment is expected to occur (there were zero occurrences for water column depths greater than 100 m for which a wave of a given wavelength was estimated to begin interaction with the sediment). The transitional zone represents the water column depths between the non-depositional and depositional zones, where irregular surface wave interaction with the underlying sediments occurs.

Sediment measurements of nutrients and organic carbon collected at the sampling locations in Lake Michigan were separated into strata based upon depth of the overlying water column including 0 to 40 m, greater than 40 m to 100 m, and greater than 100 m. A set of descriptive statistics was used to examine the distribution of sediment measurements within each stratum and also for the total collection of data points as an unstratified whole (Table 1). Descriptive statistics included calculation of the 1st, 2nd, and 3rd quartiles as well as the range of the data. For the sediment nutrients and organic carbon field-collected data, the stratum based upon water column depth provided a very distinct separation of values. For example, boxplots of the sediment measurements for total phosphorus separated into strata based upon overlying water column depth and also left unstratified are illustrated in Figure 3, and it becomes obvious that the measures of central tendency as well as the interquartile range of the boxplots are very different. A similar pattern was found for the other sediment nutrients and organic carbon.

Wilcoxon sum of ranks tests were used to examine if the distribution of values within each stratum was significantly different from the distribution within the entire data set (Table 1). Results demonstrated that for all of the sediment nutrients and organic carbon, the distribution of values within each stratum was significantly different from the

Figure 2. Frequency of overlying water column depth at which surface wave interaction with the underlying sediment begins to occur as estimated from U.S. Department of Commerce, NOAA buoy data collected in Lake Michigan, and Linear Wave Theory [30] .

Figure 3. Boxplots for total phosphorus sediment concentration (mg/g) that has been grouped into strata based upon overlying water column depth of Non-depositional (0 to 40 m), Transitional (>40 to 100 m) and Depositional (>100 m) as compared to the data left unstratified (All Data).

distribution of the entire unstratified data set (all p < 0.05, α=0.05).

For each nutrient or carbon sediment data collection, a nonparametric ANOVA, the Kruskal Wallis test, was also used to examine if the central tendency of the data collected within each stratum was the same. As would be expected after considering the values in Table 1, results demonstrate that for all of the sediment nutrients and organic carbon the central tendency within each stratum is not the same for each given data set (all p < 0.001, α = 0.05). Thus, for each nutrient or organic carbon field-collected data set, not only is the distribution of values within each strata significantly different from the distribution of values within the unstratified data set as a whole, but also there is a significant difference among the central tendencies of the strata.

Application of geospatial statistical methods for investigation of stratified and non- stratified distributions of the sediment nutrients and organic carbon further quantitatively justified the stratification based upon overlying water column depth. Results for the two spatial statistical summaries that were derived by Xia and Miller [30] including ρ(h), the average perpendicular distance from each point in the h-scatterplot to the line y = x with respect to the separation lag h; and η(h), a normalized index which represents an averaged ratio of the actual variation over the maximum possible variation for all paired samples with separation h were summarized in Table 2. Both tests resulted in a smaller quantification for stratified data sets, indicating the improvement over the unstratified data sets in terms of spatial harmony.

A cross-validation approach was applied in conjunction with the model grid of the Princeton Ocean Model [37] [38] to determine the impact on estimating lake-wide distributions of sediment concentrations for nutrients and organic carbon using a stratification into three water depth zones (0 to 40 m, greater than 40 m to 100 m, and greater than 100 m), as well as using data sets left unstratified (Table 3). Results demonstrate that significantly larger errors occur using interpolations that do not take into account the distribution of the sediment nutrients and organic carbon within the sediment of the lake (Table 3). As an example, Figure 4 illustrates a comparison of error expressed as a percentage of the prediction set (% Pi) for a spatial dependent interpolation (inverse distance squared) for total phosphorus (mg phosphorus/g) constructed using strata versus constructed without using strata. The interquartile range of error associated

Table 1. Statistical comparison of field-collected data for nutrients and carbon in Lake Michigan sediments.

*Wilcoxon Sum of Ranks test; **Kruskal-Wallis test; (ND) = Non-Depositional Zone, (T) = Transitional Zone, (D) = Depositional Zone, All Data = No Strata.

Table 2. Summary of geospatial statistical results for nutrients and carbon from Lake Michigan sediments [36] .

Figure 4. A comparison of error expressed as a proportion of the prediction set (proportion of Pi) for a spatial dependent interpolation (inverse distance squared) for total phosphorus (mg phosphorus/g) constructed using strata versus constructed without using strata.

with application of inverse distance squared for interpolation of total phosphorus constructed without using strata was 24% error to 234% error with extreme values at over 1300% error. This can be compared to the interquartile range of error associated with application of distance squared inverse for interpolation of total phosphorus constructed using strata being 8% error to 83% error. Comparison of error expressed as a percentage of the prediction set (% Pi) for inverse distance squared interpolations for the other sediment nutrients and organic carbon showed a similar trend (Table 3). Application of nonparametric statistical comparison tests revealed that for all nutrients and carbon there was a statistically significant difference between the distribution of error obtained in applying the interpolation with versus without stratification based upon overlying water column depth (all p < 0.001, α = 0.05).

As an additional example, Figure 5 illustrates a comparison of error expressed as a percentage of the prediction set (% Pi) for a spatial dependent interpolation (natural neighbor) for organic carbon (mg organic carbon/g) constructed using strata versus constructed without using strata. The interquartile range of error associated with application of natural neighbor for interpolation of organic carbon constructed without using stratification is 23% error to 192% error. This can be compared to the interquartile range of error associated with using stratification based upon overlying water column depth being 11% error to 89% error. Similar patterns for distribution of error associated with stratified versus non-stratified interpolations using natural neighbor methods existed for the sediment nutrients. Application of nonparametric statistical comparison tests revealed that for all nutrients and organic carbon there was a statistically significant difference between the distribution of error obtained in applying the interpolation with versus without stratification based upon overlying water column depth (all p < 0.05, α = 0.05 except for silica which was significant at the 10% significance level; Table 3).

4. Discussion

The field data collected and analyzed in the present study are unique in terms of the quantity and quality of samples obtained, and the spatial extent of the sampling effort in that it covers the entire Lake Michigan sediment bed [29] . When considering a large

Figure 5. A comparison of error expressed as a proportion of the prediction set (proportion of Pi) for a spatial dependent interpolation (natural neighbor) for organic carbon (mg organic carbon/g) constructed using strata versus constructed without using strata. Note that due to the range of the values, those calculated as extreme values are not illustrated on this box and whiskers plot.

Table 3. Results from application of cross-validation analysis [39] for modeling field-collected data for nutrients and carbon in Lake Michigan sediments.

*Comparison of error from each interpolation method applied to all data stratified versus compared to all data unstratified using Wilcoxon Sum of Ranks test.

scale sampling project of Lake Michigan sediment nutrients and organic carbon that encompasses a range of the entire lake bed, there is enormous value for eutrophication and bioaccumulation research projects that rely heavily on an accurate representation of lake-wide distributions of sediment nutrients and organic carbon.

Sediment distributions of nutrients are critical to understanding eutrophication in water quality studies. Nutrients from sediment sources can act as a food source for algae populations, and large algal blooms have the ability to result in hypoxic or anoxic conditions that are detrimental to fish populations. In addition, eutrophication has an important role in the recycling of toxic chemicals, such as polychlorinated biphenyls (PCBs) in Lake Michigan, that bind to organic particles (algal organic carbon). Several models have been formulated and applied to Lake Michigan to manage the fate and transport of nutrients as well as plankton dynamics on a lake-wide scale [6] [12] [15] [16] [40] - [43] . A characteristic of eutrophication models is an attempt to represent the aquatic system using an understanding of interactions between nutrients, plankton, and sediments. The more accurate the dynamics for nutrient concentrations exchanged between water and sediment, the more credible the model results. Within these modeling studies, phytoplankton growth is modeled as a function of available nutrients, temperature, and light [6] [10] [12] [44] - [47] :

(5)

where kg is growth rate, is optimum growth rate (time−1), is light growth dependency, is temperature growth dependency, and is nutrient growth dependency where (for diatom algae, for example):

(6)

where is available silica concentration, ksat-Si is half-saturation coefficient for silica uptake, NH4 is ammonia concentration, NO3 is nitrate concentration, ksat-N is half- saturation coefficient for nitrogen uptake, P is available phosphorus concentration, and ksat-P is half-saturation coefficient for phosphorus uptake. While our present study doesn’t address the importance of precisely estimating temperature or light within phytoplankton productivity equations or bioavailability, our results demonstrate that when considering sediment nutrient concentration distributions (thus nutrient concentrations available for phytoplankton uptake and growth from sediment sources) and their relationship to eutrophication, it is necessary to stratify the sediment bed of the lake. As an example, we utilized the bathymetric grid of the Princeton Ocean Model for Lake Michigan [37] [38] in conjunction with a cross-validation approach to determine the impact on estimating lake-wide distributions of sediment concentrations for nutrients and organic carbon within 2318 cells representing the lake bottom, both using a stratification into three zones based upon overlying water column depth as well as using data sets left unstratified. Without stratifying total phosphorus in the sediment, for example, the predicted sediment concentration would have significantly larger errors (as large as 1300% in using either inverse distance squared or natural neighbor interpolation methods), and these errors would then carry over into predictions of sediment and water exchange of total phosphorus (Table 2). Examination of interpolation techniques applied to the other nutrients and organic carbon exemplified a similar increase of error that would be encountered and utilized in sediment and water exchange if the modeled distributions of sodium hydroxide extractable phosphorus, organic nitrogen, biogenic silica, or organic carbon in the sediment were left unstratified (Table 2).

The current study results have important applications in understanding contaminant fate and transport processes within Lake Michigan (which are not independent of eutrophication processes) for toxic chemicals that are hydrophobic. Such chemicals, including polychlorinated hydrocarbons (PCBs) in Lake Michigan, are more likely to be sorbed or bound to organic carbon particles when they are in the water column or sediment of a lake. Organic carbon particles that are distributed in the sediment and which are continually reintroduced to the water column via resuspension processes driven by wave dynamics in the nondepositional zone and sometimes the transitional zone provide a continual exchange of toxic chemicals such as PCBs from the sediment into the water column where aquatic organisms are exposed to them. The bioavailable concentration of toxic chemicals that are sorbed or bound to organic carbon particles resulting from the exchange between the sediment and the overlying water column will depend upon physical processes such as settling, resuspension, and diffusion rates that drive sediment and water exchange, as well as processes that are biological in nature including bioturbation and irrigation. In addition, chemical processes such as partitioning and redox reactions will also affect these exposure concentrations. A generally accepted equation for the chemical exchange between a fish and its exposure environment is [20] - [23] :

(7)

where is the rate of chemical increment in fish per unit time (μg/kg/d), FW is flux of chemical uptake from water (μg/kg/d), FP is flux of chemical uptake from prey items (μg/kg/d), FE is flux of chemical elimination via respiration, and FG is flux of chemical reduction by growth dilution (μg/kg/d). For Equation (7), our study is important in estimating both the FW and FP terms. The chemical flux entering a fish from water is expressed as a product of fish ventilation rate and the chemical concentration in the water [20] - [23] :

(8)

where EC is the chemical gill transfer coefficient, KV is the gill ventilation rate (l/kg- fish/d), and CW is the chemical concentration in water (μg of chemical/l water). Further, the chemical flux absorbed by fish from diet, FP, is a result of ingestion rate of the fish and chemical concentration in the diet [20] - [23] :

(9)

where β is the chemical assimilation efficiency, KF is the food ingestion rate (i.e. g prey/ g bodyweight/d) and CP is chemical concentration in the prey (μg of chemical/g prey). It is easy to conclude from examination of Equations ((7)-(9)) that accurate estimation of both CW and CP are dependent upon an accurate estimation of concentration of organic carbon in the sediment for hydrophobic and nonionic chemicals. Further, estimation of CP, in the case of the lower food web of the lake where prey items are likely to be phytoplankton and zooplankton, will require an accurate representation of eutrophication dynamics and thus credible estimates of sediment and water exchange of nutrients.

The persistence of PCBs in the sediment of large lake systems and their bioaccumulation in fish species such as salmon and lake trout has resulted in concern for the toxic effects on humans through the consumption of fish [26] . Within Lake Michigan, PCBs have been identified as a chemical of concern related to health advisories for human consumption of fish. Additionally, there is evidence that toxic chemicals that are hydrophobic and bind to organic carbon, such as dioxins and PCBs, present an ecological risk to populations of aquatic organisms. Exposure to such chemicals has been suggested to result in declines in fish populations due to effects on reproduction and thus recruitment. For example, 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCCD) and structurally related chemicals that act through a common aryl hydrocarbon receptor (AHR)- mediated mechanism of action (including some PCB congeners) found in Lake Ontario sediments are a potential contributing factor to the lake trout population decline in Lake Ontario through lack of recruitment due to exposure [19] . Understanding how carbon is distributed within a water body (both water column and sediment) can thus aid in forecasting the patterns that would be expected for chemical exposure for a contaminant that is nonionic and hydrophobic. Also, with respect to modeling applications, water quality models are often carbon-based models in which the kinetics of carbon within these models function in determining the dynamics of hydrophobic toxic chemicals. For instance, in applying a carbon-based water quality model in Lake Michigan to investigate PCB dynamics that includes a sediment bed, specific attention should be given to congeners that exhibit a high hydrophobicity, because understanding how organic carbon is distributed in sediments will provide insights into modeling the distributions of these congeners both in the water column and the sediment of the lake over time.

5. Conclusions and Recommendations

The results presented in the current study have direct implications for managing and understanding the Great Lakes. The National Oceanic and Atmospheric Administration (NOAA), National Data Buoy Center (NDBC) maintains buoys in each of the Great Lakes that monitor for wave heights in addition to other parameters that have been utilized in analysis of Great Lakes wave patterns [48] . An estimation of stratification of the sediment into non-depositional, transitional, and depositional zones within each of the Great Lakes will provide insight into differences in distribution of nutrients and contaminants between these large water bodies. There is a very different distribution of water depths among the Great Lakes with the average water depth for Lake Erie, Lake Huron, Lake Michigan, Lake Ontario, and Lake Superior being 19, 59, 85, 86, and 147 m respectively [49] . Based upon the results of the current study, a difference in distribution of water depths would in turn result in a difference in the sizes and shapes of the sediment zones among the Great Lakes. Further, toxic chemicals that are nonionic and hydrophobic and thus bind to carbon, for example PCBs, exist in all the Great Lakes (as well as other large lake systems). The present findings demonstrate the importance of using appropriate sampling designs that take into consideration the impact of the physics of the system on particle distribution, and that can account for spatial variability, and in turn allow for more accurate estimation of site specific exposure concentration of toxic chemicals in large lake systems such as the Great Lakes. For instance, the Lake Michigan Mass Balance study utilized 3 biota zones to estimate bioaccumulation of PCB congeners [23] [29] . These biota boxes were large enough to encompass multiple strata as predicted using linear wave theory and the quantitative scores developed in the present study. Careful consideration should be taken to understand distribution of sediment measurements with respect to the boundaries of such biota zones.

In expanding upon the work of the present study, models have been developed to translate site specific exposures to toxic chemicals within a given geographic location into population level effects for various fish species, given that exposures can be linked to vital rates [50] - [53] . The ability to predict the distribution of sediment concentrations of toxic chemicals at specific geographic locations will contribute to more precise ecological risk assessments corresponding to specific contaminated sites within a water body, and therefore a better assessment of ecological risk.

Acknowledgements

This study was conducted in support of model development for the Invasive Species Initiative in the Office of Research and Development (ORD) of the U.S. EPA. The authors wish to acknowledge the Great Lakes National Program Office for their efforts in the Lake Michigan Mass Balance Study, as well as the many cooperators during the study. The authors also wish to acknowledge NOAA, GLERL for their cooperation and collaboration. We thank Doug Endicott for a valuable review of this manuscript. This paper has been reviewed according to ORD guidelines, but the statements made do not represent views of the USEPA, nor does mention of trade names indicate endorsement by the Federal Government.

Cite this paper
Miller, D. , Xia, X. , Huang, W. and Rossmann, R. (2016) Distribution of Sediment Measurements in Lake Michigan as a Case Study: Implications for Estimating Sediment and Water Interactions in Eutrophication and Bioaccumulation Models. Applied Mathematics, 7, 1846-1867. doi: 10.4236/am.2016.715153.
References
[1]   Eadie, B.J., Chambers, R.L., Gardner, W.S. and Bell, G.L. (1984) Sediment Trap Studies in Lake Michigan: Resuspension and Chemical Fluxes in the Southern Basin. Journal of Great Lakes Research, 10, 307-321.
http://dx.doi.org/10.1016/S0380-1330(84)71844-2

[2]   Eadie, B.J. and Robbins, J.A. (1987) Sources and Fates of Aquatic Pollutants. In: Hites, R. and Eisenreich, S., Eds., Advances in Chemistry Series, No. 216, American Chemical Society, Washington DC, 319-364.

[3]   Robbins, J.A. and Eadie, B.J. (1991) Seasonal Cycling of Trace Elements 137Cs, 7Be, and 239+240Pu in Lake Michigan. Journal of Geophysical Research, 96, 17081-17104.
http://dx.doi.org/10.1029/91JC01412

[4]   Brooks, A.S. and Edgington, D.N. (1994) Biogeochemical Control of Phosphorus Cycling and Primary Production in Lake Michigan. Limnology and Oceanography, 39, 961-968.
http://dx.doi.org/10.4319/lo.1994.39.4.0961

[5]   Lee, C., Schwab, D.J. and Hawley, N. (2005) Sensitivity Analysis of Sediment Resuspension Parameters in Coastal Area of Southern Lake Michigan. Journal of Geophysical Research, 110, Article ID: C03004.
http://dx.doi.org/10.1029/2004jc002326

[6]   Rodgers, P.W. and Salisbury, D. (1981) Water Quality Modeling of Lake Michigan, and Consideration of the Anomalous Ice Cover for 1976-1977. Journal of Great Lakes Research, 7, 467-480.
http://dx.doi.org/10.1016/S0380-1330(81)72072-0

[7]   Bierman, V.J. and Dolan, D.M. (1986) Modeling of Phytoplankton in Saginaw Bay: I. Calibration Phase. Journal of Environmental Engineering, 112, 400-414.
http://dx.doi.org/10.1061/(ASCE)0733-9372(1986)112:2(400)

[8]   Bierman, V.J. and Dolan, D.M. (1986) Modeling of Phytoplankton in Saginaw Bay: II. Post-Audit Phase. Journal of Environmental Engineering, 112, 415-429.
http://dx.doi.org/10.1061/(ASCE)0733-9372(1986)112:2(415)

[9]   DiToro, D.M., Thomas, N.A., Herdendorf, C.E., Winfield, R.P. and Connolly, J.P. (1987) A Post Audit of a Lake Erie Eutrophication Model. Journal of Great Lakes Research, 13, 801-825.
http://dx.doi.org/10.1016/S0380-1330(87)71692-X

[10]   Cerco, C. and Cole, T. (1993) Three-Dimensional Eutrophication Model of Chesapeake Bay. Journal of Environmental Engineering, 119, 1006-1025.
http://dx.doi.org/10.1061/(ASCE)0733-9372(1993)119:6(1006)

[11]   Jain, R. and De Pinto, J.V. (1996) Modeling as a Tool to Manage Ecosystems under Multiple Stresses: An Application to Lake Ontario. Journal of Aquatic Ecosystem Stress and Recovery, 5, 23-40.
http://dx.doi.org/10.1007/BF00691727

[12]   Chen, C., Rubao, J., Schwab, D.J., Beletsky, D., Fahnenstiel, G.L., Jiang, M. Johengen, T. H., Vanderploeg, H., Eadie, B., Budd, J.W., Bundy, M.H., Gardner, W., Cotner, J. and Lavrentyev, P.J. (2002) A Model Study of the Coupled Biological and Physical Dynamics in Lake Michigan. Ecological Modelling, 152, 145-168.
http://dx.doi.org/10.1016/S0304-3800(02)00026-1

[13]   Scavia, D., Justice, D. and Bierman Jr., V.J. (2004) Reducing Hypoxia in the Gulf of Mexico: Advice from Three Models. Estuaries and Coasts, 27, 419-425.
http://dx.doi.org/10.1007/BF02803534

[14]   Bierman, V.J., Kaur Jr., J., DePinto, J.V., Feist, T.J. and Wilks, D.W. (2005) Modeling the Role of Zebra Mussels in the Proliferation of Blue-Green Algae in Saginaw Bay, Lake Huron. Journal of Great Lakes Research, 31, 32-55.
http://dx.doi.org/10.1016/S0380-1330(05)70236-7

[15]   Pauer, J.J., Taunt, K.W., Melendez, W. and Kreis Jr., R.G. (2007) Resurrection of the Lake Michigan Eutrophication Model, MICH1. Journal of Great Lakes Research, 33, 554-565.
http://dx.doi.org/10.3394/0380-1330(2007)33[554:ROTLME]2.0.CO;2

[16]   Miller, D.H., Kreis, R.G., Huang, W. and Xia, X. (2010) Development of a Lower Food Web Ecosystem Model for Investigating Dynamics of the Invasive Species Bythotrephes longimanus in Lake Michigan. Biological Invasions, 12, 3513-3524.
http://dx.doi.org/10.1007/s10530-010-9748-1

[17]   Limno-Tech Inc (2004) LOTOX2 Model Documentation: In Support of Development of Load Reduction Strategies and a TMDL for PCBs in Lake Ontario. Technical Report to New England Interstate Water Pollution Control Commission, Lowell, and US Environmental Protection Agency-Region 2, New York.

[18]   Endicott, D.D., Richardson, W.B. and Rossmann, R. (2006) PCBs Modeling Overview. In: Rossmann, R., Ed., Results of the Lake Michigan Mass Balance Study: PCBs Modeling Report, Part 1, Chapter 2, US Environmental Protection Agency, Grosse Ile, 16-25.

[19]   Cook, P., Robbins, J.A., Endicott, D.D., Lodge, K.B., Guiney, P.D., Walker, M.K., Zabelo, E.W. and Peterson, R.E. (2003) Effects of Aryl Hydrocarbon Receptor-Mediated Early Life Stage Toxicity on Lake Trout Populations in Lake Ontario during the 20th Century. Environmental Science & Technology, 37, 3864-3877.
http://dx.doi.org/10.1021/es034045m

[20]   Gobas, F.A.P.C. (1993) A Model for Predicting the Bioaccumulation of Hydrophobic Organic Chemicals in Aquatic Food Webs: Application to Lake Ontario. Ecological Modelling, 69, 1-17.
http://dx.doi.org/10.1016/0304-3800(93)90045-T

[21]   Rudstam, L.G., Binkowski, F.P. and Miller, M.A. (1994) A Bioenergetic Model for Analysis of Food Consumption Patterns by Bloater in Lake Michigan. Transactions of the American Fisheries Society, 123, 344-357.
http://dx.doi.org/10.1577/1548-8659(1994)123<0344:ABMFAO>2.3.CO;2

[22]   Morrison, H.A., Gobas, F.A.P.C., Lazar, R. and Haffner, G.D. (1996) Development and Verification of a Bioaccumulation Model for Organic Contaminants in Benthic Invertebrates. Environmental Science & Technology, 30, 3377-3384.
http://dx.doi.org/10.1021/es960280b

[23]   Zhang, X. (2006) Part 4: LM2 Toxic. In: Rossmann, R., Ed., Results of the Lake Michigan mass Balance Study: PCBs Modeling Report, US Environmental Protection Agency, Grosse Ile, 216-452.

[24]   Buck, G.M., Mendola, P., Vena, J.E., Server, L.E., Kostyniak, P., Greizerstein, H., Olson, J. and Stephen, F.D. (1999) Paternal Lake Ontario Fish Consumption and Risk of Conception Delay, New York State Angler Cohort. Environmental Research, 80, S13-S18.
http://dx.doi.org/10.1006/enrs.1998.3926

[25]   Kostyniak, P.J., Stinson, C., Greizerstein, H.B., Vena, J., Buck, G. and Mendola, P. (1999) Relation of Lake Ontario Fish Consumption, Lifetime Lactation, and Parity to Breast Milk Polychlorobiphenyl and Pesticide Concentrations. Environmental Research, 80, S166-S174.
http://dx.doi.org/10.1006/enrs.1998.3939

[26]   US Environmental Protection Agency (1999) Polychlorinated Biphenyls (PCBs) Update: Impact on Fish Advisories. US Environmental Protection Agency, Office of Water, Washington DC.

[27]   Quinn, F.H. (1992) Hydraulic Residence Times for the Laurentian Great Lakes. Journal of Great Lakes Research, 18, 22-28.
http://dx.doi.org/10.1016/S0380-1330(92)71271-4

[28]   Robbins, J.A. and Edgington, D. (1975) Determination of Recent Sedimentation Rates in Lake Michigan Using Pb-210 and Cs-137. Geochimica et Cosmochimica Acta, 39, 285-304.
http://dx.doi.org/10.1016/0016-7037(75)90198-2

[29]   US Environmental Protection Agency (1997) Lake Michigan Mass Balance Study (LMMB) Methods Compendium, Volume 1: Sample Collection Techniques. US Environmental Protection Agency, Great Lakes National Program Office, Chicago.

[30]   Craik, A.D.D. (2004) The Origins of Water Wave Theory. Annual Review of Fluid Mechanics, 36, 1-28.
http://dx.doi.org/10.1146/annurev.fluid.36.050802.122118

[31]   US Army Coastal Engineering Research Center (1975) Shore Protection Manual, 3 Volumes. US Army Corps of Engineers, US Army Coastal Engineering Research Center, Fort Belvoir.

[32]   Rossmann, R. and Seibel, E. (1977) Surficial Sediment Redistribution by Wave Energy: Element-Grain Size Relationships. Journal of Great Lakes Research, 3, 258-262.
http://dx.doi.org/10.1016/S0380-1330(77)72257-9

[33]   Dean, R.G. and Dalrymple, R.A. (1991) Water Wave Mechanics for Engineers and Scientists. Advanced Series on Ocean Engineering, 2.
http://dx.doi.org/10.1142/1232

[34]   Rossmann, R. (2006) Part 1: Introduction. In: Rossmann, R., Ed., Results of the Lake Michigan Mass Balance Study: PCBs Modeling Report, US Environmental Protection Agency, Grosse Ile, 110-119.

[35]   Moore, D.S. and McCabe, G.P. (1999) Introduction to the Practice of Statistics. 3rd Edition, WH Freeman & Co., New York.

[36]   Xia, X. and Miller, D.H. (2011) The Stratification Analysis of Sediment Data for the LMMBP Study. Journal of Data Science, 9, 181-203.

[37]   Blumberg, A.F. and Mellor, G.L. (1987) A Description of a Three-Dimensional Coastal Ocean Circulation Model. In: Heaps, N.S., Ed., Coastal and Estuarine Sciences, Book 4, American Geophysical Union, 1-6.

[38]   Beletsky, D. and Schwab, D.J. (1998) Modeling Thermal Structure and Circulation in Lake Michigan. Proceedings of the 5th International Conference, Pine Mountain, 1-5 June 1998, 511-522.

[39]   Neter, J., Kutner, M.H., Nachtsheim, C.J. and Wasserman, W. (1996) Applied Linear Statistical Models. 4th Edition, WCB McGraw-Hill, New York.

[40]   Canale, R.P., Depalma, L.M. and Vogel, A.H. (1975) A Plankton Based Food Web Model for Lake Michigan. In: Canale, R.P., Ed., Modeling Biochemical Processes in Aquatic Ecosystems, Science Publishers, Ann Arbor, 33-74.

[41]   Chapra, S.C. and Sonzogni, W.C. (1979) Great Lakes Total Phosphorus Budget for the Mid 1970s. Journal of the Water Pollution Control Federation, 51, 2524-2533.

[42]   Scavia, D., Lang, G.A. and Kitchell, J.F. (1988) Dynamics of Lake Michigan Plankton: A Model Evaluation of Nutrient Loading, Competition, and Predation. Canadian Journal of Fisheries and Aquatic Sciences, 45, 165-177.
http://dx.doi.org/10.1139/f88-018

[43]   Lesht, B.M., Fontaine, T.D. and Dolan, D.M. (1991) Great Lakes Total Phosphorus Model: Post-Audit and Regionalized Sensitivity Analysis. Journal of Great Lakes Research, 17, 3-17.
http://dx.doi.org/10.1016/S0380-1330(91)71337-3

[44]   Monod, J. (1949) The Growth of Bacterial Cultures. Annual Review of Microbiology, 3, 371-394.
http://dx.doi.org/10.1146/annurev.mi.03.100149.002103

[45]   Odum, E. (1971) Fundamentals of Ecology. 3rd Edition, WB Saunders, Philadelphia.

[46]   Thomann, R.V. and Di Toro, D.M. (1975) Mathematical Modeling of Phytoplankton in Lake Ontario, Part 1: Model Development and Verification. US Environmental Protection Agency, Grosse Ile.

[47]   DiToro, D.M. and Connolly, J.P. (1980) Mathematical Models of Water Quality in Large Lakes, Part 2: Lake Erie. US Environmental Protection Agency, Grosse Ile.

[48]   Plattner, S., Mason, D.M., Leshkevich, G.A., Schwab, D.J. and Rutherford, E.S. (2006) Classifying and Forecasting Coastal Upwellings in Lake Michigan Using Satellite Derived Temperature Images and Buoy Data. Journal of Great Lakes Research, 32, 63-76.
http://dx.doi.org/10.3394/0380-1330(2006)32[63:CAFCUI]2.0.CO;2

[49]   US Environmental Protection Agency (2008) The Great Lakes Environmental Atlas and Resource Book. US Environmental Protection Agency and Government of Canada.
http://www.epa.gov/glnpo/atlas/index.html

[50]   Brown, A.R., Riddle, A.M., Cunningham, N.L., Kedwards, T.J., Shillabeer, N. and Hutchinson, T.H. (2003) Predicting the Effects of Endocrine Disrupting Chemicals on Fish Populations. Human and Ecological Risk Assessment, 9, 761-788.
http://dx.doi.org/10.1080/713609966

[51]   Miller, D.H. and Ankley, G.T. (2004) Modeling Impacts on Populations: Fathead Minnow (Pimephales promelas) Exposure to the Endocrine Disruptor 17B-Trenbolone as a Case Study. Ecotoxicology and Environmental Safety, 59, 1-9.
http://dx.doi.org/10.1016/j.ecoenv.2004.05.005

[52]   Miller, D.H., Jensen, K.M., Villeneuve, D.L., Kahl, M.D., Makynen, E.A. and Durhan. E.J. (2007) Linkage of Biochemical Responses to Population-Level Effects: A Case Study with Vitellogenin in the Fathead Minnow (Pimephales promelas). Environmental Toxicology and Chemistry, 26, 521-527.
http://dx.doi.org/10.1897/06-318R.1

[53]   Miller, D.H., Tietge, J.E., McMaster, M.E., Munkittrick, K.R., Xia, X. and Ankley, G.T. (2013) Assessment of Status of White Sucker (Catostomus commersoni) Populations Exposed to Bleached Kraft Pulp Mill Effluent. Environmental Toxicology and Chemistry, 32, 1592-1603.

 
 
Top