1 Introduction

With the predicted climate change, it is expected that the chances of river flooding increase. During flood events, sediments will resuspend and when sediments are polluted, contaminants can be transferred to the surrounding water (Förstner 2005; Owens 2005) and further downstream whether dissolved or attached to particles to delta areas. Most of the research is focussed on the transport of particulates themselves (e.g., (Rovira and Ibanez 2007). Only recently, the role of sediments as secondary source of pollution is recognized (Wolz et al. 2009) Mass transfer of organic compounds like Persistent Organic Pollutants (POPs) from sediment particles to the surrounding aqueous phase, and vice versa, is affecting fate and transport of these chemicals in the aqueous environment. Mass transfer and mass transfer limitation of organic contaminants in polluted sediments and soil has been a key research issue for the last decades, and hundreds of articles have been published in this area of research (Reichenberg and Mayer 2006). Roughly, a distinction can be made between research focused on equilibrium sorption, and research focused on (de)sorption kinetics. Equilibrium sorption and sorption kinetics are both part of mass transfer phenomena and have a strong relation (Rulkens et al. 2004; Sabbah et al. 2005). An overview of different model concepts presently used to describe desorption kinetics was presented by Saffron et al. (Saffron et al. 2006). All but one of these models requires at least two (mathematical) compartments to fit the experimental data. These compartments can then define a combination of an instantaneous compartment where no mass transfer limitation is assumed with a dynamic compartment where mass transfer is limited (Ball and Roberts 1991; Sabbah et al. 2005), two dynamic compartments (Brusseau et al. 1991; Weber et al. 1992; Xing and Pignatello 1996; Cornelissen et al. 1998; Ghosh et al. 2000; Shor et al. 2003; Gamst et al. 2004) or a continuum of compartments with various parameters (Culver et al. 1997; Werth et al. 2000; Werth and Hansen 2002). Although sometimes excellent fits of experimental and modeled data were demonstrated in the different papers, the physical explanation of the desorption process and its limitations is only briefly elaborated and do not include the effect of particle-size distribution.

In diffusion models, the driving force of (de)sorption is related to a concentration gradient and a sorption concept. For example, Freundlich sorption isotherms (Miller and Pedit 1992; Lin et al. 1994; Rugner et al. 1999; Braida et al. 2001; Gamst et al. 2004) or linear sorption isotherm (Wu and Gschwend 1986; Steinberg et al. 1987; Rijnaarts et al. 1990; Ball and Roberts 1991; Brusseau et al. 1991; Pignatello et al. 1993; Pedit and Miller 1994; Li and Werth 2004) were used to model desorption kinetics. In radial diffusion models, particle-size effects are generally lumped into a single-fit-parameter Da/R2, where Da is a diffusion coefficient and R the radius of the spheres. Particle-size effects or more specifically, particle-size distribution effects are generally not included. (Wu and Gschwend 1986), however, specifically included measured particle sizes in their diffusion model and concluded that the radial diffusion model was the best model to fit their experimental data using artificial contaminated soils and sediments. They reported that large particles show a slower sorption approach to equilibrium than otherwise similar smaller particles when using the same sorbate. They concluded that sorption kinetics is controlled by intraparticle diffusion.

In this paper, we discuss a numerical intraparticle diffusion model similar to the model used by Wu and Gschwend that simulates desorption of dieldrin from a suspension of contaminated porous sediment particles with a well-characterized particle-size distribution. The experimental setup and the experimental data that are used in this paper were published before as a separate paper (Smit et al. 2008).The objective of this study is to model the experimentally obtained desorption rates (flux) of dieldrin. Similar to the experiments described in the experimental paper, we base our model on a suspension of field-aged sediment. We model the desorption process using different hydraulic retention times (HRT) of the aqueous phase. Furthermore, we elaborate the effect of particle-size distribution on mass transfer.

2 Materials and methods

2.1 Samples

Sediment samples were from Broekpolder (Vlaardingen, The Netherlands) and were taken from a depth of 0–0.5 m. Sediments were contaminated with dieldrin (5.0 mg/kg dry matter) for more than 40 years. Before experiments started, the natural samples were sieved using a sieve with hole openings of 2-mm. Then, a second fractionation (sieves with hole opening of 32 μm and 125 μm) was performed to obtain a narrower particle size distribution. The dieldrin concentration of this fraction was found to be 3.0 mg/kg dry matter. The particle-size distribution of the fraction we used in the experiments measured by a laser diffraction method, showed two peaks with a log-normal size distribution around the peaks. One peak was found at a particle diameter of 10 μm and one peak at 84 μm. Detailed information of the sample characteristics were described before (Smit et al. 2008).

2.2 Desorption in SPEED reactor

Desorption kinetics were measured using the SPEED reactor described before (Smit et al. 2008). In short, a Schott bottle (500 cm3) was used as continuously stirred tank reactor (CSTR, 400 rpm). Contaminated sediment and 0.01 M CaCl2 solution were mixed to obtain a homogeneous suspension. From the start of the experiment, water is pumped at a set flow rate (Q) from the reactor vessel through stainless steel filters (Supelco, 2 μm) to a glass column containing 3 g of Tenax where the dieldrin is extracted from the aqueous solution. Before reaching the Tenax particles, the water passed a glass filter that serves two functions: retaining Tenax in the column and filtering small particles that were possibly present in the water even after the first steel filter. After the Tenax column, the water was directed through a control Tenax column to validate that all dieldrin was indeed extracted and then recycled into the CSTR. At predetermined volume intervals, the loaded Tenax column was replaced by a clean column and analyzed for pesticides. All reactor parts were made from HDPE, glass, or stainless steel to minimize sorption of dieldrin. The reactor setup is presented in Fig. 1. The experimental results were published in a separate paper (Smit et al. 2008).

Fig. 1
figure 1

Conceptual model of mass transfer from organic particles in the SPEED reactor. The offset illustrates a pore in an organic particle. In our modeling approach, n = 61. Symbols are elaborated in the text

Full size image

2.3 Modeling concept

The model is made in such a way that the input and output can be compared with the experimental data we published before (Smit et al. 2008) The model is focused on the concentration of dieldrin in the aqueous phases. We distinguish a bulk aqueous phase and an aqueous phase in the pores of organic matter. The dieldrin concentration within the pores is considered to be variable in 1 dimension (1D diffusion). Organic matter was assumed to be present as separate, porous particles that are homogeneously distributed within the sediment sample. Each particle had identical physical/chemical properties except for the particle size. This seems reasonable as the fraction of organic matter and the dieldrin concentration were similar within the subsamples described in the experimental paper (Smit et al. 2008). Inorganic particles were not considered.

2.3.1 Mass balance

The volume of the organic particles in the reactor, \( {V_{{\text{om}}}}\left[ {{\text{cm}}_{{\text{particles}}}^3} \right] \), was calculated with Eq. 1.

$$ {V_{\text{om}}} = \frac{{{X_{\text{sed}}} \times {f_{\text{om}}}}}{{{\rho_{\text{om}}}}} $$
(1)

where X sed is the sediment mass [g], f om is the mass fraction organic particles [–], and ρ om is the density of the organic particles \( \left[ {{\text{g}} \times {\text{cm}}_{{\text{particles}}}^{ - 3}} \right] \). We assumed a density of \( 1.20\,{\text{g}} \times {\text{cm}}_{{\text{particles}}}^{ - 3} \). At any location r within the pores of organic particles, we assumed local equilibrium of dieldrin between the immobile phase (organic matter) and the mobile phase (pore water). The sorption isotherm for this equilibrium was linear:

$$ s(r) = {K_{\text{om}}} \times {c_{\text{aq}}}(r) $$
(2)

where s(r) is the dieldrin concentration of the organic matter at location r \( \left[ {{\text{g}} \times {\text{cm}}_{{\text{om}}}^{ - 3}} \right] \), K om the sorption coefficient of dieldrin to organic matter \( \left[ {{\text{cm}}_{{\text{pw}}}^3 \times {\text{cm}}_{{\text{om}}}^{ - 3}} \right] \), and c aq(r) the dieldrin concentration in the pore water at location r \( \left[ {{\text{g}} \times {\text{cm}}_{{\text{pw}}}^{ - 3}} \right] \). The sorption coefficient K om was calculated by log(Koc) = 4,46 divided by the density of organic matter resulting in \( 3.48 \times {10^4}\,{\text{cm}}_{{\text{pw}}}^3 \times {\text{cm}}_{{\text{om}}}^{ - 3} \). The concentration of dieldrin in the bulk liquid is assumed to be homogeneous and changes only with time. As dieldrin is reported to be very resistant to (bio)degradation and volatization, the overall amount of dieldrin present in the system remains constant. The mass balance equation of dieldrin in the SPEED reactor will be:

$$ \frac{{d{C_{\text{aq}}}(t)}}{{dt}} = \frac{1}{V} \times \left( {J(t) - {C_{\text{aq}}}(t) \times Q} \right) $$
(3)

where C aq(t) is the bulk aqueous dieldrin concentration [g × cm–3] at time t [s], V is the volume of the bulk aqueous phase in the SPEED reactor [400 cm3], and Q is the applied flow rate [cm3 × s–1]. The dieldrin mass flow rate J(t) from the organic particles to the bulk aqueous phase [g × s–1] is calculated as:

$$ J(t) = \frac{{dS(t)}}{{dt}} \times {V_{\text{om}}} $$
(4)

where S(t) is the average dieldrin concentration of the organic particles \( \left[ {{\text{g}} \times {\text{cm}}_{{\text{particle}}}^3} \right] \). Analytical solutions of Eq. 3 are given for two boundary conditions: without desorption J(t) = 0, Eq. 5:

$$ \begin{array}{*{20}{c}} {\frac{{{C_{\text{aq}}}(t)}}{{{C_{\text{aq}}}\left( {t = 0} \right)}} = {e^{\frac{{ - Q}}{V} \times t}}} & {\text{and}} & {\frac{{S(t)}}{{S\left( {t = 0} \right)}} = 1} \\ \end{array} $$
(5)

and for an instant equilibrium between organic matter and the bulk aqueous phase Eq. 6.

$$ \begin{array}{*{20}{c}} {\frac{{{C_{\text{aq}}}(t)}}{{{C_{\text{aq}}}\left( {t = 0} \right)}} = {e^{\frac{{ - Q}}{{V + {V_{\text{om}}} \times {K_{\text{om}}}}} \times t}}} & {\text{and}} & {\frac{{S(t)}}{{S\left( {t = 0} \right)}} = \frac{{{C_{\text{aq}}}(t)}}{{{C_{\text{aq}}}\left( {t = 0} \right)}}} \\ \end{array} $$
(6)

Experimental results of a control experiment without any particles confirmed the validity of Eq. 5 down to a dieldrin concentration \( \left( {{{\text{C}}_{\text{aq}}}/{{\text{C}}_{\text{aq}}}\left( {{\text{t}} = 0} \right)} \right) \approx 0.0{\text{5}} \) (data not shown). The solutions of Eqs. 5 and 6 are plotted with dashed lines in Fig. 2a, b as function of dimensionless time \( \theta = \frac{{Q \cdot t}}{V}\left[ - \right] \).

Fig. 2
figure 2

Calculated (monodisperse, solid lines) and experimentally measured (◊) normalized dieldrin concentration as function of dimensionless time (θ) in water (a) and sediment (b) for various values of particle diameter. Dashed lines are theoretical limits Eqs. 2 and 3. The thick solid line through the symbols was calculated using the bi-disperse model with particle diameters of 10 and 84 μm

Full size image

Particle-size distribution (PSD) The measured PSD of the sediments was presented before (Smit et al. 2008). The PSD showed two distinct peaks (k = 2). The total number of organic particles, N p,tot [–], is the sum of particles in each defined particle-size class N p,i [–]:

$$ {N_{{\text{p,tot}}}} = \sum\limits_{i = 1}^k {{N_{{\text{p,}}i}}} = \sum\limits_{i = 1}^k {\frac{{{f_{{\text{p,}}i}} \times {V_{\text{om}}}}}{{4/3\pi \times R_i^3}}} $$
(7)

where k is the number of different particle sizes with radius R i [cm] and f p,i is the volume fraction as calculated from the PSD. Both diameters represent a narrow PSD as described by Cooney et al. (Cooney et al. 1983). The volume fraction associated with the peak at 10 μm was 0.27 and for the peak at 84 μm was 0.73. These particle sizes and their corresponding volume fractions were the input of the bi-disperse particle size distribution.

2.3.2 Radial diffusion model

The mass flow rate J(t) of dieldrin from the organic particles to the bulk liquid was modeled assuming that mass transfer of dieldrin is only possible through the pore liquid and that local sorption equilibrium is instantaneous. Furthermore, we assumed that because of vigorous mixing of the slurry, mass transfer limitations only occurred within the particles (intraparticle diffusion). The local total volumetric dieldrin concentration within an organic particle at location r is defined as:

$$ S\prime(r) = \varepsilon \times {c_{\text{aq}}}(r) + \left( {1 - \varepsilon } \right) \times s(r) $$
(8)

where S’(r) is the local total volumetric dieldrin concentration \( \left[ {{\text{g}} \times {\text{cm}}_{{\text{particle}}}^{ - 3}} \right] \) and ε is the particle volumetric porosity \( \left[ {{\text{cm}}_{{\text{pw}}}^3 \times {\text{cm}}_{{\text{particle}}}^{ - 3}} \right] \). The particle volumetric porosity was assumed to be 0.4. The change of local total volumetric dieldrin concentration in time as function of the concentration gradient within the particle can then be stated as:

$$ \frac{{\delta S\prime(r)}}{{\delta t}} = \varepsilon \times \frac{{\delta {c_{aq}}(r)}}{{\delta t}} + \left( {1 - \varepsilon } \right) \times \frac{{\delta s(r)}}{{\delta t}} = \varepsilon \times \frac{{{D_{\text{aq}}}}}{\kappa } \times \left( {\frac{{{\delta^2}{c_{\text{aq}}}(r)}}{{\delta {r^2}}} + \frac{2}{r} \times \frac{{\delta {c_{\text{aq}}}(r)}}{{\delta r}}} \right) $$
(9)

where D aq/κ is the matrix diffusion coefficient of dieldrin in the aqueous phase corrected with tortuosity [cm2 × s–1]. D aq/κ was optimized by fitting all experimental data and minimizing the sum of squared differences between experimental data and modeled results. The optimized value (2.0·10–7 cm2 s–1) was then used for all calculations. Substitution of Eqs. 2 and 8 into Eq. 9 leads to:

$$ \frac{{\delta S\prime (r)}}{{\delta t}} = \frac{{\varepsilon \times \frac{{{D_{\text{aq}}}}}{\kappa }}}{{\varepsilon + \left( {1 - \varepsilon } \right) \times {K_{\text{om}}}}} \times \left( {\frac{{{\delta^2}S\prime (r)}}{{\delta {r^2}}} + \frac{2}{r} \times \frac{{\delta S\prime (r)}}{{\delta r}}} \right) $$
(10)

Initial and boundary conditions for Eq. 10 are as follow:

$$ \begin{array}{*{20}{c}} {\frac{{S\prime (r)}}{{\left( {\varepsilon + \left( {1 - \varepsilon } \right) \times \rho \times {K_{\text{om}}}} \right)}} = {C_{\text{aq}}}\left( {t = 0} \right)} \hfill & {\begin{array}{*{20}{c}} {{\text{for}}\,0 < {\text{r}} < {\text{R}}} & {\left( {{\text{IC}}\,1} \right)} \\ \end{array} } \hfill \\ {\frac{{S \prime \left( {r = R} \right)}}{{\left( {\varepsilon + \left( {1 - \varepsilon } \right) \times \rho \times {K_{\text{om}}}} \right)}} = {C_{\text{aq}}}(t)} \hfill & {\left( {{\text{BC}}\,1} \right)} \hfill \\ {{{\left. {\frac{{\delta s \prime (r)}}{{\delta r}}} \right|}_{r = 0}} = 0} \hfill & {\left( {{\text{BC}}\,2} \right)} \hfill \\ \end{array} $$

The mass flux of dieldrin from the organic particles to the bulk liquid will then be:

$$ J(t) = \frac{{dS(t)}}{{dt}} \times {V_{\text{om}}} = \sum\limits_{i = 1}^k {{N_{{\text{p,}}i}}} \times \frac{d}{{dt}}\int_0^{{R_i}} {4\pi {r^2} \times S\prime (r) \times dr} $$
(11)

We used a numerical integration method similar to the method described by Rügner (Rugner et al. 1999) to approximate the solution of Eqs. 3 and 11. A Crank Nicolson discretization scheme was used to transform the differential equation into a set of linear equations that can be solved according to LU decomposition. We used 61-space nodes for sufficient resolution of the intraparticle concentration of dieldrin (Basagaoglu et al. 2002). The system of equations was programmed in Matlab®. The parameters that were measured in our experiments (Smit et al. 2008) and that are used to model the different experimental conditions are given in Table 1.

Table 1 Input parameters for SPEED model calculations
Full size table

3 Results and discussion

3.1 Particle-size effect on desorption

In Fig. 2a (bulk aqueous phase) and b (sediment phase), dieldrin concentrations are presented as a function of dimensionless time (θ) for various values of the particle size. All calculations were performed using conditions similar to the experiments performed at HRT ≈ 11 min. Lines, labeled by the particle diameter, are the model results of that single-particle diameter. These model results were not validated by experimental data as our sediment did not have these properties. As can be seen in Fig. 2b, the dieldrin concentration of the sediment phase decreases faster when particles are smaller. For these smaller particles, this leads to a higher concentration in the bulk aqueous phase at a specific time compared with the larger particles (see Fig. 2a). The continuously decreasing concentration is the result of the applied dilution of the reactor content (1/HRT ∼ 0.1 min–1). Only when particles are almost depleted with dieldrin, the concentration in the aqueous becomes lower compared with the larger particles as can be seen for particles with a diameter of 10 μm. In Fig. 2b, we see that the boundary conditions calculated with Eqs. 5 and 6 are similar to the upper and lower limits of the soil particle size. Small particles with a diameter of 10 μm already show some mass transfer limitations and a reduced desorption rate compared with the boundary condition. The largest particles (2,000 μm) demonstrate a very slow release of dieldrin, and concentrations are about equal to the boundary condition where desorption is absent. Only at a very low normalized aqueous dieldrin concentration \( \left( {{C_{\text{aq}}}/{C_{\text{aq}}}\left( {t = 0} \right) \approx {\text{6}}.{\text{1}} \times {\text{1}}{0^{ - 4}}} \right) \) desorption starts to become visible. A single set of data points measured and described in the experimental paper (HRT ≈ 11 min) is plotted as diamonds. The calculated concentrations using the measured bi-disperse particle-size distribution with particles of 10 μm and 84 μm are plotted as the bold line in Fig. 2. From the start of the experiment, the behavior of the bi-disperse model shifts from small particles to larger particles, or in other words from rapid to slow desorption. This phenomenon could not be approached with a single, monodisperse particle-size distribution.

3.2 Comparison of modeling results with experimental data

In Fig. 3, we show the normalized concentration of dieldrin in the bulk aqueous phase (3A) and the sediment phase (3B) as function of desorption time for various values of hydraulic retention time (HRT). Particle size is considered to be constant in the experimental and modeled data sets. Experimental data was published before (Smit et al. 2008) and is depicted as points with various shapes for different HRT. The results of calculations with the bi-disperse particle-size distribution using the model that is elaborated in this paper are shown as lines. As can be seen in the figures, model results agree well with the experimental data. A smaller HRT, thus a higher dilution rate of the water, results in a faster decrease of bulk aqueous concentration in time. This is the result of diffusion limitation of dieldrin from the sediment to the bulk aqueous phase. Dieldrin that is removed from the bulk aqueous phase cannot instantly be replaced by dieldrin from the sediment phase, and therefore, the concentration in the bulk aqueous phase will be lower. These lower bulk aqueous dieldrin concentration enlarge the concentration difference between sediment and aqueous phase. As other properties like diffusion distance and diffusion coefficient remain constant, the desorption rate increases (1st law of Fick). Shorter HRT or higher dilution rate therefore results in a faster dieldrin release from the sediment in a specific time. Furthermore, we can see in Fig. 3a that, with increasing HRT, the curve linearity also increases. Increasing linearity indicates an approach to the concept of instant equilibrium. A higher HRT provides a longer time period available for mass transfer per liter of recycled water. Model results of the concentration of dieldrin in the sorbed phase (see Fig. 3b) tend to overestimate the extent of desorption at longer times. This overestimation might be related to the presence of larger particles in the experimental setup then we accounted for in the model but were found to be present in the measured PSD (Smit et al. 2008).

Fig. 3
figure 3

Calculated (solid lines) and experimentally measured (symbols) normalized dieldrin concentration as function of time in water (a) and sediment (b) for various values of HRT. Experimental results are duplicates and symbols represent: filled diamond, open diamond = 430 minutes, filled square, open square = 53 min, filled triangle, open triangle = 20 min, and filled circle, open circle = 11 min

Full size image

The calculated concentration gradients within particles are presented in Fig. 4. The normalized dieldrin concentration after different desorption times (labels) are given for small (4A and 4C) and large (4B and 4D) particles and for the highest (4C and 4D) and lowest (4A and 4B) HRT. Just before desorption commences, the concentration of dieldrin has the same value at any position in the particle. When desorption starts, the concentration of dieldrin at the particle/water interface decreases. Dieldrin concentration at the particle/water interface (5 μm in 4A and 4C, 42 μm in 4B and 4D) at a given time are similar for the small and large particles and are equal to the aqueous concentration presented by the graph in Fig. 3a for the same experiment (HRT 11 or HRT 430). The concentration gradient of dieldrin within small particles at HRT = 430 is very different from the concentration gradient within large particles. The apparent absence of a concentration gradient in small particles suggests equilibrium between these particles and the surrounding aqueous phase at any time, whereas the presence of such a gradient in large particles points at intraparticle diffusion limitations. This different behavior is even more pronounced at a low HRT (11 min.).

Fig. 4
figure 4

Normalized concentration of dieldrin within pores of small and large particles after different desorption times as function of the radial position at various values of time for HRT = 430 min and 11 min

Full size image

These results stretch the importance to combine traditionally separated disciplines utilized in water research as proposed in Wolz et al. (2009). The combination of hydrodynamic and ecotoxicological technologies will provide a fundament to unravel the true environmental impact of contaminated sediments and the effect of short-“pulse” exposure, which is characteristic for flood events.

4 Conclusions

Flood events, mostly characterized with turbulent flow conditions and high aqueous flow rates, will enhance the release of contaminants. The turbulent flow conditions can suspend sediment particles that were originally not in direct contact with the flowing water, and the high aqueous flow rate (low HRT) will refresh the water and increase the mass transfer of contaminants from particles to the water. In this study, we demonstrate that intraparticle diffusion limits desorption of hydrophobic contaminants from organic particles. Including the particle size and the particle size distribution in the radial diffusion model enabled us to successfully model our experimental data. The often-observed rapid and slow release of contaminants from sediments can be explained by the particle size distribution where small particles contribute mostly to the rapid release and large particles contribute mostly to the slow release. The time that is needed to remove all dieldrin from sediment particles varied from about 10 h for small particles (10 μm) to much more than 150 h for larger particles (84 μm). Our radial diffusion model combined with particle-size distribution therefore facilitates primarily the understanding and prediction of contaminant fluxes from sediments to the aqueous phase. The subsequent fate and transport of these contaminants once present in the aqueous phase is not covered with our model.

5 Recommendations and perspectives

To our opinion, many studies overlooked the importance of particle-size distribution and the time required to reach equilibrium. Currently, particle-size distribution is often reported by ISSS and NEN protocols; however, a more refined method is required to assess the particle-size distribution as small differences in particle size have a major effect on mass transfer rate (see Fig. 3). We think that the use of particle-size distributions contributes to the understanding of the phenomenology related to bioavailability in practice. The release of contaminants from sediment particles is merely the first step in a series of processes that ultimately leads to risks. Although we focused on the effect of aqueous-phase dilution on the release of dieldrin from sediment particles, other processes like uptake in dissolved organic matter, biological or chemical degradation will also affect the release rate. These processes share the behavior that the aqueous phase concentration of a contaminant is lowered enabling a further release of contaminants from sediment particles. Furthermore, environmental conditions like pH and redox might affect sediment properties like porosity and tortuosity and the distribution coefficient of contaminants between sorbed and dissolved states. A mechanistic approach, such as described in this paper, can contribute to unravel all these different aspects that are related to risk assessment and mitigation.