Table des matières

En juillet 1996 la région du Saguenay fut marquée par un déluge exceptionnel qui causa la crue des rivières et entraîna la déposition de plusieurs millions de mètres cubes de sédiments propres dans le fjord du Saguenay. Cette couche de turbidite, composée de sédiments propres représente une barrière contre la migration de métaux lourds et de HAP de la couche contaminée sous-jacente vers la nouvelle interface sédiments-eau. Un modèle numérique a été développé pour simuler la migration verticale et la remobilisation de contaminants dissous dans la couche de recouvrement naturelle. Le modèle représente les principaux processus physiques, chimiques et biologiques qui influencent le transport de soluté dans une colonne de sédiment. Le calage du modèle a été réalisée avec les profils de concentration de l'arsenic dissous dans les sédiments qui ont été mesurés à deux stations du fjord du Saguenay. Une analyse de sensibilité détaillée, basée sur la méthode du design factoriel, montre que les paramètres du modèle associés à la bio-irrigation ont un impact important sur la réponse du modèle et donc sur l’efficacité d’une couche de recouvrement.

In July 1996, two days of intense rainfall caused severe flooding in the Saguenay region and led to the discharge by rivers of several million cubic meters of clean sediments to the Saguenay Fjord. This turbidite layer, composed of clean sediments, represents a potential barrier for the migration of heavy metals and PAHs from the underlying contaminated sediments towards the new sediment-water interface. A numerical model has been developed to simulate the vertical migration and remobilization of dissolved contaminants in the capping layer. The model includes the main physical, chemical and biological factors affecting the transport of a dissolved compound in a sediment column. Calibration of the model has been achieved through comparison with the profiles of dissolved arsenic in the sediments, measured at two sampling stations of the Saguenay Fjord. A detailed sensitivity analysis, based on factorial design, shows that the model parameters associated with bio-irrigation have the greatest impact on the model output and by extension on the effectiveness of a capping layer.

The capping of contaminated sediments with clean material is considered as a promising alternative for the isolation of contaminants from the water column and the prevention of the diffusion of pollution through the ecosystem. This approach has the advantage of being applicable to large contaminated surfaces. A capping layer is designed to physically isolate the contaminants from the benthic fauna, stabilize the sediments, and prevent or at least reduce the contaminant flux towards the new water-sediment interface. Any attempt to predict the long-term effectiveness of the cap requires an understanding of the physical, chemical and biological factors affecting the migration of contaminants in the sediments.

The TRANSCAP-1D numerical model was developed to simulate the migration of dissolved contaminants through a subaqueous capping layer. The model considers the physical, chemical and biological processes affecting the fate of a dissolved contaminant in a sedimentary environment and can be used to predict the effectiveness of a sediment cap in isolating the contaminant from the water column. The transport equation includes advection, diffusion and the effect of bio-irrigating worms.

The model has been applied to the Saguenay Fjord, where a major flood led to the deposition of a natural capping layer over contaminated sediments. Since the catastrophic flood, which occurred in 1996, the geotechnical and geochemical characteristics of the capping layer, as well as the recolonisation of the benthic fauna, have been extensively studied in the context of a multidisciplinary research project. The ultimate objective of this research project is to assess if the new capping layer effectively isolates the contaminants from the water column. A large database was available for the model calibration and the response of the numerical model was successfully matched to concentrations of dissolved arsenic measured in the Saguenay Fjord. Arsenic is not a contaminant in the Saguenay Fjord, but the calibration of the model required profiles of the dissolved concentration in the sediment column, and those data were only available for this compound.

Sensitivity analysis is commonly used in engineering to assess the influence of a model input on the output. The simplest form of sensitivity analysis consists in increasing or decreasing the value of one input parameter at a time by some fraction, or perturbation, and then to evaluate the corresponding variation in the model output. This procedure is repeated for each parameter tested and the parameter causing the greatest modification in the model output is deemed the most important. Although this approach is simple and commonly used, there is no guideline for the selection of parameters to evaluate and for the magnitude of the perturbation to impose. In addition, since parameters are evaluated one by one, any significant combined effect of two or more parameters cannot be detected.

