Mass Transfer Principles in Column Percolation Tests: Initial Conditions and Tailing in Heterogeneous Materials

Initial conditions (pre-equilibrium or after the first flooding of the column), mass transfer mechanisms and sample composition (heterogeneity) have a strong impact on leaching of less and strongly sorbing compounds in column percolation tests. Mechanistic models as used in this study provide the necessary insight to understand the complexity of column leaching tests especially when heterogeneous samples are concerned. By means of numerical experiments, we illustrate the initial concentration distribution inside the column after the first flooding and how this impacts leaching concentrations. Steep concentration gradients close to the outlet of the column have to be expected for small distribution coefficients (Kd<1 L kg−1) and longitudinal dispersion leads to smaller initial concentrations than expected under equilibrium conditions. In order to elucidate the impact of different mass transfer mechanisms, film diffusion across an external aqueous boundary layer (first order kinetics, FD) and intraparticle pore diffusion (IPD) are considered. The results show that IPD results in slow desorption kinetics due to retarded transport within the tortuous intragranular pores. Non-linear sorption has not much of an effect if compared to Kd values calculated for the appropriate concentration range (e.g., the initial equilibrium concentration). Sample heterogeneity in terms of grain size and different fractions of sorptive particles in the sample have a strong impact on leaching curves. A small fraction (<1%) of strongly sorbing particles (high Kd) carrying the contaminant may lead to very slow desorption rates (because of less surface area)—especially if mass release is limited by IPD—and thus non-equilibrium. In contrast, mixtures of less sorbing fine material (“labile” contamination with low Kd), with a small fraction of coarse particles carrying the contaminant leads to leaching close to or at equilibrium showing a step-wise concentration decline in the column effluent.

Although most leaching test procedures aim at equilibrium conditions in the column before the leaching starts, the true concentration distribution before the start of the percolation depends not only on the test procedure (contact time, pre-equilibration time, flow Materials 2021, 14 velocity during first flooding) but also on the properties of both the solid material and the pollutant of interest [26]. Equilibration and long-term leaching are further complicated if the test sample consists of a heterogeneous mixture of different material types and grain sizes [26][27][28], which is the common case if waste materials such as demolition waste or soils are tested. Finkel and Grathwohl (2017) evaluated the role of initial conditions for column leaching tests with intraparticle pore diffusion models by comparing the hypothetical scenarios, a perfectly equilibrated column vs. a column that wasn't equilibrated at all [26]. They could show that in many practical cases, peak and cumulative leachate concentrations are rather independent of the initial conditions. However, if release kinetics are slow due to large grain size or small intragranular porosity, the sensitivity to initial conditions is relevant, in particular for initial peak and early cumulative leachate concentrations.
The shortcomings of all previous studies [26,29,30], is that only uniform initial concentrations in the columns were considered in the leaching models. However, due to the specific conditions during the flooding of the column filled with initially dry material, the true initial conditions at the start of the leaching test may considerably deviate from this ideal, i.e., uniform distribution.
Against this background, the objectives of this study are (i) to illustrate the possibly non-uniform initial conditions that may be achieved after the first flooding of the column, (ii) to show the impact of these initial conditions on the temporal development of the effluent concentrations, and (iii) to investigate how heterogeneous mixtures of particles having different properties affect both the initial conditions in the column and the leaching of the solutes. To achieve that, we used numerical solutions for flow and transport in a column coupled to two kinetic models: (i) solute diffusion through an aqueous boundary layer and (ii) intraparticle pore diffusion. The implementation of the numerical models is described in detail in the appendices.

Local Equilibrium: The Advection-Dispersion Equation
In order to facilitate the understanding of mass transfer-limited cases of contaminant release in a column, we briefly introduce the equilibrium case for which the advectiondispersion model is commonly used: where is the dry bulk density of the packed bed in the column (porous media; ρ s is the solids density). For local equilibrium conditions the concentration in the solid phase (C s ) is in equilibrium with the solute concentration in water (C w ) and the distribution coefficient K d (= C s /C w ) allowing for the calculation of the respective concentrations. Under these conditions, Equation (1) can be simplified as: where R d [-] represents the retardation factor, defined as: Assuming equilibrium conditions and an initially uniform distribution of the solute in the column, leaching may be described by the analytical solution of Equation (2) [31]: where erfc denotes the complementary error function. The aqueous concentration at equilibrium, C w,eq , can be calculated from the initial solid concentration (C s,ini ) accounting for the mass balance when the contaminant mass in the water is equilibrated with the mass in the solid: The ratio n/ρ b [L 3 M −1 ] equals the liquid to solid ratio within the column, which in most cases is much smaller than in a batch leaching test (e.g., 0.25 L kg −1 for a column test with a porosity of n = 0.4 and a solid density of ρ s = 2.65 g cm −3 , compared to e.g., 10 L kg −1 in a batch test). Since leaching tests start for practical reasons with material packed more or less dry into the column, a uniform initial concentration is not necessarily achieved during the first flooding of the column. Initial conditions as assumed in Equation (4) (uniform concentration distribution), would only be achieved if the material is first mixed with water, equilibrated and then packed into the column, which is not practical. During the first flooding of the column, especially less sorbing solutes are displaced from the inlet and higher concentrations occur towards the outlet, as illustrated in Figure 1 (see also Appendix E). This may be accounted for by subtracting the distance of the solute displaced initially (x/R d with R d > 1) in Equation (4): where erfc denotes the complementary error function. The aqueous concentration at equilibrium, , , can be calculated from the initial solid concentration ( , ) accounting for the mass balance when the contaminant mass in the water is equilibrated with the mass in the solid: The ratio / [L 3 M −1 ] equals the liquid to solid ratio within the column, which in most cases is much smaller than in a batch leaching test (e.g., 0.25 L kg −1 for a column test with a porosity of = 0.4 and a solid density of = 2.65 g cm −3 , compared to e.g., 10 L kg −1 in a batch test). Since leaching tests start for practical reasons with material packed more or less dry into the column, a uniform initial concentration is not necessarily achieved during the first flooding of the column. Initial conditions as assumed in Equation (4) (uniform concentration distribution), would only be achieved if the material is first mixed with water, equilibrated and then packed into the column, which is not practical. During the first flooding of the column, especially less sorbing solutes are displaced from the inlet and higher concentrations occur towards the outlet, as illustrated in Figure 1 (see also Appendix E). This may be accounted for by subtracting the distance of the solute displaced initially ( / with > 1) in Equation (4): In this case the initial solute concentration ( , ) in the column effluent is in equilibrium with the initial concentration in the solids ( , = , / ) and higher than calculated in Equation (5), especially if values are low.  In this case the initial solute concentration (C w,eq ) in the column effluent is in equilibrium with the initial concentration in the solids (C w,eq = C s,ini /K d ) and higher than calculated in Equation (5), especially if K d values are low.
In order to explore the influence of mass transfer limitations on initial and longterm solute concentrations in column tests, two relevant mass transfer mechanisms are compared: (1) film diffusion where diffusion from the solid-water interface occurs through an aqueous boundary layer with a given thickness and (2) intraparticle pore diffusion where diffusion inside the porous particle limits mass transfer.

