A critical step in evaluating and selecting contaminated sediment remediation options is the simulation of a long-term future trend of exposure concentrations that will occur under natural conditions without active in situ remediation. This ‘natural attenuation’ forecast provides a baseline reference against which to compare alternative sediment management options. Model development and application at a number of contaminated sediment sites has indicated that one of the most important parameters for accurately simulating natural attenuation is the depth of the upper mixed (i.e., ‘active’) sediment layer. The upper mixed sediment layer essentially acts as a biologically available contaminant reservoir, determining risks to ecological and human health. Contaminant concentrations in the upper mixed sediment layer are dynamic and determined by: the sediment inventory created by historical loadings to the sediment; a balance of fate and transport processes both into and out of the mixed layer at the sediment water interface; and burial of contaminants from the mixed layer into deeply buried ‘inactive’ sediments. Model sensitivity analyses at one of our sediment management sites, the Lower Fox River, has demonstrated that the rate at which polychlorinated biphenyl exposure decreases with time under natural processes is inversely proportional and very sensitive to the depth of mixing in the surficial sediments of the system. These studies have shown that increasing the assumed mixed layer depth in the model from 10 cm to 30 cm increases the simulated recovery half-time—as measured by the rate of change of the upper layer sediment polychlorinated biphenyl concentration—from about 20 years to more than 100 years. In the case of the Lower Fox River, radioisotope and polychlorinated biphenyl concentration depth profiling supports a 10 cm mixed depth layer. Implementing a 30 cm mixed depth in the model smears the existing polychlorinated biphenyl profile, artificially mixing higher concentrations in the 10 to 30 cm layer with lower concentrations in the 0 to 10 cm layer, resulting in the increased simulated recovery times.
Numerical models of chemical fate and transport have proved to be very valuable tools for informing management decisions on contaminated sediment sites (NRC, 2001). In particular, these models have provided long-term forecasts of contaminant concentrations in water, sediments and biota in response to external forcing functions and internal processes that impact sediment and chemical transport and fate. Among the most important model applications is forecasting the recovery of the system following elimination of significant external contaminant loads while implementing no active remediation of the in situ contaminated sediment. This ‘monitored natural recovery’ (MNR) scenario not only serves as a management option in its own right, but it serves as a reference against which to compare active remedial options such as dredging or capping of contaminated sediments. The MNR scenario, using a forecast of the most likely long-term hydrometeorological forcing conditions, can also serve as a baseline of the predicted system response for comparison to the response predicted for extreme events (100-year flood, high wind events, low water levels) for which there is no actual experience.
Both the NRC (2001) report and recent draft EPA guidance on contaminated sediment remediation recommends that the evaluation of management options, including MNR, should be conducted using a risk-based approach. In most contaminated sediment remediation studies, both ecological and human health risk pathways of concern are driven either directly or indirectly by contaminants in the surface sediments of the system. This, of course, assumes that external contaminant loads have been eliminated and the system is in a recovery phase of its typical life cycle (Wolfe and DePinto, 2001). Given this situation, the recovery rate of contaminants in the surface sediments of a system will depend upon the rates of the processes depicted in Figure 1 relative to each other. These processes determine the three pathways for recovery of surface sediments: burial of sediment and associated contaminant; loss of contaminant from the surface sediments via sediment-to-water transfer processes (both resuspension and non-resuspension chemical fluxes); and chemical decay/biotransformation. The overall rate at which these processes serve to reduce the biologically available contaminant concentration in the upper mixed sediment layer determines whether MNR is a viable remediation option for a given system.
This paper presents a diagnostic analysis of the fate and transport model (FoxSim) we developed for use in the Lower Fox River RI/FS study (LTI, 2002a). The purpose of the model diagnosis was to investigate the effect of mixed layer thickness on natural recovery rates in the Lower Fox River system.
Fate and transport model development
The Lower Fox River in Wisconsin extends from Lake Winnebago to its mouth at Green Bay. Major industrialization in the Fox River valley has resulted in significant contamination of the river and its sediments. In particular, river sediments contain high levels of polychlorinated biphenyls (PCBs), which were released into the river by pulp and paper companies during the 1955 to 1970 period. The sub-reach between DePere Dam and Green Bay is a freshwater estuary impacted by wind-driven seiches in Lake Michigan and Green Bay. Considerable data collection efforts have been directed at developing a comprehensive conceptual model for this reach. Extensive sampling of sediment physical properties and PCB concentrations was conducted in 1995, with analyses conducted for the top 3 meters of the sediment column, including the following depth intervals: 0–5 cm, 5–10 cm, 10–30 cm, and 30–50 cm. Physical and chemical data indicate that a large percentage of PCB mass downstream of DePere Dam is buried deeply, with 85% of the PCB mass below 30 cm and 73% below 50 cm. The PCB concentrations in the 10–30 cm interval are 3 times higher on average than concentrations in the 0–10 cm interval (WDNR, 1999). Analysis of 137Cs and 210Pb profiles from approximately 30 radioisotope cores collected in this reach indicate that the reach is net depositional with an average burial rate of 0.5 cm y−1 (range: 0.3–2.1 cm y−1) and a sediment mixed depth ranging between 6–12 cm.
A numerical model of Lower Fox River sediment and PCB fate and transport (FoxSim) has been developed and applied for the purpose of forecasting long-term sediment recovery and risk reduction at the site (LTI, 2002a). The FoxSim model is a time-variable mass balance model of water, solids and PCBs in the water column and sediment bed. FoxSim is a modified version of the USEPA WASP4/TOXI4 mass balance model (Ambrose et al., 1988). The modifications include all relevant revisions included in WASP5 (Ambrose et al., 1993) as well as enhancements to sediment bed handling that eliminate the potential for artificial mixing between sediment bed layers. Various versions of this USEPA model have been applied to address the fate and transport of PCBs for a number of contaminated river systems around the country (e.g., Hudson River, New York; Buffalo River, New York; Saginaw River, Michigan; Kalamazoo River, Michigan; and St. Clair River, Michigan).
The FoxSim application simulates three state variables: inorganic/abiotic solids, organic/biotic solids, and total PCBs. A conceptual diagram of the FoxSim mass balance model is presented in Figure 2. FoxSim accounts for the important fate and transport processes affecting the solids and PCBs loaded into and existing within the system, including advective and dispersive transport, deposition and resuspension, mixing of surface sediments, scour and burial of subsurface sediment, partitioning to organic carbon on sediment particles, diffusive exchanges between the sediment and water column and between layers in the sediment bed, and volatilization. In addition, the temperature dependence and seasonal behavior of important processes such as partitioning, diffusion and volatilization are represented.
For the reach between DePere Dam and Green Bay, the FoxSim sediment bed is discretized into 96 units in the horizontal (Figure 3) and 2-cm layers in the vertical to provide accurate representation of sediment bed processes. Sediment PCB distribution in the Lower Fox River is characterized by a high degree of heterogeneity in both the horizontal and the vertical; therefore, the development of a fate and transport model with this level of spatial resolution in the sediment column was considered essential for proper evaluation of system response to alternative remediation options.
FoxSim model parameterization was developed based on the conceptual model constructed for the Lower Fox River based on a variety of site-specific datasets. Contaminant data were used to specify PCB depth profiles for each sediment unit, and modeled sediment transport processes were parameterized to be consistent with radioisotope data and other relevant site data. Average burial rate simulated by the model is approximately 0.5 cm y− 1, with variations on the order of 0.1–0.2 cm y− 1 for individual segments. Target burial rates were determined based on analysis of 137Cs profiles for cores collected downstream of DePere Dam.
The sediment mixed depth (z) was specified as 10 cm in FoxSim based on analysis of 137Cs and 210Pb profiles for these cores (LTI, 2002b). Mixing rates were varied from a maximum of 6 cm2 y−1 at the surface to a minimum of 0.6 cm2 y−1 at the bottom of the mixed layer. While there are several physical and biological factors that contribute to mixing in surface sediments, bioturbation is most often viewed as the dominant mixing process. The level of biologic activity in aquatic sediments is often sufficient to maintain nearly uniform concentrations over depths ranging from 3 to 10 cm (Boudreau, 1994; Thoms et al., 1995; Matisoff and Wang, 2000). In fresh water, oligochaetes (worms) are among the most important bioturbators; 90% of oligochaetes mix in the top 10–12 cm of sediments (McCall and Tevesz, 1982). The dominance of aquatic worms (class: Oligochaeta) and chironomid midges (order: Diptera) was confirmed in 1999 and in earlier studies of soft substrates of the Lower Fox River (Integrated Paper Services, 1990, 2000). Boudreau (1994, 1998) analyzed measurements of mixing depth for a large number of bioturbation sites and reported a ‘worldwide’ mean mixed depth of 9.8 cm with a standard deviation of 4.5 cm. These observations are consistent with the establishment of a 10 cm mixed layer depth in FoxSim.
Calibration of the Lower Fox River model was conducted for a 12-year period using extensive site-specific data collected during the 1989–2001 period. Major calibration datasets included sediment PCB surveys conducted in 1995 and 2001, water column total suspended solids (TSS) and PCB surveys conducted in 1989, 1995, and 2001, and lipid-normalized fish tissue PCB concentration data collected for 1971–2000. Key temporal trends for PCB concentrations in water, sediment, and biota were identified based on analysis of these datasets and used as long-term calibration targets. Primary model calibration parameters included empirical settling and resuspension coefficients and rates, and depth and rate of mixing in the sediment bed layers. Comparisons of model-predicted PCB concentration trends and data-based trends are shown in Figure 4. Water column TSS and PCB concentration data were also used to calibrate the short-term response of the model during a range of flow conditions, including high flow events, observed in the system during the 1989–2001 period.
The calibrated FoxSim model was applied to evaluate the relative effectiveness of a range of alternative remediation strategies, including MNR. Our experience in developing and applying FoxSim to evaluate these alternative options suggests that the thickness of the upper mixed layer is a key parameter controlling the rate at which a system will recover when subjected to natural processes. This stands to reason because the thickness of the upper mixed layer of sediments represents the reservoir of risk-generating contaminants subject to recovery processes. Whether the transport of contaminant out of the surface sediments is by burial, resuspension, or porewater mass transfer (Figure 1), the rate of the process can be expressed as a mass transfer velocity, u (e.g., cm y−1). Then if we divide that mass transfer rate by the thickness of the sediment layer, z (e.g., cm), to which it is being applied (the upper mixed layer), we get a specific chemical ‘washout’ rate (k = u/z, with units of y−1). The larger the value of k, the faster the recovery rate for chemicals in the upper mixed layer. For a given mass transfer rate, decreasing mixed layer thickness will produce faster recovery rates. The calibrated FoxSim model predicts ongoing natural recovery for the surface sediments (0–10 cm) in the reach downstream of DePere Dam with a half-time of approximately 20 years.
Model diagnostic approach
A diagnostic sensitivity analysis was conducted with the FoxSim model to investigate the impact of varying key sediment input and process parameters on forecasted recovery rates for surface sediment. Sediment mixing depth, initial PCB profile, and net sedimentation rate were included in the sensitivity analysis. It should be noted that the rate of mixing between sediment layers was also identified as a potentially important parameter; however, it was considered to be of less importance than the other three parameters, and therefore not included in the sensitivity analysis. A constant vertical mixing rate of 30 cm2 y−1 was applied between 2-cm sediment segments in the upper mixed layer for all diagnostic model simulations. This higher rate of mixing (relative to calibrated rates) was utilized to ensure that complete mixing occurred over the upper mixed depth, thereby eliminating the impact of mixing rate on the simulated recovery rate. The selected mixing rate of 30 cm2 y−1 is consistent with values reported by Boudreau (1998), Thoms et al. (1995), and Chan and Bentley (2004).
Mixing depths of 10, 20, and 30 cm were simulated, with the 10 cm case serving as the baseline. Net burial rates were modified by changing model settling parameters. Average burial rates of 0.10, 0.25, 0.50, and 1.00 cm y−1 were simulated, with the 0.50 cm y−1 case serving as the baseline. The sediment PCB profile was evaluated by simulating three different initial profiles with the FoxSim model (Figure 5). The baseline profile was based on the 1995 sediment contaminant data and characterized by concentrations of 4 ppm for the surface (0–10 cm) layer and 17 ppm for the 10–30 cm depth interval. A modified profile was constructed by specifying concentrations of 8 ppm and 25 ppm for the 10–20 cm and 20–30 cm layers, respectively. The average concentration for the 10–30 cm was maintained at 17 ppm, therefore providing a realistic alternative to the baseline profile developed using a single concentration for the 10–30 cm depth interval. Finally, a uniform profile was constructed by applying the 0–10 cm concentration to the 10–30 cm depth interval for each individual sediment unit. This profile was included for the purpose of eliminating the impact of vertical PCB concentration gradient in the analysis.
Model diagnostic results
Model simulations were conducted for a 50-year forecast period. Surface sediment recovery trends were grouped and analyzed for each of the three PCB profile cases. The results for the baseline PCB profile case are shown in Figure 6. The baseline forecast (10 cm mixed layer depth) demonstrates a consistent trend of recovery for the entire 50-year period, with a half-time of approximately 20 years. The projections for the 20-cm and 30-cm mixing depth cases are significantly different than the baseline case projection. The model predicts a more than 100% increase in average surface sediment PCB during the first 4–8 years of the forecast period for the two sensitivity simulations. Following the initial increase, sediment PCB concentrations decline over the remainder of the 50-year period; however, the rate of recovery is notably slower than that exhibited by the baseline case. Based on trend analyses, the 20-cm mixed depth simulation demonstrates recovery at a half-time of 75 years, while the 30 cm case demonstrates recovery at a half-time of approximately 120 years.
Recall that the baseline sediment PCB profile was characterized by a 0–10 cm concentration of 4 ppm and a 10–30 cm concentration of 17 ppm (Figure 5). For the baseline mixing depth of 10 cm, a downward trend in the average surface sediment PCB concentration is observed for the entire forecast period as a result of net burial of relatively cleaner sediments. For this case, PCB mass initially buried below 10 cm has little or no impact on surface sediment concentrations. The initial increases in average surface sediment PCB concentration observed for the 20-cm and 30-cm mixed depth cases occur as a result of PCB mass below 10 cm ‘bleeding’ into the 0–10 cm depth interval via deeper sediment mixing processes in the model. During the first 4–8 years of these simulations, the concentration gradient between the 0–10 cm and 10–20 cm layers is significant, and PCB mass is added to the top 10 cm by vertical mixing faster than it can be removed to deeper sediment via burial processes. At the end of this initial period of increase, a uniform concentration profile has been established over the 0–30 cm depth interval, burial becomes the dominant process, and surface sediments demonstrate a steady downward trend in PCB concentration over the remainder of the forecast period.
The results for the modified PCB profile case are compared to the baseline profile results in Figure 6. The forecast for the 10 cm mixed depth case is virtually identical to the baseline case because PCB mass initially in the 10–30 cm layer is not mixed into the 0–10 cm layer. The results for the 30-cm mixing depth are also similar to the results for the baseline; however, a longer period of time is required for the 0–30 cm interval to become completely mixed due to the alternate distribution of PCB mass (i.e., higher concentration/mass in the 20–30 cm layer than the 10–20 cm layer). The most significant difference between corresponding simulations for the baseline and modified PCB profiles is observed for the 20-cm mixed depth case. The modified profile simulation demonstrates a smaller initial increase and an overall faster recovery rate, with a half-time of 50 years compared to the 75-year half-time for the baseline profile, z = 20 cm case. Of course, this is because the modified profile has less PCB mass in the ‘active’ 0–20 cm interval than the baseline profile, and more PCB mass in the 20–30 cm interval that remains essentially inactive.
Results were also compared for the uniform PCB profile case; the forecast results for this set of simulations are shown in Figure 7. For this case, no initial PCB concentration increases are observed for the 20-cm and 30-cm mixing depth simulations. No gradient is present at the beginning of these simulations because the 0–30 cm depth interval is already completely mixed; therefore, only small quantities of PCB mass from below 10 cm are introduced to the surface layer during the first 5–10 years of the simulation period. The predicted recovery half-times for the 20 cm case (32 years) and the 30 cm case (40 years) are longer than that predicted for the 10 cm case (20 years). Based on the theory presented above for a uniform concentration profile, which does not include upstream loads or processes other than burial, we would expect a doubling of mixed layer thickness to double the recovery half-time. In our model application, doubling the thickness of the mixed layer only increases the recovery half-time by about 50%. This is because in the actual modeled system, watershed and upstream sources are non-zero and vertical mixing processes oppose the burial processes that are promoting recovery of the 0–10 cm surface sediment layer.
Clearly the above analysis has demonstrated that there is an important interaction between the depth of the upper mixed layer and the PCB concentration profile in the sediments being modeled. Figure 8 presents a summary of the recovery half-times for the three depth profiles and the three mixed layer depths tested with FoxSim. It is quite apparent that when the PCB concentration depth gradient is large, it becomes more important to determine how deeply the upper sediments are mixing. That is evident by comparing the response to deeper mixed depth for the uniform profile results with those for the modified profile. Systems that have been undergoing recovery from historically high chemical loads to the sediments will have a profile gradient similar to that observed in the Lower Fox River; hence, quantifying the mixed layer depth will be quite important.
Another important factor in determining the rate of natural recovery is the rate of sedimentation (accumulation of cleaner sediments after external load abatement has occurred). All of our previous analyses have been conducted using an average sedimentation rate of 0.5 cm y−1. If, however, we investigate the mixed depth sensitivity over a range of sediment burial rates (Figure 9), we can make two observations. First, the system displays greater sensitivity to changes in the sediment burial rate at deeper mixed layer depths. In reducing the sediment burial rate from 1.0 to 0.1 cm y− 1for the 10 cm mixed depth, we see a smaller increase in recovery half-time than for the case of a 30 cm mixed depth. This is because at a given burial rate, the percent contribution of ‘clean’ sediments to the upper mixed layer in a given year is a smaller fraction of the total reservoir for the 30 cm mixed depth. The second observation is that accurately determining the mixed layer depth is more important for smaller sediment burial rates. Again, the residence time of sediments in the upper mixed layer (i.e., the source of exposure in these systems) is directly proportional to mixed layer depth and indirectly proportional to sedimentation rate.
We should point out that the Lower Fox River below DePere Dam is a typical contaminated sediment site in that it is net depositional with increasing PCB concentrations in deeper layers of the sediments. However, the concept of mixed layer depth being an important factor to sediment recovery also applies to systems that are at dynamic equilibrium with respect to sedimentation or even net erosional. A system with zero net sediment burial is characterized by a deposition flux that is balanced by the resuspension flux. A system like this can still experience contaminant recovery in the upper mixed sediments if the sediments from upstream that are depositing are ‘cleaner’ than the near-surface sediments being resuspended. But of course, that process of contaminant dilution will occur much more slowly if deeper high concentration sediments are being mixed into the upper sediment-water exchange layer.
Finally, we should mention that the analysis presented in this paper does not consider a situation where concentrations at the surface are higher than concentrations at depth. This situation is rare for contaminated sediment sites; however, it can exist in areas that are net erosional (or at dynamic equilibrium) but are recovering because clean upstream sediments are mixing with historically contaminated, in-place sediments. For this type of initial profile, a larger mixing depth can actually ‘bleed’ contaminant mass into sub-surface sediment layers, resulting in a false initial decline in surface concentrations that is analogous to the false initial increase observed in Figure 6 when the initial profile is increasing with depth.
We have used an extensively calibrated PCB fate and transport model (FoxSim) for the reach of the Lower Fox River below DePere Dam to demonstrate that this system's natural recovery rate is inversely proportional to sediment mixing depth and very sensitive to model parameterization for burial rate, mixing depth, and initial PCB depth profile. Given that the Lower Fox River site is typical of contaminated sediment sites, these findings suggest that accurate characterization of mixing depth, burial rate, and contaminant depth profile is essential for producing valid forecasts for recovery. Therefore, it is critical to incorporate this information into a conceptual model of a site prior to, or concurrent with, development of a numerical model.
In developing a site conceptual model, a sampling program aimed at deriving mixing and burial parameters is recommended. However, one should be aware that these parameters are strongly scale-dependent. Therefore, the resolution of the contaminant depth profiling is also very important. In the absence of available data on mixing depth, burial rates, and contaminant profiles, the recommended sequence for data collection is as follows: 1) high resolution radioisotope-dated cores to estimate mixed depth and long-term net sediment burial rate (note that there are other methods for determining mixing depth in surface sediments); 2) contaminant depth profile sampling with vertical resolution based on results of radioisotope coring; and3) corroboration of mixing depths with (average) measured contaminant profiles. Sediment layers that are well-mixed are expected to have very similar contaminant levels.
We gratefully acknowledge the support of the Fox River Group, who funded studies on which this paper is based. We also thank colleagues at Blasland, Bouck and Lee, who conducted the 2000–2001 Fox River sampling program. Nevertheless, the views expressed here are our own, as is the responsibility for any errors or omissions.