Kabala (2001) presents an overview of sensitivity analysis, with an application to pumping tests in aquifers. The definition of sensitivity of model output to input parameters given by Kabala (2001) is more precise than the general definition given above, because it is expressed in a mathematical form as the differential in the model output with respect to selected input parameters. By using this definition, sensitivities can be calculated with more precise mathematical techniques, such as the direct method that requires differentiating the governing equations or the adjoint method that requires solving the governing equation and its mathematical adjoint.

Other approaches have been developed in chemical engineering to determine the combined effect of several parameters on model output. An example is a two-level factorial design, which is used to determine the influence of n parameters on some observed value with a minimum number of experiments (Box et al. 1978). Two levels, or values, are selected for each parameter and a total of 2n experiments are sufficient to test the impact of each parameter, as well as every possible combination of parameters, on the experimental output. By limiting the number of trials and still considering all combined effects, the method is well suited for the rapid identification of the main parameters for the model, which could then be studied or characterized in more detail. The application of this method to the sensitivity analysis of numerical models has been very limited, especially in the environmental domain, although it is well suited for that purpose. A good example of the application of the factorial design to simulation results is presented by Elliott et al. (2000), where the method was used to determine the kinetic parameters with a significant effect on the efficiency of an anaerobic treatment process.

A large variability and uncertainty associated to several physical, chemical and biological parameters was observed for the Saguenay Fjord. This situation requires a systematic approach in order to identify the dominant parameters affecting the fate of contaminants in sediments. To investigate the sensitivities of the model parameters that are either variable or uncertain, a factorial design approach is presented here for the model simulating the migration of contaminants in sediments of the Saguenay Fjord. The purpose of the sensitivity analysis presented here is to deepen our understanding of the effect of variable input values on the model output and to identify the factors that have the greatest influence on the contaminant migration towards the new sediment-water interface.

The TRANSCAP-1D model simulates advection, diffusion, chemical reactions and the effect of bio-irrigation . The consolidation of the sediments and associated fluid advection are not included. The sediments are represented as a dual porosity medium composed of sediment pores and bio-irrigated tubes. The mathematical formulation consists of two partial differential equations, corresponding to the contaminant transport in the two media. These equations are coupled via a non-local exchange term that represents the mass transfer of contaminant between the sediments and the tubes . The formulation is mathematically similar to the one describing solute transport in a dual-porosity medium, in the context of groundwater flow (for example, Zheng and Bennett 2002).

The equation describing solute transport in the sediment porewater is given by:


where CS and CT are the solute concentration [M L-3] in the sediments and the tubes, respectively, ns is the sediment porosity [-], vs is the fluid velocity in the sediments [L T-1], Rt is the solute retardation factor [-], and β is a first-order mass transfer coefficient [T-1]. The dispersion coefficient of the solute in the sediments, Ds [L2 T-1], is given by:


where Dd is the effective diffusion coefficient [L2 T-1] and αL is the sediment dispersivity [L]. Note that when the fluid velocity is small in the sediments the dispersion coefficient is approximately equal to the diffusion coefficient.

The equation describing solute transport in the tubes is given by:


where nT is the porosity of the tubes [-], Dm is the dispersion coefficient for the tubes [L2 T-1], given by an expression similar to equation (3.2), and vT is the fluid or irrigation velocity in the tubes [L T-1]. In equation (3.3), S is a general source or sink term [M L-3 T-1] representing the release of contaminant from mineral dissolution associated to bio-irrigation. This term is defined later in this section.

The mass transfer coefficient β and the porosity of the tubes nT decrease exponentially with depth, following the distribution of bio-irrigated tubes (Martin and Banta 1992). The mass transfer coefficient β can be split in two coefficients β1 [T-1] and β2 [L-1], representing the maximal value of the mass transfer coefficient at the surface of the sediments and the coefficient of exponential decrease respectively:


As illustrated by Boudreau (1997), the value of the non-local term β1 [T-1] at the surface of the sediments can be estimated using the parameters of Aller’s cylindrical diffusion model (Aller 1980):


where r1 is the inner radius of the tubes/burrows [L], r2 is the half-distance between two tubes/burrows [L] and δ is the distance [L] from the burrow axis to a point where the concentration equals the horizontally integrated value.

The porosity of the tubes n°T at the surface is calculated using cylinders to approximate their geometry:


where m represents the quantity of tubes per unit area [L-2]. The exponential decrease with depth of the tube porosity can be represented with the same coefficient β2 [L-1]:


The source term S of the equation for the tubes (3.3) represents a chemical reaction responsible for the remobilization of contaminants. The representation of the chemistry is simplified and the model only considers the dissolution of trace metals, originally adsorbed or coprecipitated with iron-monosulfides (FeS). The model was calibrated with profiles of dissolved arsenic (As), a compound that was extensively studied in the Saguenay Fjord (Mucci et al. 2000a, 2000b, Salunier and Mucci 2000). In the anoxic zone, As is associated to FeS, but the oxidation of the Fe-sulfides in the bio-irrigated tubes, can lead to the remobilization of this compound. Other trace metals coprecipitate with Fe-sulfides (Ni, Cu, Co, Hg, Pb) and could possibly be affected by the rapid exposure to an oxic environment. However, compared to As, the trace metals that contaminate the Fjord sediments (Hg, Pb, Zn, Cu) are probably less mobile, as they are more likely associated to less reactive mineral phases (e.g. organic matter).

The dissolution mechanism represented in the model is a consequence of the exposure of the anoxic sediments, containing sulfides, to oxygenated water advected in the tubes by bio-irrigation. As observed by Furukawa (2001) the irrigated burrow walls are the interface between oxic-anoxic media and are characterized by steep geochemical gradients and rapid chemical mass transfer. In the model, the dissolution term is composed of a kinetic factor γ [L-1] and a coefficient L [ML-3], representing the maximal dissolved contaminant concentration for a given FeS concentration in the sediments. The coefficient L is derived from the linear correlation between the released As and initial solid FeS that was observed during resuspension tests of Saguenay sediments (Saulnier and Mucci 2000). The source term can be expanded in the following way:



where k1 and k2 are regression parameters that can be viewed as dissolution coefficients and CFeS [ML-3] is the concentration of FeS.

Initial and boundary conditions are required to solve both equations (3.1) and (3.3). To impose the initial conditions one specifies the initial solute concentration in the two domains, sediments and tubes. Boundary conditions are required at the top and bottom of the domain and can be either a prescribed concentration or prescribed mass flux for both domains.

The two governing equations (3.1) and (3.3) are discretized with the finite volume method (Versteeg and Malalasekera 1995). The equations are coupled through the mass exchange term and a simultaneous solution for the concentration in the sediments and tube is obtained for each finite volume. The system of equations is thus solved in a fully-coupled fashion, avoiding the use of iteration that is necessary when the equations are decoupled. The block tridiagonal matrix resulting from the discretization and assembly of the terms is solved using Thomas algorithm adapted to block matrices. The numerical solution was tested by comparing to the semi-analytical solution of Neville et al. (2000), which has been developed for one-dimensional solute transport in a dual porosity medium with multiple non-equilibrium processes. The numerical model reproduced almost perfectly the concentrations computed with the analytical solution (results not presented).

The first step of the two-level factorial design is to determine the number of factors, or parameters, that have to be tested and define the two levels, which are the minimum and maximum value for each input parameter. Each factor is assumed independent of the other factors, but a factor can be composed of a set of dependent input values. For n factors at two levels each, a total of 2n simulations are necessary for a full factorial design.