Desorption Kinetics Limited by Film Diffusion
The simplest model for the kinetic release of a solute from solids considers diffusion through an aqueous boundary layer surrounding spherical particles (Figure 2). Such film diffusion models are also widely used for the dissolution of minerals with high solubilities (e.g., salts). The solute release from the solid surface into the bulk water phase can be described by a linear driving force with constant mass transfer coefficient k = D aq /δ f : and d [L] denotes the thickness of the external film, the volume of water, the dry mass of the solids in the column and the particle diameter, respectively. A o (= 6 m d /(V w ρ s d)) is the specific surface area of the particles per unit volume of water in the column [m 2 m −3 = m −1 ] (the term 6/ρ s d) represents the specific surface area of spherical particles per dry mass, e.g., in m 2 g −1 ). C w is the concentration at the solid-water interface where local equilibrium conditions apply (C w = C s /K d ). The external film thickness (δ f ) can be estimated from empirical Sherwood numbers (Sh) and the particle diameter (d): Materials 2021, 14, x FOR PEER REVIEW 5 of 33 Figure 2. Scheme of mass transfer limited by film diffusion during the first flooding with fixed concentration at the interface (because the infiltrating water is always contacting fresh material as it advances).

Desorption Limited by Intraparticle Pore Diffusion
If the release of compounds from the solid phase is governed by intra-granular diffusion, e.g., within a porous grain (Figure 3), then mass transfer is described by Fick's second law in radial coordinates: with the boundary conditions , ( = , ) = (12) For an overview on empirical relationships for the estimation of Sherwood numbers see Appendix A. The mass balance in such two-phase systems expressed by their respective rates is: Thus, the solute mass gained (or lost) by the water phase equals the solute mass lost (or gained) from the solid phase. If V w and m d in a packed bed (porous media) are replaced by n and ρ b , the density of the solids (ρ s ) in Equation (7) drops out. Film diffusion coupled to the one-dimensional advection-dispersion equation (Equation (1)) yields: Using the finite volume method, the column is spatially discretized by a number of cells (see Figure A1) and the governing equation (Equation (10)) is solved iteratively by employing the Newton-Raphson scheme. Details of the numerical solution of the film diffusion model are presented in Appendix B.

Desorption Limited by Intraparticle Pore Diffusion
If the release of compounds from the solid phase is governed by intra-granular diffusion, e.g., within a porous grain (Figure 3), then mass transfer is described by Fick's second law in radial coordinates: with the boundary conditions is the radial coordinate in the sphere and D e [L 2 T −1 ] the effective diffusion coefficient.
is the concentration of solute in the intra-granular pore water. ε [-] denotes the intraparticle porosity. R [L] and ρ p [M L −3 ] (= ρ s (1 − ε)) denote the radius and bulk density of the particle (sphere).
Materials 2021, 14, x FOR PEER REVIEW 5 of 33 Figure 2. Scheme of mass transfer limited by film diffusion during the first flooding with fixed concentration at the interface (because the infiltrating water is always contacting fresh material as it advances).

Desorption Limited by Intraparticle Pore Diffusion
If the release of compounds from the solid phase is governed by intra-granular diffusion, e.g., within a porous grain (Figure 3), then mass transfer is described by Fick's second law in radial coordinates: with the boundary conditions [L] and [M L −3 ] (= (1 − )) denote the radius and bulk density of the particle (sphere).  For linear sorption with concentration independent distribution coefficients (i.e., C s = K d C w,intra ) Equation (11) becomes: where D a [L 2 T −1 ] is the apparent diffusion coefficient, defined as: Empirical studies showed that D e increases approximately with the square of the porosity [32]. This corresponds to a tortuosity τ [-] of the intra-granular pores-if expressed as a function of intra-granular porosity-of τ ≈ 1/ε.
Considering intraparticle diffusion, the advection-dispersion model (Equation (1)) can be rewritten as: The corresponding equilibrium concentration (C w,eq ) in water after first flooding can be rewritten as C w,eq = C s,ini /(ε/ρ p + K d ), which is slightly lower than for non-porous solids (C w,eq = C s,ini /K d ) because the intragranular pore space is assumed to be initially free of water. The deviation becomes insignificant with the increase of K d (ε/ρ p + K d ≈ K d when K d ε/ρ p ). Coupling the intraparticle pore diffusion model (Equation (11)) to the one-dimensional advection-dispersion equation (Equation (16)) allows for the expression of the change of the solute concentration in the bulk water: The intraparticle pore diffusion model (Equations (11)-(13)) was implemented numerically using a finite volume method where the spherical particle was discretized by a number of spherical shells of equal volume (based on the method of Jäger and Liedl [28]). The column was spatially discretized by a number of cells (see Figure A1) and all the governing equations (Equations (11)- (13) and (17)) were solved iteratively applying the Newton-Raphson scheme. Compared to the 1D film diffusion case, the intraparticle pore diffusion case is more complicated and becomes a 2D problem. Details of the numerical solution of the intraparticle pore diffusion model are given in Appendix C.