The choice of the parameters to be tested is based on the available information about the processes affecting the migration of contaminants in the sediments of the Saguenay Fjord. At the study site the sediments show a low permeability but are intensively bio-irrigated. Thus bio-irrigation, dissolution and retardation are supposed to influence the migration of contaminant. Consequently, a factorial design with five factors was selected: three factors related to bio-irrigation (number of tubes per square meter, depth of the tubes and bio-irrigation velocity) one factor related to the reactions (dissolution coefficients) and the retardation factor. The fluid velocity in the sediments vs was not considered for the sensitivity analysis because in the sediments of the Saguenay Fjord there is no evidence for such advection. This parameter was therefore set equal to zero for all simulations.

The minimum and maximum values of each tested parameter were determined from the available information, which includes data measured at the Saguenay Fjord and data derived from the literature. The values of the parameters are discussed in the following paragraphs and are presented in Table 3.1.

The depth of bio-irrigation is represented in the model by two parameters: the maximum depth attained by the tubes, biox, and the corresponding coefficient β2, where β2 describes the exponential decrease of both the mass transfer between tubes and sediments β and the tube porosity nT. The two parameters, biox and β2, are related since the maximum depth corresponds to the location where the mass transfer β and the tube porosity nT exponentially tend to zero. According to field observations and to measured bioturbation intensities (De Montety et al. 2000), the minimal value of the bio-irrigation depth biox was set to 0.06 m and the maximal value to 0.26 m. This corresponds to values of β2 equal to 70 m-1 and 16 m-1, respectively.

The number and the dimension of the tubes in the upper sediment layer are an indicator of the intensity of bio-irrigation in this zone. This factor affects two input values of the model: the tube porosity at the surface nºT and the mass transfer coefficient β1. To calculate these parameters, one must first set the minimum and maximum values of the inner radius of the tubes r1, as well as the number of tubes per square meter, m. The value of m also defines the half-distance between two tubes r2 according to:


The value of r1 results from observation on sediments sampled with a box-corer in the Saguenay Fjord during summer 1999 and 2000 and varies between 0.0005 m and 0.001 m. The minimum and maximum values of the number of tubes m were visually determined from photographs of the undisturbed bottom sediments, taken during the summer of 2001, and are equal to 500 and 1500 tubes per square meter, respectively. Thus, the mass transfer coefficient β1 may vary between 0.0157 d-1 and 0.119 d-1, whereas the tube porosity nºT varies between 0.00039 and 0.00471.

The value for the retardation factor, Rt, has been estimated from the model calibration on two measured arsenic concentration profiles and taking into account the average of the range of values presented by Fuller (1978) for arsenic. The minimum and maximum values of R have thus been set to 5 and 45, respectively.

The values for the bio-irrigation velocity vT were derived from data measured during laboratory tests carried out on polychaete families living in shallow water (Riisgard 1989, Riisgard 1991). In the calibrated model, the irrigation velocity is equal to 1 m d-1. The minimum and maximum values for the factorial design were determined by modifying this value by a factor of five. Thus, we obtain a minimum bio-irrigation velocity of 0.2 m d-1 and a maximum value of 5 m d-1.

The minimum and maximum values of the dissolution coefficients k1 and k2 were derived from the linear correlation between the dissolved As and solid FeS measured at two sampling stations of the Saguenay Fjord in 1998 (Saulnier and Mucci 2000). To define the levels of these parameters, it has been assumed that the regression coefficients, determined by means of the measured data, correspond to the minimum and maximum values of the dissolution coefficients. Thus, we assume that k1 varies between 1 and 1.7 and that k2 varies between 1 and 1.3.

The input parameters that were not tested in the present sensitivity analysis were set to constant values according to the calibrated model (Dueri et al. in press). Thus, the sediment porosity nS was set to 0.7, whereas the effective diffusion coefficient Dd and the molecular diffusion coefficient Dm were set to 1.6 x 10-5 m2 d-1 and 2.5 x 10-5 m2 d-1, respectively. Moreover, the initial concentration profile of contaminant corresponds to the one used for calibration, which was extrapolated from the arsenic concentration measured in the sediments of the Saguenay Fjord in 1998.