Set-Up of "Numerical" Column Tests
The boundary conditions of the numerical experiments are based on the set-up of column tests in daily practice in Germany [21,33,34]. Table 1 lists the relevant material properties and the parameter ranges applied. The saturation time for the first pore volume of the column and the contact time (after the first flooding period) were set to 5 h. Initially, experiments with uniform materials were simulated where the intraparticle porosity, distribution coefficient, aqueous diffusion coefficient and tortuosity were set the same for each grain size fraction. The Sherwood number in packed beds was estimated based on the empirical formula proposed by Liu et al. (2014) (Equation (A3)) [35]. In order to illustrate the influence of longitudinal dispersion on the initial concentration distribution in the column after the first flooding, we used fine particles (d p, f ine = 63 µm) where kinetics are very fast, and the local equilibrium assumption is valid. The numerical model did not consider dispersion beyond the outlet of the column. Non-linear sorption was simulated using the Freundlich model (C s = K f r C 1 n w where K f r and 1/n denote the Freundlich coefficient and Freundlich exponent, respectively).
Many factors may contribute to sample heterogeneity, such as grain size distribution and particle properties (sorption, porosities, etc.). To highlight the impact of particle size and properties we focused on two grain size classes and different fractions of sorptive/reactive particles in the sample. Distribution coefficients were varied in a range of 0.1-100 L kg −1 . Lower K d values (<0.1 L kg −1 ) were not considered here (this would have resulted in very high initial aqueous concentrations). If the K d values become large (K d > 100 L kg −1 ), then the differences between the pre-equilibrated case and the "first flooding" scenario vanish and effluent concentrations are constant over long time periods. The K d range chosen covers many frequent environmental contaminants, such as per-and polyfluoroalkyl substances (PFAS), chlorinated solvents, polycyclic aromatic hydrocarbons and some heavy metals.