Four series of simulations were carried out for the sensitivity analysis (Table 3.2). The first two series of simulations considered five variable factors, including a variable dissolution coefficient, and were thus composed of 25 = 32 simulations. These simulations use dissolution coefficients determined for arsenic in the Saguenay sediments and thus represent a compound with a higher remobilization potential compared to the inorganic contaminants of the Fjord. Series 3 and 4 simulated the response of the model for a non-reactive contaminant and the dissolution coefficient was therefore equal to zero. This means that there were only four factors to test and the factorial design was composed of 24 = 16 simulations. To evaluate the influence of the thickness of the capping layer, two of the series were performed with a cap thickness of 0.2 m (series 1 and 3) and the other two with a cap thickness of 0.1 m (series 2 and 4).

Each series was composed of simulations representing every possible combination of the minimal and maximal values of each factor (Table 3.3). The results of the simulation of each series were sorted in a standard order and thereafter the Yates algorithm was applied in order to calculate the effect of each factor and of the combination of factors (Box et al. 1978). Although these effects are calculated from the results of simulations, which are values with a physical meaning, like fluxes or concentrations, after the application of the algorithm they represent the significance of the tested factors and are no longer directly linked to physical parameters. The effects can also be described as the change in the response as we move from the maximum value to the minimum value of a factor, and thus the significance of a factor or a combination of factors for the simulation results.

The sensitivity analysis was performed over a simulation time of 10 years after the capping event and considers two types of model output:

  1. The total cumulative mass that leaves the sediment column from the upper boundary and reaches the water column during the simulation.

  2. The final concentration at the top of the sediment column, to which the benthic fauna may be exposed.

The results of the sensitivity analysis were represented in a normal probability plot. This type of plot is generally used to test the normality of a distribution, but in our case it allows a straightforward identification of the factors with a significant influence on the output of the model. Each point on the graph corresponds to the effect of one or a combination of factors. The effects do not have units since they represent the average variation of the response for a changing factor over all conditions of the other factors (Box et al. 1978). The points that fall on or near the straight line are normally distributed and thus these effects cannot be distinguished from the normally distributed error. On the opposite, the effects that deviate from the straight line and therefore from the normal error distribution point out significant factors. The sum of squares indicates the deviation of the effect from the normal distribution and thus this value can be used to quantify the importance of a factor for the response.

The result of the first set of simulations, representing a reactive compound migrating through a cap of 0.2 m, is displayed in Figure 3.1. The graph illustrates that the factors with the greatest influence on the released mass are those associated to bio-irrigation (irrigation velocity V, number and dimension of tubes B2 and depth of the tubes B1). There is also a second group of factors affecting the output of the model. This group is composed of the combination of bio-irrigation factors with the factor of dissolution D. As the second group of effects shows a smaller deviation from the normal distribution, represented by the straight line, the influence on the model response is less significant.

Figure 3.1: Effect of tested factors on the mass release of a reactive compound. Cap thickness of 0.2 m. Significant effects are identified by the variable name.

The result of the second series of simulations, representing a reactive compound migrating through a cap of 0.1 m, shows a very similar trend to the previous series. The similarity is best illustrated by comparing the sums of squares, presented in Table 3.4.

Figure 3.2 shows the results of the third series of simulations, representing a non-reactive compound migrating through sediments covered by a cap of 0.2 m. In this case the release is strongly affected by the bio-irrigation depth B1, followed by the combination of bio-irrigation factors B1, B2 and V. The retardation factor Rt also exhibits a significant deviation from the normally distributed values. Figure 3.3 illustrates the results of series 4, representing the effects for a non-reactive compound migrating through a cap of 0.1 m. In this case, the depth of the tubes B1 is the only important factor identified on the graph.

Figure 3.2: Effect of tested factors on the mass release of a non-reactive compound. Cap thickness of 0.2 m. Significant effects are identified by the variable name.

Figure 3.3: Effect of tested factors on the mass release of a non-reactive compound. Cap thickness of 0.1 m. Significant effect is identified by the variable name.

In this section, the sensitivity analysis was used to determine the significant factors affecting the concentration at the top of the sediment column from the four series of simulations described previously. Figure 3.4 shows the result of the first series of simulations, considering a reactive compound and a cap thickness of 0.2 m. The graph shows the significant influence of the dissolution coefficients D, the bio-irrigation depth B1, the number and dimensions of the tubes B2 as well as the retardation factor Rt.

Shown in Figure 3.5 are the results of the factorial design for a reactive compound migrating through a cap of 0.1 m. In this set of simulations the most significant factor is the bio-irrigation depth B1, followed by the dissolution coefficients D. The combination of the depth and the number and dimensions of the tubes B1B2 is also an important factor, whereas the retardation factor Rt does not play a significant role.

Figure 3.4: Effect of tested factors on the surface concentration of a reactive contaminant. Cap thickness of 0.2 m. Significant effects are listed.

Figure 3.5: Effect of tested factors on the surface concentration of a reactive contaminant. Cap thickness of 0.1 m. Significant effects are listed.

The series representing a non-reactive compound are illustrated in Figures 3.6 and 3.7. The graphs do not display any outliers and none of the factors shows a significant deviation from the normally distributed population, represented by the straight line. Thus, they do not reveal any significant factor affecting the concentration in the surface layer.

Table 3.5 summarizes the significant factors highlighted by the sensitivity analysis. The second and third columns define the variable characteristics of the physical system (cap thickness) and of the contaminant (possibility of dissolution). The fourth and the fifth columns review the significant factors determined by the factorial design, depending on the considered output variable, i.e., the released mass or the surface concentration.

Figure 3.6: Effect of tested factors on the surface concentration of a non-reactive contaminant. Cap thickness of 0.2 m.

Figure 3.7: Effects of tested factors on the surface concentration of a non-reactive contaminant. Cap thickness of 0.1 m.

For a reactive compound (Series 1 and 2), the mass released over a period of 10 years is mainly affected by bio-irrigation factors (B1, B2, V). The cap thickness variation does not lead to relevant differences in the results of these two series. The dissolution factor D is not listed between the significant factors, thus the variations of the dissolution coefficients tested in this study do not have a significant effect on this type of response. However, if we consider a non-reactive contaminant (Series 3 and 4) and set the dissolution coefficients to zero, the result of the factorial design are different from those of the first and second series. In this case the influence of the cap thickness is evident and for a different thickness we obtain different significant factors.

The results for the surface concentration 10 years after capping show that in the case of a reactive compound, the cap thickness slightly influences the results of the factorial design. For a cap thickness of 0.2 m, the most important factor is the dissolution coefficient D, whereas if the cap is only 0.1 m thick, the depth of bio-irrigation B1 becomes the most important factor. In both cases, the bio-irrigation velocity V does not affect the results. Comparing the results of both output variables, one observes that for a reactive contaminant the bio-irrigation velocity V dominates the mass release, whereas the surface concentration is regulated by the dissolution coefficients D. The model is generally very sensitive to bio-irrigation parameters, except for the surface concentration of a non-reactive contaminant.

Factorial design has been used to determine the most significant parameters for the theoretical transport of reactive or non-reactive inorganic contaminants (trace metals) through the natural capping layer of the Saguenay Fjord. Two types of responses were considered for the sensitivity analysis: the total mass released from the system to the water column in the first 10 years after capping and the final concentration at the top of the sediment column. Both responses represent a potential harm for the aquatic fauna, since they indicate an exposure pathway for the organisms living in the marine environment. The results of the study pointed out that for the simulated conditions the significant parameters are a function of the response variable, the contaminant type (reactive or non-reactive) and to a lesser extent the cap thickness. The sensitivity analysis showed that the bio-irrigation parameters are generally very important for the prediction of the released mass, especially for reactive compounds that are associated to sulfides under anoxic condition and can be released if exposed to bio-irrigated water. It also demonstrated that the concentration of reactive compounds in the upper sediment layer depends on the bio-irrigation parameters as well as the dissolution coefficient.

The present study identifies the input values with a significant effect on the result, thus highlighting the parameters that should be studied carefully in order to reduce the uncertainty and get a better estimate of the performance of the capping layer. However, the most accurate characterization will not be able to eliminate the variability of the input values. Therefore, an uncertainty analysis will be performed to determine the distribution of the response considering the variability of the input parameters.