Impact of Initial Conditions on Leaching
In order to investigate the impact of initial conditions on leaching behavior, we compared two scenarios: (1) a column filled with pre-equilibrated material where the initial concentration distribution in the column was uniform (C w,eq = C s,ini / K d + n ρ b ) and (2) columns with non-uniform concentration distributions after first flooding where concentrations increased towards the outlet (to a maximum of C w,eq = C s,ini /K d ) while at the inlet the solute was already depleted. To illustrate this, we used the film diffusion model with fine grain sizes, and thus, fast kinetics (local equilibrium conditions). Figure 4 shows the initial aqueous concentration distribution in the up-flow column test after the first flooding of the column compared to the pre-equilibrium case. The results show that the differences in the initial concentration profiles became smaller with increasing sorption. At high K d values, the deviation of the initial concentration profiles only occurred at the inlet of the column. At low K d values, very high concentrations are expected at the column outlet; in extreme cases this may lead to a density increase in the leachate and-especially if flow is stopped-to density driven flow within the column. This would cause dilution and lower leachate concentrations when the flow is restarted as was potentially observed by Naka et al. [36]. Density driven mixing may be caused by soluble materials, e.g., sulphate or other salts and not necessarily the target compounds. This phenomenon is quite similar to the case where the dispersion is taken into account (see bottom curves of Figure 4 and also Figures S1-S8 in SM), which leads to more "mixing" in the column and thus lower initial concentrations at the outlet, especially for low K d values.
concentration expected (e.g., = 0.1 L kg −1 ). The effect of initial conditions on normalized concentrations looks like a phase shift (see Figure 5, 1st row). This would lead to an underestimation of values derived from the pre-equilibrium analytical solution (Equation (4)) if the conditions in the column after the first flooding are not appropriately considered. The lower the , the earlier the cumulative leachate concentration reaches its maximum value ( , = 1000 µg kg −1 ). Dispersion shifts this point to later times (see  Figure 5 shows how the initial conditions (pre-equilibrated column and column after first flooding) influence the leaching curves. Compared to the pre-equilibrated case, a higher equilibrium concentration appeared at the outlet of the column after the first flooding period and more contaminant mass was released from the column at early times, followed by a rapidly decreasing concentration (see Figure 5, 2nd row). The deviations vanished with increasing K d and became almost insignificant for K d ≥ 10 L kg −1 . Dispersion also reduces differences between the pre-equilibrated and the first flooding case. At high K d values, the maximum concentrations were still achieved but the tailings became smoother. With the decrease of the K d value, the concentration gradients at the inlet became steeper and the "back" dispersion fluxes towards the outlet increased as well. In extreme cases, the peak concentration at the column outlet was smaller than the maximum concentration expected (e.g., K d = 0.1 L kg −1 ). The effect of initial conditions on normalized concentrations looks like a phase shift (see Figure 5, 1st row). This would lead to an underestimation of K d values derived from the pre-equilibrium analytical solution (Equation (4)) if the conditions in the column after the first flooding are not appropriately considered. The lower the K d , the earlier the cumulative leachate concentration reaches its maximum

Initial Conditions and Leaching with Mass Transfer Limitations
Mass-transfer limitations may change the picture considerably, with respect to initial conditions and the development of leachate concentrations over time. Figure 6 shows the influence of film diffusion (FD) compared to intraparticle pore diffusion (IPD) limited desorption on the initial concentration distribution in the column after the first flooding period. For large values, equilibration is achieved after shorter distances in the column because of the retardation of the clean water front. The deviations between FD and IPD are due to different mass transfer zone lengths, , . % (see Appendix D for a discussion

Initial Conditions and Leaching with Mass Transfer Limitations
Mass-transfer limitations may change the picture considerably, with respect to initial conditions and the development of leachate concentrations over time. Figure 6 shows the influence of film diffusion (FD) compared to intraparticle pore diffusion (IPD) limited desorption on the initial concentration distribution in the column after the first flooding period. For large K d values, equilibration is achieved after shorter distances in the column because of the retardation of the clean water front. The deviations between FD and IPD are due to different mass transfer zone lengths, X s,63.2% (see Appendix D for a discussion of the concept and calculation of this length for FD and IPD). For FD, the mass transfer zone length is independent of the K d value and proportional to the particle size (Equation (A28)). In contrast, for IPD the length of the mass transfer zone increases with particle size to the power of 3/2 (d 3/2 ) and decreases with √ K d (Equation (A32)) (e.g., X s,63.2% = 10 cm, 3.5 cm, and 1.1 cm for K d values of 0.1 L kg −1 , 1 L kg −1 and 10 L kg −1 (see Figure A2)). The length of the mass transfer zone for IPD is much longer than for FD, but differences decrease with increasing K d values. For K d values of 1 L kg −1 and 10 L kg −1 , the mass transfer zone lengths for IPD are much shorter than the column length (X col = 30 cm), which indicates that the equilibrium concentration is achieved at the outlet of the column after the first flooding. For small K d values (e.g., K d = 0.1 L kg −1 ), the equilibrium concentrations are not achieved at the outlet if dispersion is considered (see Figure 6, lower panel) although the mass transfer zone length (X s,63.2% = 10 cm) is still shorter than the column length. This is because the "clean" water front is close to the column outlet and dispersion "dilutes" the steep concentration gradients ("back dispersion"). The deviations between FD and fast kinetics almost vanish when dispersion is considered, indicating that with film diffusion, equilibrium is almost achieved. The development of the concentration distribution for IPD is also illustrated in animated graphs provided in the Supplementary Material (SM).   shows concentrations in aqueous leachates that correspond to the initial conditions shown in Figure 6. If leaching is limited by IPD, then the leaching process is slower and concentrations at early times are considerably lower than in the FD or the equilibrium model. This is due to the retarded diffusive transport within the tortuous intragranular pores and the correspondingly small effective diffusion coefficients (D e ). The contaminant release rate becomes lower and lower over time. Leachate concentrations decrease first with the square root of time (typical for transient diffusion) and then exponentially (see Figure 7 without dispersion, and Appendix D for details about the development of the internal mass transfer resistance over time). Note, the cumulative concentration curves confirm the mass conservation of the numerical solutions.

Nonlinear Sorption Isotherms
For many of the environmental contaminants and solid materials that are typically analyzed in column leaching tests, non-linear sorption isotherms describe the distribution of solutes between the aqueous and solid phases. The significance of this non-linearity for the development of the conditions in the column before the leaching starts has been analyzed exemplarily by defining Freundlich isotherms that result in the same "effective" for the aqueous concentration at equilibrium: = / , . Figure 8 shows the influence of nonlinear sorption on both the initial concentrations in the column and the leaching curves for the example of = 1 L kg −1 when no dispersion is considered. The differences in the concentration distribution before percolation starts are moderate. Concentration profiles tend to be smoother with nonlinear sorption with a slightly lower maximum concentration at the column outlet for low to mid values if dispersion is considered (see SM, Figure S1). Differences become more obvious in the tailing part of the leaching curves. Freundlich exponents smaller than 1 result in a longer tailing as is expected. The effect of nonlinear sorption looks similar to the dispersion effect, in both cases the leaching curves show more tailing (see SM, Figure S2). Nonlinearity of

Nonlinear Sorption Isotherms
For many of the environmental contaminants and solid materials that are typically analyzed in column leaching tests, non-linear sorption isotherms describe the distribution of solutes between the aqueous and solid phases. The significance of this non-linearity for the development of the conditions in the column before the leaching starts has been analyzed exemplarily by defining Freundlich isotherms that result in the same "effective" K d for the aqueous concentration at equilibrium: K f r = K d /C 1 n −1 w,eq . Figure 8 shows the influence of nonlinear sorption on both the initial concentrations in the column and the leaching curves for the example of K d = 1 L kg −1 when no dispersion is considered. The differences in the concentration distribution before percolation starts are moderate. Concentration profiles tend to be smoother with nonlinear sorption with a slightly lower maximum concentration at the column outlet for low to mid K d values if dispersion is considered (see SM, Figure S1). Differences become more obvious in the tailing part of the leaching curves. Freundlich exponents smaller than 1 result in a longer tailing as is expected. The effect of nonlinear sorption looks similar to the dispersion effect, in both cases the leaching curves show more tailing (see SM, Figure S2). Nonlinearity of sorption is notably less significant than kinetic limitations in the mass transfer mechanisms.

Impact of Heterogeneous Sample Composition
Real world materials that are typically investigated in leaching tests are not always homogeneous. Although the sample might be addressed as 'one material', its individual grains have different sizes and differ most likely also in other properties such as porosity, tortuosity, sorption capacity, etc., and may contain different amounts of the contaminants of interest. In order to illustrate the impact of material heterogeneity, we have carried out several numerical leaching experiments with hypothetical material compositions.
First, we consider three bi-modal material compositions. Each of these compositions consist of a fraction of contaminated particles (e.g., particulate organic carbon particles with high ) and another fraction of particles that neither contain contaminants nor possess any sorption capacity. If equilibrium conditions prevail during the first flooding and leaching, the heterogeneity of the sample does not matter, it is simply the average ( , ) that rules. The situation changes if mass transfer between the solid and the aqueous phases is limited due to diffusion (FD or IPD). If only a small fraction of the particles in the sample carries the compounds of interest, the volume of particles released by the compound and thus the surface area available for mass transfer becomes much smaller. This may lead to pronounced non-equilibrium conditions after first flooding (see, e.g., Equations (A26) and (A30)) and during leaching. Figure 9 shows a comparison of the initial concentration profiles in the column after the first flooding, as well as the corresponding leaching curves that would develop for the three bi-modal material compositions (100/10/1% of the material is contaminated at a = 1/10/100 L kg −1 , respectively). A small fraction of strong sorbents showed lower desorption rates compared to a large fraction of the weak sorbents. For this "exotic" case where only 1% of the particles carries all the contamination, initial nonequilibrium and long tailing was observed. This effect was

Impact of Heterogeneous Sample Composition
Real world materials that are typically investigated in leaching tests are not always homogeneous. Although the sample might be addressed as 'one material', its individual grains have different sizes and differ most likely also in other properties such as porosity, tortuosity, sorption capacity, etc., and may contain different amounts of the contaminants of interest. In order to illustrate the impact of material heterogeneity, we have carried out several numerical leaching experiments with hypothetical material compositions.
First, we consider three bi-modal material compositions. Each of these compositions consist of a fraction of contaminated particles (e.g., particulate organic carbon particles with high K d ) and another fraction of particles that neither contain contaminants nor possess any sorption capacity. If equilibrium conditions prevail during the first flooding and leaching, the heterogeneity of the sample does not matter, it is simply the average K d (K d,av ) that rules. The situation changes if mass transfer between the solid and the aqueous phases is limited due to diffusion (FD or IPD). If only a small fraction of the particles in the sample carries the compounds of interest, the volume of particles released by the compound and thus the surface area available for mass transfer becomes much smaller. This may lead to pronounced non-equilibrium conditions after first flooding (see, e.g., Equations (A26) and (A30)) and during leaching. Figure 9 shows a comparison of the initial concentration profiles in the column after the first flooding, as well as the corresponding leaching curves that would develop for the three bi-modal material compositions (100/10/1% of the material is contaminated at a K d = 1/10/100 L kg −1 , respectively). A small fraction of strong sorbents showed lower desorption rates compared to a large fraction of the weak sorbents. For this "exotic" case where only 1% of the particles carries all the contamination, initial nonequilibrium and long tailing was observed. This effect was very pronounced for intraparticle pore diffusion; the concentrations initially started on a plateau ("like equilibrium"), but then rapidly declined and showed a pronounced tailing and decrease with the square root of time (or LS). It may be noted, that longitudinal dispersion becomes less relevant if non-equilibrium conditions prevail at high K d values (see Figures S3 and S4 in SM). If such pronounced initial nonequilibrium is observed, then extended periods of time would be needed to equilibrate the water in the column with the solids (e.g., a manifold of the contact time of 5 h).  Samples consisting of mixtures of different particle sizes represent another typical and frequently occurring case of material heterogeneity. To illustrate the impact of such grain size heterogeneity, two bi-modal grain size distributions are considered here, introducing two grain sizes, coarse particles having a diameter of d p,coarse = 2000 µm, and fine particles with d p, f ine = 63 µm. The 1st hypothetical grain size distribution consists of 10% fine particles and 90% coarse particles, the 2nd distribution of 90% fine particles and 10% coarse particles.
If mass transfer is limited by film diffusion, the establishment of the initial conditions as well as leaching ( Figure 10) occurs under conditions close to equilibrium for both grain size distributions at all K d values (0.1, 1, and 10 L kg −1 ). While the shapes of all leaching curves are very similar, their locations are shifted in time according to the different K d values by a factor of 10. If intraparticle pore diffusion is considered, tailing is observed if coarse particles predominate. This applies to both, the development of initial conditions in the column and leaching. If fine particles predominate, the leaching is close to equilibrium at early times; at later times, tailing is observed with the typical square root of time behavior. Considering the dispersion effect, non-equilibrium concentrations can be seen at the column effluent after first flooding especially at low K d values (K d =0.1 L kg −1 ). Initial non-equilibrium conditions become more salient for intraparticle pore diffusion if coarse particles predominate (see Figures S5 and S6 in SM). coarse particles predominate. This applies to both, the development of initial conditions in the column and leaching. If fine particles predominate, the leaching is close to equilibrium at early times; at later times, tailing is observed with the typical square root of time behavior. Considering the dispersion effect, non-equilibrium concentrations can be seen at the column effluent after first flooding especially at low values ( = 0.1 L kg −1 ). Initial non-equilibrium conditions become more salient for intraparticle pore diffusion if coarse particles predominate (see Figures S5 and S6 in SM). Combining the heterogeneity of particle size ( ) and sorption capacity ( ), we consider three material compositions in the third case, which aims at showing circumstances where strong non-equilibrium conditions may be expected. For many materials this is probably not very realistic, but it may occur in material mixtures where a small particle Combining the heterogeneity of particle size (d) and sorption capacity (K d ), we consider three material compositions in the third case, which aims at showing circumstances where strong non-equilibrium conditions may be expected. For many materials this is probably not very realistic, but it may occur in material mixtures where a small particle fraction carries a "labile" contamination with a low K d vs. just a few large particles carrying the contaminant with a large K d . A hypothetical mixture containing 10% of fine particles with low sorption capacity (K d = 10 L kg −1 ) and 90% of coarse particles with high sorption capacity (K d = 100 L kg −1 ) is compared with two extreme cases where a hypothetical sample only contains pure fine particles with low sorption capacity and another hypothetical sample contains pure coarse particles with high sorption capacity. Figure 11 shows the initial concentration distribution for these three compositions after the first flooding period as well as the corresponding leaching curves. Sorption equilibrium is achieved rapidly if the sample consists of only fine particles with a small K d or only coarse particles with a high K d . Pure coarse material with a high K d shows a low equilibrium concentration (C w,eq = C s,ini /K d = 1000 µg kg −1 /100 L kg −1 = 10 µg L −1 ) while pure fine material with a low K d presents a much higher equilibrium concentration (C w,eq = C s,ini /K d = 1000 µg kg −1 /10 L kg −1 = 100 µg L −1 ) after a short flow distance. Interestingly, the mixed case where 10% of the column is fine material caused a high concentration which would be sorbed by the coarse materials leading to a slightly higher plateau concentration compared to pure coarse materials. The pollutants were redistributed between fine and coarse materials during the first flooding of the column. The concentration increase towards the outlet of the column in the mixed case is due to fast desorption from the fine material followed by slow sorption by the coarse material. The redistribution is almost complete at the inlet of the column because of the long residence time (t c = 5 h). Since the fine particles make up only to 10% of the total mass, they are already depleted in contaminant concentrations inside the column and in equilibrium with the coarse particles (reflecting both extreme cases). The front of the high concentration caused by the fine particles is already close to the outlet, while the rest is in equilibrium with the 90% coarse particle fraction.
The leaching curve of the mixed case (red lines) reflects the properties of the two pure (homogeneous) cases with either fine or coarse particles. Ten percent of the fine particles with low sorption capacity led to a peak effluent concentration which was only slightly lower than the equilibrium concentration of the pure fine particles with low sorption capacity. However, because the fine particles made up only 10% of the total mass, this peak concentration leached out rapidly and the eluate concentrations followed the coarse particles with a high sorption capacity for long time periods (blue curves). Although this may appear to indicate non-equilibrium conditions (because of the rapid initial decline of the concentrations followed by a plateau or "tailing"), leaching from fine and coarse particles occurs at, or close to equilibrium. Compared to FD, the IPD in the mixed sample showed a slower concentration decline because of the desorption kinetics of the IPD of the coarse particles was slower than the case of FD and on the long-term control release kinetics. For the cumulative mass release the mixed case is close to the coarse material for both the FD and the IPD, whereas the fine-grained particles showed a much higher and faster release (see Figure 11: bottom panel).

Summary and Conclusions
We conducted numerical simulations to investigate the release characteristics of low to strongly sorbing compounds ( = 0.1-100 L kg −1 ) in column leaching tests. Two different scenarios for the establishment of the initial conditions before the start of the leaching phase were considered: a fully pre-equilibrated column and a more realistic scenario where a column is flooded with water from the bottom. In order to highlight the effect of mass transfer limitations, two mechanisms are compared: film diffusion and intra-particle diffusion. Cases without and with dispersion illustrate how dispersive mixing may mask diffusion limited mass transfer. Furthermore, we looked into the impact of heterogeneous sample compositions in terms of reactive particle fractions and particle sizes. Since possible parameter combinations amount to almost infinite numbers, we have limited our analysis to just a few exemplary cases that illustrate the role of individual material properties. These few cases already show that virtually any leaching behavior can be produced with

Summary and Conclusions
We conducted numerical simulations to investigate the release characteristics of low to strongly sorbing compounds (K d = 0.1-100 L kg −1 ) in column leaching tests. Two different scenarios for the establishment of the initial conditions before the start of the leaching phase were considered: a fully pre-equilibrated column and a more realistic scenario where a column is flooded with water from the bottom. In order to highlight the effect of mass transfer limitations, two mechanisms are compared: film diffusion and intra-particle diffusion. Cases without and with dispersion illustrate how dispersive mixing may mask diffusion limited mass transfer. Furthermore, we looked into the impact of heterogeneous sample compositions in terms of reactive particle fractions and particle sizes. Since possible parameter combinations amount to almost infinite numbers, we have limited our analysis to just a few exemplary cases that illustrate the role of individual material properties. These few cases already show that virtually any leaching behavior can be produced with highly heterogeneous samples (depending on the mixing of different materials). The most important conclusions are: Initial conditions have a significant impact on leaching at low K d values (K d < 1 L kg −1 ). With increasing K d , the differences between pre-equilibrium and non-equilibrium conditions gradually vanish for K d > 10 L kg −1 (see Figure 5). Compounds with very low K d ("salts") would reach extremely high concentrations (K d << 1 L kg −1 ) at the column outlet (see Figure 4) potentially leading to enhanced dispersion due to density fingering. The K d values derived from retardation factors (R d in Equation (4)) would be underestimated if the conditions in the column after the first flooding are not appropriately considered, due to a "phase shift" in normalized concentrations curves (C w /C w,eq vs. LS).
Dispersion generally causes "smoothing" of concentrations gradients in the column and tends to "mask" film and intraparticle diffusion characteristics due to enhanced "mixing" of the solute within the column. It may lead to smaller initial concentrations at the column outlet after the first flooding period than expected for equilibrium; this is pronounced especially at low K d values (see Figure 7 and Figure S2), which may be interpreted as non-equilibrium, but is just a consequence of dilution by dispersive mixing.
Intraparticle pore diffusion (IPD) generally shows slower desorption kinetics than film diffusion (FD) through an aqueous boundary layer. This is due to the much smaller effective diffusion coefficient in the intraparticle pores and the large diffusion distance that develops inside the particle over time, resulting in the typical square root of time decrease of concentrations (a slope of 1/2 is observed in log-log plots of leaching curves, see Figures 7,9 and 10). IPD is more sensitive to the variation of particle sizes than FD (see Figure 10). Mass transfer limitations in an aqueous boundary layer commonly exists for surface adsorbed compounds and easily soluble solids ("salts"). Elements such as heavy metals, which are slowly released from the solid phase, would require much lower solid state diffusion coefficients; if reaction fronts propagate into the particle releasing metals, intraparticle (solid) diffusion models apply again (shrinking core), which are very similar to the IPD approach used here.
Non-linear sorption has little influence on the leaching test results if the "right" effective K d value is calculated for the proper concentration range (since for the nonlinear sorption the K d depends on the concentration, large deviations may occur if just the K f r is determined far away of the sample's concentration is used as "K d "); nevertheless, as concentrations decrease nonlinear sorption causes more tailing (see Figure 8).
Heterogeneous samples with only a small fraction of strongly sorbing particles lead to much slower desorption rates (because of less surface area), especially if mass release is limited by intraparticle pore diffusion (see Figure 9). In extreme cases (just 1% of the material is contaminated at K d = 100 × K d,av ), leaching may start at a plateau (suggesting equilibrium), but far below equilibrium concentrations (C w,peak C w,eq ) and concentrations later decrease further; The K d values derived from the initial aqueous concentration (K d = C s,ini /C w,peak ) would be overestimated while the K d values calculated from retardation factors would be underestimated.
In contrast to that, already relatively small amounts of fine particles lead to initial equilibrium, but long-term tailing occurs and is dominated by the coarse particle fraction, especially for intraparticle pore diffusion. Since our FD simulations are close to equilibrium, results are not very affected by grain size heterogeneity (see Figure 10). Material mixtures of small amounts of fine particles (10%) with low sorption capacity (K d =10 L kg −1 ) and large amounts of coarse particles with high sorption capacity (K d = 100 L kg −1 ), exhibit the respective characteristics of each of the individual components in different time periods (see Figure 11). Small amounts of fine particles with low sorption capacity dominate short term behavior of the mixtures and lead to a peak effluent concentration (C w,peak ) which approaches the equilibrium concentration expected for fine particles (see Figure 11). Since the mass fraction of fine particles is small (10%), the leachate concentrations drop rapidly and reach slightly higher equilibrium levels of 100% pure coarse particles due to the redistribution of pollutants between fine and coarse particles. Ten percent of fine particles with low sorption capacity causes a high equilibrium concentration which are sorbed by the coarse particles with high sorption capacity. K d values derived from the initial aqueous concentration (K d = C s,ini /C w,peak ) would be underestimated, while K d values derived from the following plateau concentration would be overestimated. Cumulative mass release, however, is often quite insensitive to mass transfer mechanisms (FD or IPD) especially for LS < 5 (see Figure 11).  Figure S4. Normalized concentrations (C w /C w,eq ) as well as cumulative concentrations (m cum ) in the column effluent vs. time (expressed as liquid to solid ratio: LS) for different combinations of sorbing particles and distribution coefficients (initial conditions depicted in Figure S3); left: without dispersion; right: with dispersion; solid lines: film diffusion cases, dashed lines: intraparticle diffusion cases, Figure S5. Initial concentration distribution in the column after the first flooding (up-flow) for two different bi-modal grain size distributions of fine and coarse particles; solid lines: fine particle mass fraction 10%; dashed lines: fine particle mass fraction 90%. (n = 0.45, v = 1.67 × 10 −5 m s −1 , α/x = 0 or 0.1, C s,ini = 1000 µg kg −1 , t c = 5 h, D aq = 1 × 10 −9 m 2 s −1 , ε = 0.05, d p,coarse = 2000 µm, d p, f ine = 63 µm); top panel: without dispersion; bottom panel: with dispersion, Figure S6. Influence of different grain size fractions and distribution coefficients on normalized concentrations (C w /C w,eq ) as well as cumulative concentrations (m cum ) in the column effluent vs. time (expressed as liquid to solid ratio LS); left: without dispersion; right: with dispersion; solid lines: fine particle mass fraction 10%; dashed lines: fine particle mass fraction 90%; kinetic parameters are the same as Figure S5, Figure S7. Initial concentration distribution in the column after the first flooding (up-flow) for different bi-modal material compositions of fine particles with low sorption capacity (K d = 10 L kg −1 ) and coarse particles with high sorption capacity; left: 100% coarse particles (K d = 100 L kg −1 ); middle: mixed sample with 10% fine particles; right: 100% fine particles; solid lines: film diffusion (FD), dashed lines: intraparticle diffusion cases (IPD); n = 0.45, v = 1.67 × 10 −5 m s −1 , α/x = 0 or 0.1, C s,ini = 1000 µg kg −1 , t c = 5 h, D aq = 1 × 10 −9 m 2 s −1 , ε = 0.05, d p,coarse = 2000 µm, d p, f ine = 63 µm; top panel: without dispersion; bottom panel: with dispersion, Figure S8. Leachate concentrations (C w ) as well as cumulative concentrations (m cum ) in the column effluent vs. time (expressed as liquid to solid ratio: LS) for different combinations of fine particles with low sorption capacity (K d = 10 L kg −1 ) and coarse particles with high sorption capacity (K d = 100 L kg −1 ); left: without dispersion; right: with dispersion; solid lines: film diffusion cases, dashed lines: intraparticle diffusion cases; kinetic parameters are the same as Figure S7.  Institutional Review Board Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available in this article.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Empirical Relationships for the Estimation of Sherwood Numbers
There are many studies available in the literature in which solid-liquid mass transfer in fluidized beds and flow through systems are investigated over a wide range of Reynolds numbers. Most of these correlations can be adequately described by the following equation: where Sh is the Sherwood number. A is a constant (theoretically = 2 for spherical particles in a stagnant infinite medium) and B is a constant to be determined by regression analysis of experimental data. Re and Sc denote the Reynolds and Schmidt number which are defined as: Equation (A3) is equivalent to Equation (A1) for θ = γ = 1/2. This Sherwood number correlation was applied in the numerical models. Sh numbers obtained for the chosen column setup were close to 2 indicating slow mass transfer close to the theoretical limit (2).

Appendix B. Film Diffusion Coupled to Advective-Dispersive Transport
Equation (10) in the main text shows the governing equations of film diffusion coupled to advective-dispersive transport. These partial differential equations are solved numerically using the finite volume method (as illustrated in Figure A1a).
Discretizing the transport operator in space while keeping the time derivative yields the following system of ordinary differential equations: where the subscripts + 0.5 and − 0.5 represent the corresponding parameter value between shells ( and + 1) and ( and − 1), respectively. Subscript j denotes the corresponding parameter value in volume j. Subscripts and + 1 denote the "old" and "new" time levels. Based on the boundary conditions (Equations (12) and (13), main text), the innermost shell and the outermost shell are treated specially. The solute concentration in the intragranular water phase at the new time step ( , ) can be expressed as: The approximation of the time derivative of Equation (A4) can be expressed as the concentration difference between the new and the previous time, divided by the time interval ∆t. A time weighting factor ϕ was used to navigate between implicit and explicit time integration. For ϕ = 0.5, the Crank-Nicholson-scheme is realized, whereas for ϕ = 0 and ϕ = 1, the fully implicit and explicit scheme is used, respectively.
where the indices k and k + 1 denote the corresponding concentration values at the previous time step and at the new time step.
In order to solve this system of equations, we may merge the two concentration vectors into a single one (C = [C w ; C s ]; with the semicolon being a line delimiter): with j [1, 2, . . . , N].
A standard method of solving non-linear ordinary equations is the Newton-Raphson scheme [45]. It is based on linearizing the residual function f (C k+1 ) at the current guess C k+1 guess of C k+1 . The residual function f (C k+1 ) is defined as: The residual function vector can be expressed as: The residual function vector becomes a zero vector if C k+1 is chosen right and a single step of the Newton-Raphson method can be denoted as: where J denotes the Jacobian matrix, which is the matrix of derivatives of all values of f C k+1 with respective to all values of C k+1 . The residual f C k+1 is reevaluated after updating C k+1 . If the resulting residual is not sufficiently close to zero, C k+1 guess is set to the last solution of C k+1 and Equation (A9) is reapplied. In our case, the Jacobian matrix can be derived analytically: In order to ensure the accuracy of the model, error control was employed at each time step and an error vector (∆C k+1 ) was used to monitor the difference between the old and new guess values, which is defined as: 2j, guess,old . . .
The iteration process stops if the maximum value of error vector ∆C k+1 max is smaller than the tolerable error e (e.g., e = 10 −15 ). for shell p = 2 to shell p = L − 1 and 1 ∆t εC k+1 1,j + ρ p K f r C k+1 After transformation : for shell 1 (or the innermost shell, p = 1) and: After transformation : Based on the mass balance, solute mass change in the external water phase (M w ) equals the solute mass change in the spherical particles; for better understanding, the simple case of particles with uniform size is shown: where F [M L −2 T −1 ] denotes the solute flux density into the external water phase. R and N p denote the radius and the total number of the spherical particles. The latter can be calculated by: The solute flux density into the external water phase is given by: Substituting F and N p with Equations (A18) and (A17) in Equation (A16) and taking advection and dispersion into account, the solute concentration in the external water phase at the new time step C k+1 w,j can be expressed by: After transformation : In order to solve this system of equations, we may merge the two concentration vectors to a single one (C = [C w ; C p ]; with the semicolon being a line delimiter): Using the Newton-Raphson scheme, the following residual function f (C k+1 ) is linearized at the current guess C k+1 guess of C k+1 [45]: The residual function vector can be organized as: Similar to the film diffusion case (see Appendix B), the C k+1 vector can be determined by Equation (A9) as well and the Jacobian matrix of intraparticle pore diffusion case can be expressed as: In order to ensure the accuracy of the model, an error vector (∆C k+1 ) was used as described above for Equation (A11): The iteration processes stop when the maximum value of the error vector ∆C k+1 max is smaller than the tolerable error e (e.g., e = 10 −15 ). which upon integration yields the following analytical solution for the initial condition C w(t=0) = 0 (desorption): The contact time in Equation (A27) can be substituted with the ratio of the travel distance (x) and the flow velocity (v). The length of the mass transfer zone is defined by setting the argument of the exponential function to −1, referring to the location x where the solute concentration in the groundwater reaches 63.2% of the equilibrium concentration.
Equation (A28) shows that the length of mass transfer zone depends on the flow velocity, inter-granular porosity as well as particle size, but is independent of the distribution coefficient.
If the length of the mass transfer zone is shorter than the column length (X col ), a concentration higher than 63.2% of the equilibrium concentration will be observed in the column effluent until the mass transfer zone arrives at the column outlet. The time needed to reach 63.2% equilibrium concentration at the column outlet equals: Considering fast kinetics ( X s → 0 ), t 63.2% (≈ X col /(v/R d )) is mainly dominated by the retarded seepage velocity (v/R d ).

Appendix D.2. Analytical Solution Based on the Intraparticle Pore Diffusion Model
Expressing internal mass transfer resistance by means of intraparticle pore diffusion, with mass transfer coefficient k = D e /δ p , where D e is the effective intraparticle diffusion coefficient (D e = D aq ε/τ ≈ D aq ε 2 ) and the mean square displacement δ p (δ p = √ πD a t c ) representing the diffusion distance, which grows with the square root of contact time between particles and water (t c ) at early times, leads to: The contact time between water and dry particles can be estimated by the ratio of particle size and flow velocity (t c = d/v).
For the initial condition C w(t=0) = 0, integration of Equation (A30) yields the following analytical solution: (A32) X s,63.2% based on intraparticle pore diffusion increases with particle size to the power of 3/2 (d 3/2 ) and decreases with the square root of the distribution coefficient ( √ K d ). The time required to observe the corresponding concentration of 63.2% of the equilibrium concentration at the column outlet can be determined with Equation (A29). In Figure A2, analytical solutions for the increase of the concentration in the first water parcel over distance (which here represents time: t = x/v) during the first flooding of the column are shown for FD (Equation (A28)) and IPD (Equation (A31)). The numerical solutions are described in Appendices B and C.
A comparison reveals that the analytical solutions and the numerical solutions overlap almost perfectly for both FD and IPD (see Figure A2). This verifies the accuracy of the numerical model. The length of the mass transfer zone for FD is 0.35 cm and independent of K d , and much shorter than for IPD with X s = 10 cm, 3.5 cm and 1.1 cm for K d values of 0.1 L kg −1 , 1 L kg −1 and 10 L kg −1 , respectively. The deviations between FD and IPD gradually vanish with increasing K d values. If the initial concentration in the column leachate is close to equilibrium, it may be used for the determination of K d (K d = C s,ini /C w,peak ); K d is overestimated if the initial effluent concentration does not reach equilibrium (C w,peak < C w,eq ). The length of the mass transfer zone (Equations (A28) and (A32)) may be used to assess equilibrium at the beginning of the column test. In Figure A2, analytical solutions for the increase of the concentration in the first water parcel over distance (which here represents time: = / ) during the first flooding of the column are shown for FD (Equation (A28)) and IPD (Equation (A31)). The numerical solutions are described in section A2 and A3.
A comparison reveals that the analytical solutions and the numerical solutions overlap almost perfectly for both FD and IPD (see Figure A2). This verifies the accuracy of the numerical model. The length of the mass transfer zone for FD is 0.35 cm and independent of , and much shorter than for IPD with = 10 cm, 3.5 cm and 1.1 cm for values of 0.1 L kg −1 , 1 L kg −1 and 10 L kg −1 , respectively. The deviations between FD and IPD gradually vanish with increasing values. If the initial concentration in the column leachate is close to equilibrium, it may be used for the determination of ( = , / , ); is overestimated if the initial effluent concentration does not reach equilibrium ( , < , ). The length of the mass transfer zone (Equations (A28) and (A32)) may be used to assess equilibrium at the beginning of the column test.

Appendix E. Comparison of Analytical and Numerical Solution (Code Verification)
In order to further confirm the accuracy of the numerical solution, the initial concentration distribution in the column after first flooding, as well as leaching curves are compared with the analytical solution (Equation (6)). The analytical solution is only valid for equilibrium sorption conditions and to compare it with the numerical solution, fine particles ( , = 63 μm) are used to get close to equilibrium (to fast FD kinetics). Figures

Appendix E. Comparison of Analytical and Numerical Solution (Code Verification)
In order to further confirm the accuracy of the numerical solution, the initial concentration distribution in the column after first flooding, as well as leaching curves are compared with the analytical solution (Equation (6)). The analytical solution is only valid for equilibrium sorption conditions and to compare it with the numerical solution, fine particles (d p, f ine = 63 µm) are used to get close to equilibrium (to fast FD kinetics). Figures A3 and A4 show the good agreement of both solutions. The slight deviations between analytical and numerical solutions, especially at low K d values, are due to kinetics in the numerical solution. Deviations gradually vanish with the increase of K d .