Aller, R. C. 1980. Quantifying Solute Distributions in the Biotubated Zone of Marine Sediments by Defining an Average Microenvironment. Geochim. Cosmochim. Acta. 44 : 1955-1965.

Box, G. E. P., Hunter, W. G., Hunter, J. S. 1978. Statistics for Experimenters: an Introduction to Design, Data Analysis, and Model Building. John Wiley, New York.

Boudreau, B. P. 1997. Diagenetic Models and Their Implementation: Modelling Transport and Reactions in Aquatic Sediments. Springer-Verlag, Berlin Heidelberg.

De Montety, L., Long, B., Desrosiers, G., Crémer, J.-F. and Locat, J. 2000.Quantification des structures biogènes en fonction d’un gradient de perturbation dans la baie des Ha! Ha! à l’aide de la tomodensitométrie axiale. Proceedings of the 53th Canadian Geotechnical Conference, Montréal, 15-18 Oct. 2000, Vol. 1, pp. 131-135.

Dueri, S., Therrien, R. and Locat, J. (In press). Numerical modeling of the migration of dissolved contaminants through a subaqueous capping layer. J. Environ. Eng. Sci.

Elliott, M., Zheng, Y. and Bagley, D. M. 2000. Determining Significant Anaerobic Kinetic Parameters Using Simulation. Environ. Technol. 21 : 1181-1191.

Fuller, W. H. 1978. Investigation of Landfill Leachate Pollutant Attenuation by Soils, US EPA, Municipal Environmental Research Laboratory, Cincinnati, OH.

Furukawa, Y., Bentley, S. J. and Lavoie, D. L. 2001. Bioirrigation Modeling in Experimental Benthic Mesocosms. J. Mar. Res. 59 : 417-452.

Kabala, Z. J. 2001. Sensitivity Analysis of a Pumping Test on a Well with Wellbore Storage and Skin. Adv. Wat. Res. 24 : 483-504.

Martin, W. R. and Banta, G. T. 1992. The Measurement of Sediment Irrigation Rates :A Comparison of the BR- Tracer and 222RN/226RA Disequilibrium Techniques. J. Mar. Res. 50 : 125-154.

Mucci, A., Guignard, C. and Olejczuyk, P. 2000a. Mobility of metals and As in sediments following a large scale episodic sedimentation event. Proceedings of the 53th Canadian Geotechnical Conference, Montréal, 15-18 Oct. 2000, Vol. 1, pp. 169-175.

Mucci, A., Richard, L.-F., Lucotte, M. and Guignard, C. 2000b. The differential geochemical behaviour of arsenic and phosphorous in the water column and sediments of the Saguenay fjord estuary, Canada. Aquatic Geochem. 6 : 293-324.

Neville C. J., Ibaraki, M. and Sudicky, E. A. 2000. Solute Transport with Multiprocess Nonequilibrium: a Semi-Analytical Solution Approach. J. Cont. Hydrol. 44 : 141-159.

Riisgard, H. U. 1989. Properties and Energy Cost of the Muscular Piston Pump in the Suspension Feeding Polychaete Chaetopterus Variopedatus. Mar. Ecol. Progr. Ser., 56 : 157-168.

Riisgard, H. U. 1991. Suspension Feeding in the Polychaete Nereis Diversicolor. Mar. Ecol. Progr. Ser. 70 : 29-37.

Saulnier, I. and Mucci, A. 2000. Trace Metal Remobilization Following the Resuspension of Estuarine Sediments : Saguenay Fjord, Canada. Appl. Geochem. 15 : 191-210.

Versteeg, H.K. and Malalasekera, W. 1995. An introduction to computational fluid dynamics: the finite volume method. Longman Scientific & Technical, Burnt Mill, Harlow, Essex.

Zheng, C. and Bennett, G. D. 2002. Applied Contaminant Transport Modeling. 2nd Edition, John Wiley, New York.