**Leaching of Titanium Dioxide Nanomaterials from Agricultural Soil Amended with Sewage Sludge Incineration Ash: Comparison of a Pilot Scale Simulation with Standard Laboratory Column Elution Experiments**

**Boris Meisterjahn 1 , Nicola Schröder 1 , Jürgen Oischinger 2 , Dieter Hennecke 1, \* , Karlheinz Weinfurtner 1 and Kerstin Hund-Rinke 1**


**Abstract:** Nanoscale titanium dioxide (nTiO<sup>2</sup> (Hombikat UV 100 WP)) was applied to sewage sludge that was incinerated in a large-scale waste treatment plant. The incineration ash produced was applied to soil as fertilizer at a realistic rate of 5% and investigated in pilot plant simulations regarding its leaching behavior for nTiO<sup>2</sup> . In parallel, the applied soil material was subject to standard column leaching (DIN 19528) in order to test the suitability of the standard to predict the leaching of nanoscale contaminants from treated soil material. Relative to the reference material (similar composition but without nTiO<sup>2</sup> application before incineration) the test material had a total TiO<sup>2</sup> concentration, increased by a factor of two or 3.8 g/kg, respectively. In contrast, the TiO<sup>2</sup> concentration in the respective leachates of the simulation experiment differed by a factor of around 25 (maximum 91.24 mg), indicating that the added nTiO<sup>2</sup> might be significantly mobilisable. Nanoparticle specific analysis of the leachates (spICP-MS) confirmed this finding. In the standard column elution experiment the released amount of TiO<sup>2</sup> in the percolates between test and reference material differed by a factor of 4 to 6. This was also confirmed for the nTiO<sup>2</sup> concentrations in the percolates. Results demonstrate that the standard column leaching, developed and validated for leaching prediction of dissolved contaminants, might be also capable to indicate increased mobility of nTiO<sup>2</sup> in soil materials. However, experiments with further soils are needed to verify those findings.

**Keywords:** nano titanium dioxide (nTiO<sup>2</sup> ); engineered nanomaterial (ENM); sewage sludge incineration (SSI); ENM containing sewage sludge ash (SSA); leaching; column elution; agricultural use

#### **1. Introduction**

Engineered nanomaterials (ENMs) applied, e.g., in consumer products, can be released to the environment during their use (e.g., release of silver nanoparticles (AgNPs) from facade painting [1] or TiO2-nanoparticles (TiO2-NPs) from sunscreens [2]), while, after use, a major fraction of the ENMs is supposed to be released to wastewater streams [3] and, furthermore, becomes attached to sewage sludge during wastewater treatment [4,5]. The majority of sewage sludge is incinerated and ends up in landfills [6].

Sewage sludge incineration ash (SSA) can be used for the production of phosphorous fertilizer (e.g., according German fertilizer ordinance [7]). This recycling route will be increased in the future because the recovery of phosphorus from sewage sludge and sewage sludge ash is required by law in Germany and other countries of the EU.

**Citation:** Meisterjahn, B.; Schröder, N.; Oischinger, J.; Hennecke, D.; Weinfurtner, K.; Hund-Rinke, K. Leaching of Titanium Dioxide Nanomaterials from Agricultural Soil Amended with Sewage Sludge Incineration Ash: Comparison of a Pilot Scale Simulation with Standard Laboratory Column Elution Experiments. *Materials* **2022**, *15*, 1853. https://doi.org/10.3390/ ma15051853

Academic Editor: Stefano Lettieri

Received: 16 December 2021 Accepted: 25 February 2022 Published: 1 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The use of materials in soils always requires a specific assessment to prevent soil contamination by re-use of waste materials. Integral part of the assessment is the determination of the source strength from those materials for leaching of contaminants. This is performed, e.g., by a standard column elution according to DIN 19528 for dissolved heavy metals and a set of organic substances.

So far, very little has been published regarding the release of nanoparticles from SSA. In a very recent study, Wielinski et al. [8] investigated the release of different nanomaterials from SSA in column experiments, but the column experiments described therein are very different to the standard column elution used for evaluation of source strength in the scope of regulation. Further, laboratory column experiments need references from the "real world" to assess the robustness of the laboratory results. However, we could not find any study regarding the release of nanomaterials from SSA after use in agricultural soils as fertilizer. Respective studies mostly focus on heavy metal contaminations (not in the form of NPs) or availability of phosphorous from SSA. Thus, there were two knowledge gaps identified for our study: (i) the question of leaching of NPs that might enter agricultural soils by fertilization with NP-contaminated SSA under realistic and environmentally relevant conditions and (ii) the predictive power of a column elution procedure actually used in the German soil protection ordinance [9] for NP leaching.

These gaps were addressed by the present study. In order to investigate the potential leaching of a representative ENM nano titanium dioxide (nTiO2) from SSA applied to soils, combined elution experiments in pilot and laboratory scale were performed. nTiO<sup>2</sup> was selected as test material, as it is among the most produced nanomaterials [10,11] and was used, e.g., in suncreens [2], as photocatalyst [12] or food additives [13]. Ash material was produced in a large-scale waste incinerator from incineration of sewage sludge amended with nTiO2. The resulting SSA was mixed with a reference soil at a rate of about 5%. This is about a factor of 10 greater than the maximum allowed application rate [14] and was considered to be a worst case scenario. The SSA amended soil was used for a pilot scale bioreactor trial in order to investigate nTiO<sup>2</sup> leaching from the material under controlled realistic conditions. In a timelapse experiment, three annual summer/winter cycles were simulated within 250 days to verify the influence of seasons on the leaching behavior. Those pilot scale simulations were accompanied by standard laboratory column elution experiments. Based on the results obtained in both tests, the suitability of the standard column elution with regard to a prediction of nTiO<sup>2</sup> leaching from soil treated with SSA should be determined. In case of a significant correlation, the standard laboratory column elution could be extended for risk assessment of nanomaterial leaching from SSA amended soils. So far, this has not been considered in soil protection, since transport of nanomaterials in soil cannot be described by common leaching models, as the sorption theory used for dissolved chemicals does not apply for nanomaterials [15].

#### **2. Materials and Methods**

#### *2.1. Preparation of ENM Containing Sewage Sludge Incineration Ash (SSA)*

For the production of the ENM-containing sewage sludge ash (SSA), the product Hombikat UV 100 WP from Co. VENATOR was used. The aqueous Hombikat UV 100 WP dispersion consists of 42% (*w*/*w*) nTiO<sup>2</sup> and about 6.5% (*w*/*w*) of poly acrylate acting as stabilizer. As a basis for the formulation of Hombikat UV 100 WP, nTiO2, in the form of the product Hombikat UV 100 with a primary particles size of <10 nm, was employed. Further information can be found in Börner et al. [16] and Oischinger et al. [17]. The application of Hombikat UV 100 WP lies mainly in the photocatalytic area, as the nTiO<sup>2</sup> is present in the anatase modification.

The sewage sludge incineration was performed at the Sewage Sludge Incineration Plant (SSIP) of the waste water treatment plant ZVK Steinhäule at Neu-Ulm. Annually a wastewater amount of 440,000 population equivalents is cleaned in the plant, resulting in 10,000 t of sewage sludge (dry matter (DM)) and 2500 t SSA [18]. The SSIP consists of a centrifuge, a dryer, a fluidised-bed incinerator with selective non-catalytic NO<sup>x</sup> reduction, an electrostatic precipitator, a 2-stage scrubber and an activated-carbon reactor with fabric filter. A schematic diagram of the plant is depicted in [16].

For the production of the reference ash material on the 1st day, the average amount of sewage sludge mounted up to 2171 kg h−<sup>1</sup> (DM) and, for the day with nTiO<sup>2</sup> injection, on the 2nd day, up to 2120 kg h−<sup>1</sup> (DM). During the measurement with nTiO<sup>2</sup> injection, 630 kg of nTiO<sup>2</sup> dispersion was added to the sewage sludge with a peristaltic pump over 6 h. On both days, the measured background concentration of titanium was <0.1 wt%. Hence, an augmentation of titanium in the sewage sludge of about 1.25 wt% was achieved by the injection of the dispersion. The SSA was sampled on the basis of the recommendations of LAGA PN98 [19] for moving waste, as far as they were applicable to the large scale plant [20].

#### *2.2. Simulation Experiments*

2.2.1. Leaching in Pilot Scale Simulation Reactors

For investigation of the release of nanoparticles from SSA applied in soil-related applications, experiments in pilot-scale simulation reactors were conducted. For this purpose, the incineration residue (with and without treatment with nTiO2) was mixed with an agricultural soil (refesol 04-A, details of the soil characteristics see Table 1).

**Table 1.** Refesol 04-A soil characteristics (loamy sand).


The sandy soil with higher organic carbon content was selected in order to simulate worst case conditions regarding potential nanoparticle release and mobility [21,22]. Per reactor, 700 kg of soil material (dry weight basis) were mixed with an amount of SSA, representing 5% with regard to the dry weight of the soil material, corresponding to approximately 190 t sludge ash per ha, taking into account an assumed plough share depth of 25 cm. At mixing, the soil was adjusted to a water content representing 50% of its water holding capacity (WHC), being optimal for microbial activity. The mixture was then subjected to an experiment simulating at least three annual seasonal cycles during a total incubation time of 250 days. Thus, the soil/ash mixture was incubated at 20–22 ◦C simulating a summer phase, followed by freezing of the mixture to −10 to −15 ◦C final temperature. The soil in the reactor was kept at this temperature for 14 days (simulated winter period). After each winter period and thawing of the soil/ash mixture, the soil was dug up and samples were taken for the column elution experiments (seasonal cycles and sampling dates see Figure 1). The soil was irrigated over a period of 12 days with 10 L portions until approximately 5 L of seepage/leachate water was collected. The leachate was analysed for total titanium concentration by ICP-OES after microwave assisted digestion and, in addition, was subjected to nanoparticle specific analysis by single particle (sp)ICP-MS. For the next cycle, the soil was again dried to a water content of 50% WHC by digging the soil up several times and leaving the reactor surface open while incubating at 20 ± 2 ◦C. After reaching the desired water content, the reactor was closed again and the next cycle started.

**Figure 1.** Temperature curve over the entire test period with three annual cycles, measured at 20, 40 and 60 cm distance from the edge of the reactor with treated soil ash/mixture. Dotted line marks the three sampling points after the winter phases (23 January, 24 April and 27 June).

#### 2.2.2. Column Elution Tests

At three time points, samples of the soil/ash mixture in the pilot-scale reactors were sampled and filled into glass columns (diameter 5.5 cm) according to the instructions given by the standard DIN 19528 [23] for soil column leaching. At both ends of the glass columns, quartz sand was used as a filtration layer. The column was then subjected to an elution according to DIN 19528 [23]. The column was saturated with deionized water from the bottom to the top layer. After saturation, the column was eluted with deionized water and eluate samples taken at five different solid/water ratios (0.3, 1, 2, 4 and 10 L/kg). The titanium concentration in the eluate was determined as described in Section 2.3.1 without any previous eluate filtration or centrifugation. The amount of soil per soil column and the percolation rates for saturation and percolation are presented in Table 2. The percolation rates were calculated according to DIN 19528:

$$q = \frac{l \times \pi \times r^2 \times n}{t \times 60} \tag{1}$$

with *q* = percolation rate (mL min−<sup>1</sup> ); *l* = length of soil column; *r* = inner radius of column; *n* = porosity (0.43); *t* = time (2 h for saturation, 5 h for percolation).

**Table 2.** Mass of soil per soil column and percolation rates for saturation and percolation.


#### *2.3. Chemical Analysis*

2.3.1. TiO<sup>2</sup> Concentrations in Reactor Leachates and Column Eluates

For determination of the total concentration of titanium in the leachates of the reactors, as well as for the column eluates, 10 mL of the respective sample was filled into a

TeflonTM digestion vessel (MLS, Leutkirch, Germany) and evaporated at 105 ◦C to dryness. The residue was then mixed with 4.8 mL of HNO<sup>3</sup> (69%, supra pure, Roth) and 0.2 mL hydrofluoric acid (HF, 40%, supra pure, Roth) and digested in a microwave (Ultraclave II, MLS, Leutkirch, Germany; digestion parameters: 220 ◦C, 30 min, 100 bar). After digestion, 1 mL boric acid was added in order to complex the remaining HF; the mixture filled up to a final volume of 15 mL with ultrapure water. The concentration of Ti in the solution was then determined by ICP-OES (wavelength 334.941 nm, Instrument: Agilent 5110, Agilent technologies). The obtained results were converted into TiO<sup>2</sup> concentrations under the assumption that no other titanium-containing phase was present in the samples.

#### 2.3.2. TiO<sup>2</sup> Concentrations in Soil/Ash-Mixed Samples

For determination of total titanium concentrations in mixtures of soil and sewage, sludge ash was applied for simulation in the pilot-scale reactors; in addition, column elution experiments, 5 aliquots of 5 g (fresh weight) each, were taken and mixed to obtain a representative sample. The samples were dried at 105 ◦C until constant weight. Aliquots of approximately 200 mg of the dried residues were then weighted in TeflonTM digestion vessels and 1 mL HF (40% supra pure, Roth) and 4 mL HNO<sup>3</sup> (69%, supra pure, Roth) added. The samples were subjected to microwave-assisted digestion (Ultraclave II, MLS, digestion parameters: 220 ◦C, 30 min, 100 bar). After digestion, 5 mL of boric acid was added to complex remaining HF and the digestate filled up to a final volume of 15 mL with ultrapure water. The solution was analysed for titanium by ICP-OES at a wavelength of 334.941 nm.

#### 2.3.3. Single Particle (sp)ICP-MS Analysis of Reactor Leachates and Column Eluates

In addition to the determination of the total titanium concentration, number-based particle size distributions were determined by single particle inductively coupled plasma mass spectrometry (spICP-MS) [24,25]. The aqueous samples were directly measured without any further sample preparation, except for dilution with ultrapure water. The analyses were performed using a triple-quadrupole ICP-MS instrument (ICP-QQQ-MS, Agilent 8900, Agilent Technologies, Waldbronn, Germany). The dwell time in the single particle measurement mode of the ICP-MS was set to 100 µs and time-resolved signals were recorded on the selected *m*/*z* for 60 s. Peak detection and integration was conducted automatically by the Agilent MassHuntersoftware. Conversion of signal heights of particle spikes into particle sizes were performed by application of a calibration with a dissolved Ti standard. To apply a dissolved calibration for size calculation of nanoparticles from their signal spikes, the nebulization or transport efficiency in the interface was determined by analyzing a gold nanoparticle standard of known concentration and size [24]. Dispersions of 60 nm gold nanoparticles (AuNPs 60 nm, BBI solutions, Kent, UK) were used for the determination of the nebulization efficiency and were prepared freshly on the day of measurement. The samples were diluted in ultrapure water by a factor of 102–10<sup>5</sup> for measurement in order to reach a particle concentration of 200–2000 particle events per minute.

Due to possible interferences for the most abundant titanium isotope <sup>48</sup>Ti caused by the calcium isotope <sup>48</sup>Ca present in the leachate and eluate matrices, titanium was measured in the MS/MS mode with ammonia as reaction gas (10% NH<sup>3</sup> + 1 mL/min He). Thus, titanium was measured as [48TiNH]<sup>+</sup> with a *m*/*z* of 63. The threshold between background and particle signals was defined based on visual inspection of the measured signal distributions. The conversion of signal distributions into number-based size distributions, as well as particle number concentrations, was performed by the Agilent MassHunter software.

#### **3. Results**

#### *3.1. Characterisation of the Produced SSA*

The two SSA produced differed in their titanium content as expected. For the SSA, from the reference day, an average titanium concentration of about 0.3 wt% (dry mass

(DM)) with a standard deviation (SD) of 0.04 was determined, whereas the SSA treated with nTiO<sup>2</sup> had a mean value of 2.9 wt% (DM) with an SD of 0.19 [26]. Hence, the titanium content in the nTiO2-SSA was increased compared to the reference measurement by a factor of about 9.7 [26]. Further information on characterisation of the used SSA can be found elsewhere [26].

#### *3.2. Leaching from Pilot Scale Simulation Reactors*

#### 3.2.1. Determination of Total Titanium Content in Mixtures of Soil/SSA and in Leachates

Figure 2 shows the determined total titanium concentrations in both mixtures of soil with SSA with additional nTiO<sup>2</sup> or without nTiO<sup>2</sup> treatment (reference samples) which were filled into the pilot scale reactors. The SSA/soil mixture containing nTiO2-treated ash shows a titanium concentration elevated by a factor of approximately two compared to the reference sample. Thus, at least this ratio should also be found in the leachates.

**Figure 2.** Total TiO<sup>2</sup> concentrations (converted Ti-concentrations) in sludge ash/soil mixtures. left column: Concentration in samples using sludge ash without treatment with nTiO<sup>2</sup> , right column: Concentration in samples using SSA treated with nTiO<sup>2</sup> . Error bars refer to standard deviation of four subsamples that were digested and analysed.

However, as shown in Figure 3, the titanium/converted TiO<sup>2</sup> content in the leachates of treatment and reference reactors collected at the end of the three simulated annual seasonal cycles was found to be higher for the treated samples compared to the references by a factor of ~25. Titanium-containing particles are released from the treatment reactor in a significantly higher amount than expected from the total titanium content of the respective mixtures, thus, indicating that the added nTiO<sup>2</sup> (or titanium containing particles) might be more mobilisable than the titanium-bearing particles (most likely also being TiO2) already present in the soil and sludge ash.

**Figure 3.** Total TiO<sup>2</sup> content (derived from Ti-concentrations) determined in leachates from pilot scale simulation reactors collected after watering of the soil at end of the respective preceding winter phase.

#### 3.2.2. SpICP-MS/MS Analyses of Leachates

The significantly higher release of titanium in the leachates from the treatment reactor is also visible in the nanoparticle specific analysis of the leachates by spICP-MS, as shown in Figure 4. The higher TiO<sup>2</sup> particle discharge in leachates after the first winter period (23 January) is 3.7 times higher in the treated reactor and rises to 640 times after the second simulated winter phase (17 April) and drops to 8.2 times after the third winter phase (27 June) (Figure 4). Again, the differences in discharge are much higher than could be expected from the total Ti-content determined for the soil/ash mixtures applied to the two reactors.

**Figure 4.** TiO<sup>2</sup> particle concentrations (Particles/L) determined in leachates from pilot scale simulation reactors collected after watering of the soil at end of the respective preceding winter phase.

#### *3.3. Column Elution Experiments*

The percolate samples obtained from standard column eluate experiments were analysed analogously to the leachates from the pilot scale experiments for their total TiO<sup>2</sup> content (ICP-OES analysis after HF digestion) and particle size distribution (spICP-MS).

Figure 5 shows the TiO<sup>2</sup> content in the percolates of soil/ash-material taken at end of the three simulated annual seasonal cycles from the treatment and reference reactor to investigate the behaviour during soil-related use. In the percolates, an increased TiO<sup>2</sup> load was detected for the treated SSA/soil mixture compared to the reference mixture. In contrast to the analyses of the reactor leachates, the ratio between TiO<sup>2</sup> released from nTiO2-treated material and the reference material was found to be more or less stable during all three sampling times. Already, percolates from the first sampling (17 January) showed increased TiO<sup>2</sup> contents in the nTiO<sup>2</sup> treated SSA/soil mixture. This was most likely due to the very high water:solid ratios used in the column elution compared to the pilot scale experiment. However, the maximum TiO<sup>2</sup> concentrations determined are in the same range as the maximum concentration found in the pilot scale experiment.

**Figure 5.** Total TiO<sup>2</sup> content (µg) in column percolates from reference and nTiO<sup>2</sup> treated SSA/soil mixtures out of pilot scale simulation reactors collected at three different sampling times.

The spICP-MS/MS analyses of the column percolates (reference and nTiO<sup>2</sup> treatment) show a higher nTiO<sup>2</sup> particle discharge from the nTiO<sup>2</sup> treated SSA/soil mixture (Figure 6). The particle number in the percolates of the reference and TiO<sup>2</sup> treated SSA/soil mixtures increases in both cases with the water/solid ratio (W/S). In the TiO2-treated percolates, the particle concentration at a W/S of 10 L/kg at the first sampling (17 January) extend concentration in the reference percolates by a factor 10. During the following samplings (17 April and 21 June), particle counts in the nTiO<sup>2</sup> percolates extend the reference percolates by a factor of 4. In general, the nTiO<sup>2</sup> analysis confirms the findings of the total TiO<sup>2</sup> analysis, indicating that most of the TiO<sup>2</sup> determined was due to leached nTiO2.

**Figure 6.** nTiO<sup>2</sup> particle concentrations (Particles/L) determined in column percolates from reference and nTiO<sup>2</sup> treated SSA/soil mixtures out of pilot scale simulation reactors collected at three different sampling times.

#### **4. Discussion**

In this study, nTiO<sup>2</sup> was used as an exemplary ENM for the preparation of ENMcontaining SSA produced for the soil simulation experiments. In the simulation of agricultural use of this nTiO<sup>2</sup> amended SSA, a clearly increased release of nTiO<sup>2</sup> was found in the leachates of the pilot scale simulation reactors after a certain time compared to a reference ash from non-amended sewage sludge. The amount of release was significantly higher than could be expected due to the differences in the total TiO<sup>2</sup> contents in the soil-ash mixture relative to the reference ash. A cause for the comparatively high releases could not be determined within the project. After three freezing and thawing cycles, the amount of TiO<sup>2</sup> released from the reference was increased significantly, indicating mobilization by changes of the soil structure due to frost wedging. This is also observed in increasing amounts of TiO<sup>2</sup> in the column experiments throughout the three samplings. A difference between the leaching from the simulation reactors and standard column elution is the significantly higher ratio between TiO<sup>2</sup> released from reference and treated SSA in the reactors compared to the ratios observed for the column experiments, which were nearer to the ratio of total Ti contents of both materials (Figure 2) at the end of the study. For instance, while, for the reactors, the release of TiO2, expressed in terms of total concentrations, from treated samples compared to the reference, was elevated by a factor of approx. 25, the analogous ratio for the standard column experiments was found to be approx. 3 (Figures 3 and 5). The differences are not that pronounced if the particle counts from spICP-MS analyses are assessed with factors of approx. 10 (leachates of treated-against-reference reactors) and 4 (column percolates of treated-against-reference soil/ash-mixtures). Thus, the particle counts agree with the total titanium on a relative scale for the column experiments, whereas, there seems to be an underestimation for the reactor leachates. This might be explained by losses of, e.g., bigger particles, which would significantly contribute to the total titanium mass in the sample, due to sedimentation prior to measurement or the underestimation

of the number of bigger particles due to the applied analysis time of 60 s. In contrast, the total Ti determination after digestion captures the complete amount of Ti in the sample. However, in principle, the spICP-MS analyses give similar results to the determination of the total concentrations and can possibly be a fast alternative analytical method that, additionally, provides also information on particle sizes.

Looking further and in more detail to the obtained data, some differences can be noticed not only between the elution behavior of the test materials but also between the elution procedures itself. The elution procedure in the column and the simulation experiment were very different. In the column elution, within 7 days, a water to solid ratio of 10 L per kg soil material was applied. In the simulation, the maximum water to solid ratio used was 0.5 L per kg of soil material over a time of 250 days. This results in very different flow velocities in both approaches; it is not surprising that these huge differences in the elution conditions give different results in eluate analysis.

This significant influence of different elution conditions is strikingly visible if the released TiO<sup>2</sup> amounts are related to the total amount of soil/ash mixture in the respective experiments. For instance, for the last sampling date, after three seasonal cycles, around 130 µg/kg (released amount (Figure 3) referred to 700 kg of soil/ash (dry weight) material per reactor) was released from the reactor, while, for the column elution experiments, 5172 µg/kg (TiO<sup>2</sup> content in percolate referred to 1.1 kg soil/ash per column) was released. Thus, the column leaching experiments are overestimating the leaching for this kind of particle in comparison to the more realistic simulation experiments.

The higher relative release of TiO<sup>2</sup> from reference samples in the column experiments indicates that the natural "background" TiO<sup>2</sup> might be attached differently to the soil matrix, as it is more easily mobilized under the experimental conditions compared to the simulation experiments.

Therefore, the main finding of the study is that the standard column experiment according to DIN 19528 might be suitable to show increased mobility and emissions of nTiO<sup>2</sup> during soil-related recycling. This finding is confirmed by simulation experiments under more realistic conditions. It proves the general suitability of the standardized column elution method for predicting an increased release of nTiO<sup>2</sup> from a soil-ash mixture.

However, due to the limitations of a pilot scale experiment, which does not allow to test a reasonable number of different soils that would be needed for a more general statement, results cannot be generalized at this point. However, they indicate that increased leaching of ENM measured in a standard soil column approach relative to a reference material might serve as reasonable prediction of increased leaching of ENM under realistic environmental conditions. As a next step, a larger data set should be determined on both SSAs using a set of different soils for standard column elution to obtain the influence of different soil matrices.

#### **5. Conclusions**

The results of this study demonstrate that a standard soil column experiment might be suitable to reliably indicate an increased mobility of nTiO<sup>2</sup> from SSA mixed with agricultural soil relative to a respective reference material containing only "background" TiO2. However, due to the limited data set of so far just one soil general conclusions cannot be made at this point. Still, the result is important for future exposure assessment, as it might offer an experimental laboratory tool to determine the risk of leaching from a recycling material.

The use of an accepted standard (DIN 19528) for the determination of the mobility offers the chance to get a new standard for testing ENM in soil material just by extension of the scope of an existing standard. In order to establish effective thresholds and to make them available for regulatory practice, respective standards for laboratory use are required. However, further research is necessary to verify current results and to base findings on a broader database before standardization processes might start. Laboratory elution experiments could give a useful contribution to the classification and safer use of incineration ashes, which contain not only potential pollutants but also many components that are too valuable for landfill disposal. Even if no direct transfer from the laboratory elution tests on the simulation experiments can be derived, the general concept is confirmed.

**Author Contributions:** Conceptualization, B.M., D.H. and K.H.-R.; methodology, B.M., D.H., K.H.-R.; validation, B.M., N.S. and D.H.; formal analysis, N.S. and B.M.; investigation, N.S. and K.W; resources, K.W., J.O.; data curation, B.M. and N.S.; writing—original draft preparation, B.M. and N.S.; writing review and editing, B.M., D.H. and K.W.; visualization, N.S.; supervision, D.H. and B.M.; project administration, D.H. and B.M.; funding acquisition, D.H., J.O., K.H.-R. and B.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the German Federal Environment Agency (UBA) and was part of the project "Investigations on the possible release of nanoparticles in the deposition and soil-related application of mineral wastes" (FKZ 3712 33 327).

**Data Availability Statement:** Restrictions apply to the availability of these data. Data was obtained from the project funded by the German Environment Agency (FKZ 3712 33 327) and are available from the corresponding authors with the permission of the German Environment Agency.

**Acknowledgments:** The authors acknowledge the plant operator at ZVK Steinhäule for his support during the experiments.

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

#### **References**


**Binlong Liu, Michael Finkel and Peter Grathwohl \***

Center for Applied Geoscience, University of Tübingen, Schnarrenbergstraße 94-96, 72076 Tübingen, Germany; binlong.liu@uni-tuebingen.de (B.L.); michael.finkel@uni-tuebingen.de (M.F.)

**\*** Correspondence: grathwohl@uni-tuebingen.de

**Abstract:** 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 (*K<sup>d</sup>* < 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 *K<sup>d</sup>* 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 *K<sup>d</sup>* ) 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 *K<sup>d</sup>* ), 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.

**Keywords:** leaching test; equilibrium condition; non-equilibrium condition; modelling; sorption kinetics; non-linear sorption; heterogeneity

#### **1. Introduction**

Leaching tests are widely used for the determination of contaminant release rates from soils [1–4], recycling materials [5–11], construction products [12–16], radioactive and other waste materials [17,18]. Compared to traditional batch shaking tests, column tests are preferred for assessing the risk of release of potential pollutants into groundwater or surface waters because they are closer to natural conditions [19,20]. Different mechanisms controlling desorption kinetics may result in complex leaching behaviors. While initially the observed column effluent concentration often reflects equilibrium conditions between the solid phase (incl. intraparticle pores) and the mobile aqueous phase [21,22], the concentrations decrease and often an extended tailing is observed due to slow desorption processes such as intraparticle diffusion [23–25].

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

**Citation:** Liu, B.; Finkel, M.; Grathwohl, P. Mass Transfer Principles in Column Percolation Tests: Initial Conditions and Tailing in Heterogeneous Materials. *Materials* **2021**, *14*, 4708. https://doi.org/ 10.3390/ma14164708

Academic Editor: Saeed Chehreh Chelgani

Received: 18 July 2021 Accepted: 16 August 2021 Published: 20 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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–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.

#### **2. Theory and Background**

#### *2.1. 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:

$$\frac{\partial}{\partial t}(n\,\mathcal{C}\_{w} + \rho\_{b}\,\mathcal{C}\_{s}) + \frac{\partial}{\partial \mathbf{x}}\left(n\,\boldsymbol{\upsilon}\,\mathcal{C}\_{w} - n\,D\_{L}\frac{\partial \mathcal{C}\_{w}}{\partial \mathbf{x}}\right) = \mathbf{0} \tag{1}$$

where *v* [L T−<sup>1</sup> ], *n* [-] and *D<sup>L</sup>* (= *αv* + *Dp*) [L<sup>2</sup> T −1 ] denote the seepage velocity of the water, the intergranular porosity and the longitudinal dispersion coefficient. *α* [L], *D<sup>p</sup>* (=*nDaq*) [L<sup>2</sup> T −1 ] and *Daq* [L<sup>2</sup> T −1 ] denote the dispersivity, the pore diffusion coefficient and the aqueous diffusion coefficient of the solute. *x* [L] and *t* [T] are the length of the column and time. *ρ<sup>b</sup>* (<sup>=</sup> (<sup>1</sup> <sup>−</sup> *<sup>n</sup>*)*ρs*) [M L−<sup>3</sup> ] is the dry bulk density of the packed bed in the column (porous media; *ρ<sup>s</sup>* is the solids density). For local equilibrium conditions the concentration in the solid phase (*Cs*) is in equilibrium with the solute concentration in water (*Cw*) and the distribution coefficient *K<sup>d</sup>* (= *Cs*/*Cw*) allowing for the calculation of the respective concentrations. Under these conditions, Equation (1) can be simplified as:

$$\frac{\partial \mathbb{C}\_w}{\partial t} = \frac{D\_L}{R\_d} \frac{\partial^2 \mathbb{C}\_w}{\partial x^2} - \frac{v}{R\_d} \frac{\partial \mathbb{C}\_w}{\partial x} \tag{2}$$

where *R<sup>d</sup>* [-] represents the retardation factor, defined as:

$$R\_d = 1 + K\_d \frac{\rho\_b}{n} \tag{3}$$

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]: ௪, = 1 − 0.5 ⎣ ⎢ ⎢ erfc ⎝ ⎛ ௗ 2ට ௗ ⎠ <sup>⎞</sup> + exp ൬ ൰ erfc ⎝ ⎛ ௗ 2ට ௗ ⎠ ⎞ ⎦ ⎥ ⎥

<sup>௪</sup>

⎡

 − 

$$\frac{\mathcal{C}\_{w}}{\mathcal{C}\_{w,aq}} = 1 - 0.5 \left[ \text{erfc} \left( \frac{\mathbf{x} - \frac{\mathcal{v}}{\mathcal{R}\_d} t}{2 \sqrt{\frac{D\_L}{\mathcal{R}\_d} t}} \right) + \exp \left( \frac{\mathbf{v} \cdot \mathbf{x}}{D\_L} \right) \text{erfc} \left( \frac{\mathbf{x} + \frac{\mathcal{v}}{\mathcal{R}\_d} t}{2 \sqrt{\frac{D\_L}{\mathcal{R}\_d} t}} \right) \right] \tag{4}$$

where erfc denotes the complementary error function. The aqueous concentration at equilibrium, *Cw*,*eq*, can be calculated from the initial solid concentration (*Cs*,*ini*) accounting for the mass balance when the contaminant mass in the water is equilibrated with the mass in the solid: ௪, = ௦, ௗ + / − −

$$\mathcal{C}\_{w,\epsilon q} = \frac{\mathcal{C}\_{s,ini}}{\mathcal{K}\_d + \frac{\mathcal{U}}{\rho\_b}} \tag{5}$$

+

⎤

The ratio *n*/*ρ<sup>b</sup>* [L<sup>3</sup> M−<sup>1</sup> ] 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−<sup>1</sup> for a column test with a porosity of *n* = 0.4 and a solid density of *ρ<sup>s</sup>* = 2.65 g cm−<sup>3</sup> , compared to e.g., 10 L kg−<sup>1</sup> 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<sup>d</sup>* with *R<sup>d</sup>* > 1) in Equation (4): /ௗ ௗ <sup>௪</sup> ௪, = 1 − 0.5 ⎣ ⎢ ⎢ ⎡ erfc ⎝ ⎛ ቀ − ௗ ቁ − ௗ 2ට ௗ ⎠ ⎞

$$\begin{array}{c} \frac{\mathbf{C}\_{w}}{\mathbf{C}\_{w,q}} = 1 - 0.5 \quad \left[ \text{erfc} \left( \frac{\left( x - \frac{x}{R\_d} \right) - \frac{x}{R\_d} t}{2 \sqrt{\frac{D\_d}{R\_d} t}} t \right) \begin{array}{c} \\ \int \\ \hline \\ \frac{\left( x - \frac{x}{R\_d} \right)}{D\_d} \end{array} \right] \\ + \exp \left( \frac{\mathbf{c} \left( x - \frac{x}{R\_d} \right)}{D\_d} \right) \text{erfc} \left( \frac{\left( x - \frac{x}{R\_d} \right) + \frac{x}{R\_d} t}{2 \sqrt{\frac{D\_d}{R\_d} t}} \right) \end{array} \tag{6}$$

 ௗ = 2 **Figure 1.** Initial concentration distribution in a column of length *x* for the "pre-equilibrated" case (dashed line) and after the first flooding of an initially dry column from the bottom (solid line); no dispersion, *R<sup>d</sup>* = 2, after Grathwohl and Susset, 2009 [21].

In this case the initial solute concentration (*Cw*,*eq*) in the column effluent is in equilibrium with the initial concentration in the solids (*Cw*,*eq* = *Cs*,*ini*/*K<sup>d</sup>* ) and higher than calculated in Equation (5), especially if *K<sup>d</sup>* 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.

#### *2.2. 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* = *Daq*/*δ<sup>f</sup>* :

$$\frac{\partial \mathbb{C}\_w}{\partial t} = k \, A^o \left( \mathbb{C}'\_w - \mathbb{C}\_w \right) = \frac{D\_{aq}}{\delta\_f} \frac{m\_d \, \mathsf{6}}{\rho\_s \, d \, V\_w} \left( \mathbb{C}'\_w - \mathbb{C}\_w \right) \tag{7}$$

where *δ<sup>f</sup>* [L], *V<sup>w</sup>* [L<sup>3</sup> ], *m<sup>d</sup>* [M] 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 *md*/(*V<sup>w</sup> ρ<sup>s</sup> d*)) is the specific surface area of the particles per unit volume of water in the column [m<sup>2</sup> m−<sup>3</sup> = m−<sup>1</sup> ] (the term 6/*ρ<sup>s</sup> d*) represents the specific surface area of spherical particles per dry mass, e.g., in m<sup>2</sup> g −1 ). *C* ′ *<sup>w</sup>* is the concentration at the solid-water interface where local equilibrium conditions apply (*C* ′ *<sup>w</sup>* = *Cs*/*K<sup>d</sup>* ). The external film thickness (*δ<sup>f</sup>* ) can be estimated from empirical Sherwood numbers (*Sh*) and the particle diameter (*d*):

$$\text{Sh} = \frac{\text{kd}}{\text{D}\_{aq}} = \frac{d}{\delta\_f} \quad \rightarrow \quad \delta\_f = \frac{d}{\text{Sh}} \tag{8}$$

**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).

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:

ቈ

௪,௧( = , ) = <sup>௪</sup>

( = 0, ) = 0

௪,௧ <sup>−</sup>

<sup>ଶ</sup> <sup>+</sup>

−

<sup>−</sup> = <sup>௦</sup>

$$V\_w \frac{\partial \mathcal{C}\_w}{\partial t} = -m\_d \frac{\partial \mathcal{C}\_s}{\partial t} \tag{9}$$

(1−))

௪,௧ 

Thus, the solute mass gained (or lost) by the water phase equals the solute mass lost (or gained) from the solid phase. If *V<sup>w</sup>* and *m<sup>d</sup>* in a packed bed (porous media) are replaced by *n* and *ρ<sup>b</sup>* , the density of the solids (*ρs*) in Equation (7) drops out. Film diffusion coupled to the one-dimensional advection-dispersion equation (Equation (1)) yields:

$$\frac{\partial \mathcal{C}\_w}{\partial t} = D\_L \frac{\partial^2 \mathcal{C}\_w}{\partial \mathbf{x}^2} - v \frac{\partial \mathcal{C}\_w}{\partial \mathbf{x}} + \frac{D\_{aq}}{\delta\_f} \frac{6}{n \, d} \frac{(1 - n)}{n \, d} \left(\frac{\mathcal{C}\_s}{K\_d} - \mathcal{C}\_w\right) \tag{10}$$

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.

#### *2.3. 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:

ቈ

$$\frac{\partial}{\partial t} \left( \varepsilon \, \mathbf{C}\_{w,intra} + \rho\_p \mathbf{C}\_s \right) = D\_\varepsilon \left[ \frac{\partial^2 \mathbf{C}\_{w,intra}}{\partial r^2} + \frac{2}{r} \frac{\partial \mathbf{C}\_{w,intra}}{\partial r} \right] \tag{11}$$

<sup>ଶ</sup> <sup>+</sup>

with the boundary conditions ൫ ௪,௧ + ௦൯=

$$\mathcal{C}\_{w,intra}(r = R, t) = \mathcal{C}\_{w} \tag{12}$$

$$\frac{\partial \mathbb{C}\_{w,intra}}{\partial r}(r=0,t) = 0 \tag{13}$$

*r* [L] is the radial coordinate in the sphere and *D<sup>e</sup>* [L<sup>2</sup> T −1 ] the effective diffusion coefficient. *Cw*,*intra* [M L−<sup>3</sup> ] is the concentration of solute in the intra-granular pore water. *ε* [-] denotes the intraparticle porosity. *R* [L] and *ρ<sup>p</sup>* [M L−<sup>3</sup> ] (= *ρs*(1 − *ε*)) denote the radius and bulk density of the particle (sphere). − ௪,௧ <sup>−</sup> <sup>−</sup> = <sup>௦</sup> (1−))

**Figure 3.** Scheme of mass transfer limited by intraparticle pore diffusion.

For linear sorption with concentration independent distribution coefficients (i.e., *C<sup>s</sup>* = *K<sup>d</sup> Cw*,*intra*) Equation (11) becomes:

$$\frac{\partial \mathbb{C}\_{w,intra}}{\partial t} = D\_d \left[ \frac{\partial^2 \mathbb{C}\_{w,intra}}{\partial r^2} + \frac{2}{r} \frac{\partial \mathbb{C}\_{w,intra}}{\partial r} \right] \tag{14}$$

where *D<sup>a</sup>* [L<sup>2</sup> T −1 ] is the apparent diffusion coefficient, defined as:

$$D\_{\mathfrak{a}} = \frac{D\_{\mathfrak{e}}}{\varepsilon + \mathcal{K}\_d \rho\_p} = \frac{D\_{\mathfrak{a}q}\varepsilon}{\pi \left(\varepsilon + \mathcal{K}\_d \rho\_p\right)} \approx \frac{D\_{\mathfrak{a}q}\varepsilon^2}{\varepsilon + \mathcal{K}\_d (1 - \varepsilon)\rho\_s} \tag{15}$$

Empirical studies showed that *D<sup>e</sup>* 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:

$$\frac{\partial}{\partial t} \left( n\mathbb{C}\_{w} + (1 - n) \left( \varepsilon \mathbb{C}\_{w, \text{intra}} + \rho\_p \mathbb{C}\_s \right) \right) + \frac{\partial}{\partial x} \left( nv\mathbb{C}\_w - nD\_L \frac{\partial \mathbb{C}\_w}{\partial x} \right) = 0 \tag{16}$$

The corresponding equilibrium concentration (*Cw*,*eq*) in water after first flooding can be rewritten as *Cw*,*eq* = *Cs*,*ini*/(*ε*/*ρ<sup>p</sup>* + *K<sup>d</sup>* ), which is slightly lower than for non-porous solids (*Cw*,*eq* = *Cs*,*ini*/*K<sup>d</sup>* ) because the intragranular pore space is assumed to be initially free of water. The deviation becomes insignificant with the increase of *K<sup>d</sup>* (*ε*/*ρ<sup>p</sup>* + *K<sup>d</sup>* ≈ *K<sup>d</sup>* when *K<sup>d</sup>* ≫ *ε*/*ρ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:

$$\frac{\partial \mathbb{C}\_{w}}{\partial t} = D\_{\text{L}} \frac{\partial^{2} \mathbb{C}\_{w}}{\partial \mathbf{x}^{2}} - v \frac{\partial \mathbb{C}\_{w}}{\partial \mathbf{x}} - \frac{(1 - n)}{n} D\_{\text{e}} \left[ \frac{\partial^{2} \mathbb{C}\_{w, \text{intra}}}{\partial r^{2}} + \frac{2}{r} \frac{\partial \mathbb{C}\_{w, \text{intra}}}{\partial r} \right] \tag{17}$$

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.

#### *2.4. 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 (*dp*, *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 1 *n*

Freundlich model (*C<sup>s</sup>* = *Kf rC <sup>w</sup>* where *Kf 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−<sup>1</sup> . Lower *K<sup>d</sup>* values (<0.1 L kg−<sup>1</sup> ) were not considered here (this would

have resulted in very high initial aqueous concentrations). If the *K<sup>d</sup>* values become large (*K<sup>d</sup>* > 100 L kg−<sup>1</sup> ), 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<sup>d</sup>* range chosen covers many frequent environmental contaminants, such as per- and polyfluoroalkyl substances (PFAS), chlorinated solvents, polycyclic aromatic hydrocarbons and some heavy metals.

**Table 1.** Summary of parameter values and ranges used to set up the numerical experiments.


#### **3. Results and Discussion**

#### *3.1. 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 (*Cw*,*eq* = *Cs*,*ini*/ *K<sup>d</sup>* + *<sup>n</sup> ρb* ) and (2) columns with non-uniform concentration distributions after first flooding where concentrations increased towards the outlet (to a maximum of *Cw*,*eq* = *Cs*,*ini*/*K<sup>d</sup>* ) 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<sup>d</sup>* values, the deviation of the initial concentration profiles only occurred at the inlet of the column. At low *K<sup>d</sup>* 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<sup>d</sup>* values.

ௗ

ௗ

ௗ

௨,௫ = 1000 <sup>−</sup>

ௗ

ௗ ௗ ≥ <sup>−</sup>

ௗ

ௗ = <sup>−</sup>

**Figure 4.** Initial aqueous concentration distributions in the column after the first flooding (solid lines) if mass transfer is controlled by film diffusion for three different distribution coefficients, *K<sup>d</sup>* , compared to the pre-equilibrated case (dashed lines). **Top panel**: without dispersion; **bottom panel**: with dispersion; *<sup>n</sup>* = 0.45, *<sup>v</sup>* = 1.67 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m s−<sup>1</sup> , *α*/*x* = 0 or 0.1, *Cs*,*ini* = 1000 µg kg−<sup>1</sup> , *t<sup>c</sup>* = 5 h, *dp*, *f ine* = 63 µm.

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 *<sup>K</sup><sup>d</sup>* and became almost insignificant for *<sup>K</sup><sup>d</sup>* <sup>≥</sup> 10 L kg−<sup>1</sup> . Dispersion also reduces differences between the pre-equilibrated and the first flooding case. At high *K<sup>d</sup>* values, the maximum concentrations were still achieved but the tailings became smoother. With the decrease of the *K<sup>d</sup>* 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<sup>d</sup>* = 0.1 L kg−<sup>1</sup> ). 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<sup>d</sup>* 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<sup>d</sup>* , the earlier the cumulative leachate concentration reaches its maximum μ

value (*mcum*, *max* = 1000 µg kg−<sup>1</sup> ). Dispersion shifts this point to later times (see Figure 5, 3rd row). <sup>−</sup> <sup>−</sup> ௦, μ <sup>−</sup> ,

ௗ

௪/௪,, <sup>௪</sup> ௨ **Figure 5.** Normalized and absolute concentration (*Cw*/*Cw*,*eq*, *Cw*) as well as cumulative concentration (*mcum* ) in the column effluent vs. time (expressed as liquid to solid ratio: *LS*) for the initial conditions (depicted in Figure 4) established after the first flooding of the column (solid lines) compared to the pre-equilibrated case (dashed lines). **Left column**: without dispersion; **right column**: with dispersion.

#### *3.2. 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<sup>d</sup>* 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, *Xs*,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<sup>d</sup>* 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 √ *Kd* (Equation (A32)) (e.g., *Xs*,63.2% = 10 cm, 3.5 cm, and 1.1 cm for *K<sup>d</sup>* values of 0.1 L kg−<sup>1</sup> , 1 L kg−<sup>1</sup> and 10 L kg−<sup>1</sup> (see Figure A2)). The length of the mass transfer zone for IPD is much longer than for FD, but differences decrease with increasing *K<sup>d</sup>* values. For *K<sup>d</sup>* values of 1 L kg−<sup>1</sup> and 10 L kg−<sup>1</sup> , the mass transfer zone lengths for IPD are much shorter than the column length (*Xcol* = 30 cm), which indicates that the equilibrium concentration is achieved at the outlet of the column after the first flooding. For small *K<sup>d</sup>* values (e.g., *K<sup>d</sup>* = 0.1 L kg−<sup>1</sup> ), the equilibrium concentrations are not achieved at the outlet if dispersion is considered (see Figure 6, lower panel) although the mass transfer zone length (*Xs*,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).

<sup>−</sup> <sup>−</sup> ௦, μ <sup>−</sup> <sup>−</sup> <sup>−</sup> ,௦ μ , μ **Figure 6.** Initial aqueous concentration distributions in the column after the first flooding depending on the mass transfer limitation; dotted lines: film diffusion (FD), dashed lines: intraparticle diffusion (IPD); solid lines: fast kinetics (equilibrium, fine particles). **Top panel**: without dispersion; **bottom panel**: with dispersion; *<sup>n</sup>* = 0.45, *<sup>v</sup>* = 1.67 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m s−<sup>1</sup> , *α*/*x* = 0 or 0.1, *Cs*,*ini* = 1000 µg kg−<sup>1</sup> , *<sup>t</sup><sup>c</sup>* = 5 h, *<sup>D</sup>aq* = 1 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup><sup>2</sup> s −1 , *ε* = 0.05, *dp*,*coarse* = 2000 µm, *dp*, *f ine* = 63 µm.

Figure 7 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 (*De*). 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.

௪/௪, ௨ **Figure 7.** Normalized concentrations (*Cw*/*Cw*,*eq*) as well as cumulative concentrations (*mcum* ) in the column effluent vs. time (expressed as liquid to solid ratio: *LS*) for different mass-transfer processes, given the initial conditions depicted in Figure 6. **Left column**: without dispersion; **Right column**: with dispersion.

#### *3.3. 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" 1

ௗ = ௗ/௪, ି ଵ *Kd* for the aqueous concentration at equilibrium: *Kf r* = *Kd*/*C <sup>n</sup>* −1 *<sup>w</sup>*,*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<sup>d</sup>* = 1 L kg−<sup>1</sup> 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<sup>d</sup>* 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.

ௗ <sup>−</sup> 1/ <sup>−</sup> <sup>−</sup> ௦, μ <sup>−</sup> **Figure 8.** Influence of sorption non-linearity: initial aqueous concentration distribution in the column after the first flooding (**left graph**) and column effluent concentration (normalized: **mid graph**, cumulative: **right graph**) vs. time (expressed as liquid to solid ratio: *LS*); solid lines: linear sorption (*K<sup>d</sup>* = 1 L kg−<sup>1</sup> ); dashed lines: non-linear sorption case (Freundlich coefficient *<sup>K</sup>f r* = 7.94, exponent 1/*<sup>n</sup>* = 0.7); *<sup>n</sup>* = 0.45, *<sup>v</sup>* = 1.67 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m s−<sup>1</sup> , *α*/*x* = 0 or 0.1, *Cs*,*ini* = 1000 µg kg−<sup>1</sup> , *t<sup>c</sup>* = 5 h, *<sup>D</sup>aq* = 1 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup><sup>2</sup> s −1 , *ε* = 0.05, *dp*,*coarse* = 2000 µm.

#### <sup>−</sup> <sup>−</sup> ,௦ μ *3.4. 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.

ௗ ௗ = 1/10/100 <sup>−</sup> 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<sup>d</sup>* ) 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<sup>d</sup>* (*Kd*,*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<sup>d</sup>* = 1/10/100 L kg−<sup>1</sup> , 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<sup>d</sup>* 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).

ௗ ௗ,௩ <sup>−</sup> ௗ 10 × ௗ,௩ ௗ 100 × ௗ,௩ ௗ,௩ <sup>−</sup> <sup>−</sup> ௦, μ <sup>−</sup> − <sup>−</sup> ,௦ μ **Figure 9.** Behavior of bi-modal material compositions of sorbing and non-sorbing particles: initial concentration distribution in the column after the first flooding (**top panel**) and column effluent concentration (normalized: **mid panel**, cumulative: **bottom panel**) vs. time (expressed as liquid to solid ratio: LS). Left column: homogeneous case with average *K<sup>d</sup>* (= *Kd*,*av* = 1 L kg−<sup>1</sup> ); mid column: only 10% of the particles carry the contaminant at *K<sup>d</sup>* = 10 × *Kd*,*av*; right column: only 1% of the particles carry the contaminant at *K<sup>d</sup>* = 100 × *Kd*,*av*; the average *Kd*,*av* of the entire material is the same for all compositions; solid lines: film diffusion cases, dashed lines: intraparticle diffusion case; *n* = 0.45, *<sup>v</sup>* = 1.67 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m s−<sup>1</sup> , *α* = 0 (no dispersion), *Cs*,*ini* = 1000 µg kg−<sup>1</sup> , *<sup>t</sup><sup>c</sup>* = 5 h, *<sup>D</sup>aq* = 1 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup><sup>2</sup> s −1 , *ε* = 0.05, *dp*,*coarse* = 2000 µm.

,௦ = μ , = μ 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 *dp*,*coarse* = 2000 µm, and fine particles with *dp*, *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<sup>d</sup>* values (0.1, 1, and 10 L kg−<sup>1</sup> ). While the shapes of all leaching curves are very similar, their locations are shifted in time according to the different *K<sup>d</sup>* 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<sup>d</sup>* values (*K<sup>d</sup>* =0.1 L kg−<sup>1</sup> ). Initial non-equilibrium conditions become more salient for intraparticle pore diffusion if coarse particles predominate (see Figures S5 and S6 in SM). ௗ ௗ = <sup>−</sup>

<sup>−</sup> <sup>−</sup> ௦, μ <sup>−</sup> − <sup>−</sup> ,௦ μ , μ **Figure 10.** Behavior of the bi-modal material compositions of fine and coarse particles: initial concentration distribution in the column after the first flooding (**top panel**) and column effluent concentration (normalized: **mid panel**, cumulative: **bottom panel**) vs. time (expressed as liquid to solid ratio: *LS*); solid lines: fine particle mass fraction 10%; dashed lines: fine particle mass fraction 90%. (*<sup>n</sup>* = 0.45, *<sup>v</sup>* = 1.67 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m s−<sup>1</sup> , *α* = 0 (no dispersion), *Cs*,*ini* = 1000 µg kg−<sup>1</sup> , *t<sup>c</sup>* = 5 h, *<sup>D</sup>aq* = 1 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup><sup>2</sup> s −1 , *ε* = 0.05, *dp*,*coarse* = 2000 µm, *dp*, *f ine* = 63 µm).

ௗ

Combining the heterogeneity of particle size (*d*) and sorption capacity (*K<sup>d</sup>* ), 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<sup>d</sup>* vs. just a few large particles carrying the contaminant with a large *K<sup>d</sup>* . A hypothetical mixture containing 10% of fine particles with low sorption capacity (*K<sup>d</sup>* = 10 L kg−<sup>1</sup> ) and 90% of coarse particles with high sorption capacity (*K<sup>d</sup>* = 100 L kg−<sup>1</sup> ) 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<sup>d</sup>* or only coarse particles with a high *K<sup>d</sup>* . Pure coarse material with a high *K<sup>d</sup>* shows a low equilibrium concentration (*Cw*,*eq* = *Cs*,*ini*/*K<sup>d</sup>* = 1000 µg kg−1/100 L kg−<sup>1</sup> = 10 µg L−<sup>1</sup> ) while pure fine material with a low *K<sup>d</sup>* presents a much higher equilibrium concentration (*Cw*,*eq* = *Cs*,*ini*/*K<sup>d</sup>* = 1000 µg kg−1/10 L kg−<sup>1</sup> = 100 µg L−<sup>1</sup> ) 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<sup>c</sup>* = 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).

ௗ = <sup>−</sup> ௗ = <sup>−</sup> <sup>−</sup> <sup>−</sup> ௦, μ <sup>−</sup> <sup>−</sup> <sup>−</sup> ,௦ μ , μ **Figure 11.** Behavior of bi-modal material compositions of fine particles with low sorption capacity (*K<sup>d</sup>* =10 L kg−<sup>1</sup> ) and coarse particles with high sorption capacity (*K<sup>d</sup>* = 100 L kg−<sup>1</sup> ): initial concentration distributions in the column after the first flooding (**top panel**) and the column effluent concentration (normalized: **mid panel**, cumulative: **bottom panel**) vs. time (expressed as liquid to solid ratio: *LS*). Left column: 100% coarse particles; mid column: mixed sample with 10% fine particles; right column: 100% fine particles; *<sup>n</sup>* = 0.45, *<sup>v</sup>* = 1.67 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m s−<sup>1</sup> , *α* = 0 (no dispersion), *Cs*,*ini* = 1000 µg kg−<sup>1</sup> , *<sup>t</sup><sup>c</sup>* = 5 h, *<sup>D</sup>aq* = 1 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup><sup>2</sup> s −1 , *ε* = 0.05, *dp*,*coarse* = 2000 µm, *dp*, *f ine* = 63 µm.

#### **4. Summary and Conclusions**

ௗ = 0.1 100 <sup>−</sup> We conducted numerical simulations to investigate the release characteristics of low to strongly sorbing compounds (*K<sup>d</sup>* = 0.1–100 L kg−<sup>1</sup> ) 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<sup>d</sup>* values (*K<sup>d</sup>* < 1 L kg−<sup>1</sup> ). With increasing *K<sup>d</sup>* , the differences between pre-equilibrium and non-equilibrium conditions gradually vanish for *K<sup>d</sup>* > 10 L kg−<sup>1</sup> (see Figure 5). Compounds with very low *K<sup>d</sup>* ("salts") would reach extremely high concentrations (*K<sup>d</sup>* << 1 L kg−<sup>1</sup> ) at the column outlet (see Figure 4) potentially leading to enhanced dispersion due to density fingering. The *K<sup>d</sup>* values derived from retardation factors (*R<sup>d</sup>* 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 (*Cw*/*Cw*,*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<sup>d</sup>* 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<sup>d</sup>* value is calculated for the proper concentration range (since for the nonlinear sorption the *K<sup>d</sup>* depends on the concentration, large deviations may occur if just the *Kf r* is determined far away of the sample's concentration is used as "*Kd*"); 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<sup>d</sup>* = 100 × *Kd*,*av*), leaching may start at a plateau (suggesting equilibrium), but far below equilibrium concentrations (*Cw*,*peak* ≪ *Cw*,*eq*) and concentrations later decrease further; The *K<sup>d</sup>* values derived from the initial aqueous concentration (*K<sup>d</sup>* = *Cs*,*ini*/*Cw*,*peak*) would be overestimated while the *K<sup>d</sup>* 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<sup>d</sup>* =10 L kg−<sup>1</sup> ) and large amounts of coarse particles with high sorption capacity (*K<sup>d</sup>* = 100 L kg−<sup>1</sup> ), 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 (*Cw*,*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<sup>d</sup>* values derived from the initial aqueous concentration (*K<sup>d</sup>* = *Cs*,*ini*/*Cw*,*peak*) would be underestimated, while *K<sup>d</sup>* 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).

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/ma14164708/s1, Figure S1. Initial concentration profiles in the column after the first flooding (up-flow) without and with dispersion (top and bottom panel); solid lines: linear sorption; dashed lines: non-linear sorption cases (based on a Freundlich exponent 1/*n* = 0.7), *n* = 0.45, *<sup>v</sup>* = 1.67 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m s−<sup>1</sup> , *α*/*x* = 0 or 0.1, *Cs*,*ini* = 1000 µg kg−<sup>1</sup> , *<sup>t</sup><sup>c</sup>* = 5 h, *<sup>D</sup>aq* = 1 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup><sup>2</sup> s −1 , *ε* = 0.05, *dp*,*coarse* = 2000 µm, Figure S2. Normalized concentrations (*Cw*/*Cw*,*eq*) as well as cumulative concentrations (*mcum*) in the column effluent vs. time (expressed as liquid to solid ratio: *LS*) for different initial conditions depicted in Figure S1; solid lines: linear sorption; dashed lines: nonlinear sorption. Left column: without dispersion; right column: with dispersion, Figure S3. Initial concentration distribution in the column after the first flooding (up-flow) for different bi-modal compositions of sorbing and non-sorbing particles; left column: homogeneous case with average *K<sup>d</sup>* (=*Kd*,*av* = 1 L kg−<sup>1</sup> ); mid column: only 10% of the particles carry the contaminant at *K<sup>d</sup>* = 10 × *Kd*,*av*; right column: only 1% of the particles carry the contaminant at *K<sup>d</sup>* = 100 × *Kd*,*av*; the average *Kd*,*av* of the entire material is the same for all compositions; solid lines: film diffusion case, dashed lines: intraparticle diffusion case. Top panel: without dispersion; bottom panel: with dispersion; *<sup>n</sup>* = 0.45, *<sup>v</sup>* = 1.67 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m s−<sup>1</sup> , *α*/*x* = 0 or 0.1, *Cs*,*ini* = 1000 µg kg−<sup>1</sup> , *<sup>t</sup><sup>c</sup>* = 5 h, *<sup>D</sup>aq* = 1 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup><sup>2</sup> s −1 , *ε* = 0.05, *dp*,*coarse* = 2000 µm, Figure S4. Normalized concentrations (*Cw*/*Cw*,*eq*) as well as cumulative concentrations (*mcum*) 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%. (*<sup>n</sup>* = 0.45, *<sup>v</sup>* = 1.67 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m s−<sup>1</sup> , *α*/*x* = 0 or 0.1, *Cs*,*ini* = 1000 µg kg−<sup>1</sup> , *<sup>t</sup><sup>c</sup>* = 5 h, *<sup>D</sup>aq* = 1 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup><sup>2</sup> s −1 , *ε* = 0.05, *dp*,*coarse* = 2000 µm, *dp*, *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 (*Cw*/*Cw*,*eq*) as well as cumulative concentrations (*mcum*) 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<sup>d</sup>* = 10 L kg−<sup>1</sup> ) and coarse particles with high sorption capacity; left: 100% coarse particles (*K<sup>d</sup>* = 100 L kg−<sup>1</sup> ); middle: mixed sample with 10% fine particles; right: 100% fine particles; solid lines: film diffusion (FD), dashed lines: intraparticle diffusion cases (IPD); *<sup>n</sup>* = 0.45, *<sup>v</sup>* = 1.67 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m s−<sup>1</sup> , *α*/*x* = 0 or 0.1, *Cs*,*ini* = 1000 µg kg−<sup>1</sup> , *<sup>t</sup><sup>c</sup>* = 5 h, *<sup>D</sup>aq* = 1 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup><sup>2</sup> s −1 , *ε* = 0.05, *dp*,*coarse* = 2000 µm, *dp*, *f ine* = 63 µm; top panel: without dispersion; bottom panel: with dispersion, Figure S8. Leachate concentrations (*Cw*) as well as cumulative concentrations (*mcum*) in the column effluent vs. time (expressed as liquid to solid ratio: *LS*) for different combinations of fine particles with low sorption capacity (*K<sup>d</sup>* = 10 L kg−<sup>1</sup> ) and coarse particles with high sorption capacity (*K<sup>d</sup>* = 100 L kg−<sup>1</sup> ); left: without dispersion; right: with dispersion; solid lines: film diffusion cases, dashed lines: intraparticle diffusion cases; kinetic parameters are the same as Figure S7.

**Author Contributions:** Conceptualization, B.L., M.F., P.G.; validation, B.L., M.F., P.G.; investigation, B.L., M.F., P.G.; writing—original draft preparation, B.L.; writing—review and editing, B.L., M.F., P.G.; visualization, B.L.; supervision, M.F., P.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partly supported by the Collaborative Research Center 1253 CAMPOS, funded by the German Research Foundation (DFG, Grant Agreement SFB 1253/1 2017).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent 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:

$$Sh = A + BRe^{\theta} Sc^{\gamma} \tag{A1}$$

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:

$$Sc = \frac{\eta}{\rho\_{\rm L} D\_{aq}} \quad Re = \frac{\rho\_{\rm L} d \, v}{\eta} \tag{A2}$$

where *η* [M L−<sup>1</sup> T −1 ] denotes the dynamic viscosity of the fluid. *ρ<sup>L</sup>* [M L−<sup>3</sup> ] is the density of the fluid and *v* [L T−<sup>1</sup> ] denotes the flow velocity.

The empirical exponents *θ* and *γ* in Equation (A1) may be determined experimentally or from theory. The Blasius (1908) solved the Navier-Stokes equation and continuity equation for laminar flow over a sharp leading edge and found that the ratio of fluid velocity boundary layer thickness to concentration boundary layer thickness is proportional to the Schmidt number with a power of 1/3 (= *γ* in Equation (A1)) which is widely used in literature [37–43]. Liu et al. (2014) showed a higher empirical exponent *γ* of 1/2 based on penetration theory [35,44]. *θ* values depend on the experimental setup and are generally adapted from experimental data. Most of the empirical relationships show that *θ* values lie in the range of 0.5–0.75 [35,37,39–43].

Liu et al. (2014) proposed an empirical relationship for mass transfer in packed beds only based on the Peclet number (*Pe* = *Re* × *Sc*) [35]:

$$Sh = 2 + 0.1Pe^{1/2} \tag{A3}$$

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:

$$\begin{array}{l} \frac{\partial \mathbb{C}\_{w,j}}{\partial t} = D\_L \frac{(\mathbb{C}\_{w,j-1} - \mathbb{C}\_{w,j} + \mathbb{C}\_{w,j+1})}{\Delta x^2} - \upsilon \frac{(\mathbb{C}\_{w,j} - \mathbb{C}\_{w,j-1})}{\Delta x} \\ + \frac{D\_{aq}}{\delta\_f} \frac{6(1-n)}{nd} \left(\frac{\mathbb{C}\_{s,j}}{\mathbb{K}\_d} - \mathbb{C}\_{w,j}\right) \\ \frac{\partial \mathbb{C}\_{s,j}}{\partial t} = -\frac{D\_{aq}}{\delta\_f} \frac{6}{\rho\_s \, d} \left(\frac{\mathbb{C}\_{s,j}}{\mathbb{K}\_d} - \mathbb{C}\_{w,j}\right) \end{array} \tag{A4}$$

where *Cw*,*<sup>j</sup>* [M L−<sup>3</sup> ], *Cw*,*j*−<sup>1</sup> [M L−<sup>3</sup> ] and *Cw*,*j*+<sup>1</sup> [M L−<sup>3</sup> ] denote the solute concentration in the water phase in volume *j*, *j* − 1 and *j* + 1, respectively. *Cs*,*<sup>j</sup>* [M M−<sup>1</sup> ] denotes the solute concentration in the solid phase in volume *j*.

**Figure A1.** Discretization of the column into *N* parts (**a**). Representation of the solid phase as a composition of grains having different sizes and properties, each discretized by a number of *L* shells (**b**).

(, ାଵ ቈ + ൫, ାଵ൯ ଵ ିଵ + ∆ ଶ 1− ା.ହ − ି.ହ <sup>ቆ</sup> ା.ହ ଶ ାଵ − + ି.ହ ଶ − ିଵቇ , ାଵ 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.

= ቈ∆ ଶ 1− ା.ହ − ି.ହ ା.ହ ଶ ାଵ − ାଵ, ାଵ + ቈ∆ ଶ 1− ା.ହ − ି.ହ ି.ହ ଶ − ିଵ ିଵ, ାଵ + ቈ∆ ଶ ା.ହ − ି.ହ ା.ହ ଶ ାଵ − ାଵ, + ቈ∆ ଶ ା.ହ − ି.ହ ି.ହ ଶ − ିଵ ିଵ, + ቈ + ൫, ൯ ଵ ିଵ − ∆ ଶ ା.ହ − ି.ହ <sup>ቆ</sup> ା.ହ ଶ ାଵ − + ି.ହ ଶ − ିଵቇ , =2 =−1 nd *C k*+1 *<sup>w</sup>*,*<sup>j</sup>* −*C k w*,*j* <sup>∆</sup>*<sup>t</sup>* = (1 − *ϕ*) *D<sup>L</sup> C k*+1 *<sup>w</sup>*,*j*−1−2*<sup>C</sup> k*+1 *<sup>w</sup>*,*<sup>j</sup>* +*C k*+1 *w*,*j*+1 ∆*x* <sup>2</sup> − *v C k*+1 *<sup>w</sup>*,*<sup>j</sup>* −*C k*+1 *w*,*j*−1 <sup>∆</sup>*<sup>x</sup>* + *Daq δf* 6(1−*n*) *nd C k*+1 *s*,*j Kd* − *C k*+1 *w*,*j* ! +*ϕ D<sup>L</sup> C k <sup>w</sup>*,*j*−1−2*<sup>C</sup> k <sup>w</sup>*,*j*+*C k w*,*j*+1 ∆*x* <sup>2</sup> − *v C k <sup>w</sup>*,*j*−*C k w*,*j*−1 <sup>∆</sup>*<sup>x</sup>* + *Daq δf* 6(1−*n*) *nd C k s*,*j Kd* − *C k w*,*j* ! *C k*+1 *<sup>s</sup>*,*<sup>j</sup>* −*C k s*,*j* <sup>∆</sup>*<sup>t</sup>* = −(1 − *ϕ*) *Daq δf* 6 *ρs d C k*+1 *s*,*j Kd* − *C k*+1 *w*,*j* −*ϕ Daq δf* 6 *ρs d C k s*,*j Kd* − *C k w*,*j* (A5)

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* = [*Cw*; *Cs*]; with the semicolon being a line delimiter):

$$\mathbf{C} = \begin{bmatrix} \mathbf{C}\_{w,1} \\ \mathbf{C}\_{s,1} \\ \vdots \\ \mathbf{C}\_{w,j} \\ \mathbf{C}\_{s,j} \\ \vdots \\ \mathbf{C}\_{w,N} \\ \mathbf{C}\_{s,N} \end{bmatrix}\_{2N \times 1} = \begin{bmatrix} \mathbf{C}\_{1} \\ \mathbf{C}\_{2} \\ \vdots \\ \mathbf{C}\_{2j-1} \\ \mathbf{C}\_{2j} \\ \vdots \\ \mathbf{C}\_{2N-1} \\ \mathbf{C}\_{2N} \end{bmatrix}\_{2N \times 1} \tag{A6}$$

with *jǫ*[1, 2, . . . , N].

1 ∆ ቆ,

= ଶ

+ ଶ ୩ାଵ + ൫,

1− ା.ହ − ି.ହ ቆା.ହ

 ା.ହ − ି.ହ ቆା.ହ

ାଵ൯ ଵ − ,

> ାଵ, ାଵ − ,

ାଵ, − , 

ାଵ −

ାଵ −

ଶ

ଶ

 + 0.5 − 0.5 +1 −1

− ൫,

ାଵ

 ൯ ଵ ቇ

,

,

ାଵ − ିଵ, ାଵ − ିଵ <sup>ቇ</sup>

 − ିଵ, − ିଵ <sup>ቇ</sup>

+1

− ି.ହ ଶ

− ି.ହ ଶ

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:

*<sup>f</sup>*2*j*−<sup>1</sup> = *C k*+1 *<sup>w</sup>*,*<sup>j</sup>* −*C k w*,*j* ∆*t* −(1 − *ϕ*) *D<sup>L</sup> C k*+1 *<sup>w</sup>*,*j*−1−2*<sup>C</sup> k*+1 *<sup>w</sup>*,*<sup>j</sup>* +*C k*+1 *w*,*j*+1 ∆*x* <sup>2</sup> − *v C k*+1 *<sup>w</sup>*,*<sup>j</sup>* −*C k*+1 *w*,*j*−1 ∆*x* + *Daq δf* 6(1−*n*) *nd C k*+1 *s*,*j Kd* − *C k*+1 *w*,*j* −*ϕ D<sup>L</sup> C k <sup>w</sup>*,*j*−1−2*<sup>C</sup> k <sup>w</sup>*,*j*+*C k w*,*j*+1 ∆*x* <sup>2</sup> − *v C k <sup>w</sup>*,*j*−*C k w*,*j*−1 ∆*x* + *Daq δf* 6(1−*n*) *nd C k s*,*j Kd* − *C k w*,*j f*2*<sup>j</sup>* = *C k*+1 *<sup>s</sup>*,*<sup>j</sup>* −*C k s*,*j* <sup>∆</sup>*<sup>t</sup>* +(1 − *ϕ*) *Daq δf* 6 *ρs d C k*+1 *s*,*j Kd* − *C k*+1 *w*,*j* +*ϕ Daq δf* 6 *ρs d C k s*,*j Kd* − *C k w*,*j* (A7)

The residual function vector can be expressed as:

$$f\left(\mathbb{C}^{k+1}\right) = \begin{bmatrix} f\_1 \\ f\_2 \\ \vdots \\ f\_{2j-1} \\ f\_{2j} \\ \vdots \\ f\_{2N-1} \\ f\_{2N} \end{bmatrix}\_{2N \times 1} \tag{A8}$$

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:

$$\begin{split} f\left(\mathsf{C}^{k+1}\right) &\approx f\left(\mathsf{C}^{k+1}\_{\text{gauss}}\right) + \left(\mathsf{C}^{k+1} - \mathsf{C}^{k+1}\_{\text{gauss}}\right) \frac{\partial f\left(\mathsf{C}^{k+1}\right)}{\partial \mathsf{C}^{k+1}}\Big|\_{\mathsf{C}^{k+1} = \mathsf{C}^{k+1}\_{\text{gauss}}} \\ \mathsf{C}^{k+1} &= \mathsf{C}^{k+1}\_{\text{gauss}} - \frac{f\left(\mathsf{C}^{k+1}\_{\text{gress}}\right)}{f} = \mathsf{C}^{k+1}\_{\text{gauss}} - \frac{f\left(\mathsf{C}^{k+1}\_{\text{gress}}\right)}{\frac{\partial f\left(\mathsf{C}^{k+1}\right)}{\partial \mathsf{C}^{k+1}}\Big|\_{\mathsf{C}^{k+1} = \mathsf{C}^{k+1}\_{\text{gress}}}} \end{split} \tag{A9}$$

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 <sup>k</sup>*+<sup>1</sup> and Equation (A9) is reapplied. In our case, the Jacobian matrix can be derived analytically:

$$J = \begin{bmatrix} \frac{\partial f\_1}{\partial \mathbf{C}\_1} & \frac{\partial f\_1}{\partial \mathbf{C}\_2} & \cdots & \frac{\partial f\_1}{\partial \mathbf{C}\_{2N-1}} & \frac{\partial f\_1}{\partial \mathbf{C}\_{2N}}\\ \frac{\partial f\_2}{\partial \mathbf{C}\_1} & \frac{\partial f\_2}{\partial \mathbf{C}\_2} & \cdots & \frac{\partial f\_2}{\partial \mathbf{C}\_{2N-1}} & \frac{\partial f\_2}{\partial \mathbf{C}\_{2N}}\\ \vdots & \vdots & \cdots & \vdots & \vdots\\ \frac{\partial f\_{2N-1}}{\partial \mathbf{C}\_1} & \frac{\partial f\_{2N-1}}{\partial \mathbf{C}\_2} & \cdots & \frac{\partial f\_{2N-1}}{\partial \mathbf{C}\_{2N-1}} & \frac{\partial f\_{2N-1}}{\partial \mathbf{C}\_{2N}}\\ \frac{\partial f\_{2N}}{\partial \mathbf{C}\_1} & \frac{\partial f\_{2N}}{\partial \mathbf{C}\_2} & \cdots & \frac{\partial f\_{2N}}{\partial \mathbf{C}\_{2N-1}} & \frac{\partial f\_{2N}}{\partial \mathbf{C}\_{2N}} \end{bmatrix} \tag{A10}$$

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:

$$
\Delta \mathbf{C}^{k+1} = \begin{bmatrix}
\begin{bmatrix}
\mathbf{c}\_{1,\text{gauss},\text{new}}^{k+1} - \mathbf{C}\_{1,\text{gases},\text{old}}^{k+1} \\
\mathbf{c}\_{2,\text{gases},\text{new}}^{k+1} - \mathbf{C}\_{2,\text{gases},\text{old}}^{k+1} \\
\vdots \\
\mathbf{c}\_{2,\text{j.},\text{gases},\text{new}}^{k+1} - \mathbf{C}\_{2,\text{j.},\text{gases},\text{old}}^{k+1} \\
\mathbf{c}\_{2,\text{j.},\text{gases},\text{new}}^{k+1} - \mathbf{C}\_{2,\text{j.},\text{gases},\text{old}}^{k+1}
\end{bmatrix}
\tag{A11}
$$

$$
\begin{bmatrix}
\mathbf{c}\_{2N-1,\text{gases},\text{new}}^{k+1} - \mathbf{C}\_{2N-1,\text{gases},\text{old}}^{k+1} \\
\vdots \\
\mathbf{c}\_{2N-1,\text{gases},\text{new}}^{k+1} - \mathbf{C}\_{2N,\text{gases},\text{old}}^{k+1}
\end{bmatrix}\_{2N \times 1}
$$

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).

#### **Appendix C. Intraparticle Pore Diffusion Coupled to Advective-Dispersive Transport**

Intraparticle pore diffusion is widely used to describe the sorptive uptake of pollutants in porous materials such as activated carbon, zeolites and many technical materials. Equations (11) and (17) describe intraparticle diffusion coupled to advective-dispersive transport. The intraparticle diffusion model approximates the solid grains as spherical particles. These spherical particles are discretized into a number of shells of equal volume. Mass transfer between solid and intra-granular water phases is assumed to be fast and local equilibrium is assumed. For sorption, the Freundlich isotherm model is employed for nonlinear and linear (exponent = 1) cases. Figure A1b shows the numerical grain model where the spherical grains are divided into *L* shells. For a specific shell *p* in volume *j* the corresponding difference-equations were used [46]:

$$\begin{split} &\frac{1}{\Delta t} \left( \varepsilon \mathbf{C}\_{p,j}^{k+1} + \rho\_p \mathbf{K}\_{fr} \left( \mathbf{C}\_{p,j}^{k+1} \right)^{\frac{1}{n}} - \varepsilon \mathbf{C}\_{p,j}^{k} - \rho\_p \mathbf{K}\_{fr} \left( \mathbf{C}\_{p,j}^{k} \right)^{\frac{1}{n}} \right) \\ &= \frac{D\_\varepsilon}{r\_p^2} \frac{1-\varrho}{r\_{p+0.5}-r\_{p-0.5}} \left( r\_{p+0.5}^2 \frac{\mathbf{C}\_{p+1,j}^{k+1} - \mathbf{C}\_{p,j}^{k+1}}{r\_{p+1}-r\_p} - r\_{p-0.5}^2 \frac{\mathbf{C}\_{p,j}^{k+1} - \mathbf{C}\_{p-1,j}^{k+1}}{r\_p - r\_{p-1}} \right) \\ &+ \frac{D\_\varepsilon}{r\_p^2} \frac{\varrho}{r\_{p+0.5}-r\_{p-0.5}} \left( r\_{p+0.5}^2 \frac{\mathbf{C}\_{p+1,j}^k - \mathbf{C}\_{p,j}^k}{r\_{p+1}-r\_p} - r\_{p-0.5}^2 \frac{\mathbf{C}\_{p,j}^k - \mathbf{C}\_{p-1,j}^k}{r\_p - r\_{p-1}} \right) \end{split} \tag{A12}$$

where the subscripts *p* + 0.5 and *p* − 0.5 represent the corresponding parameter value between shells (*p* and *p* + 1) and (*p* and *p* − 1), respectively. Subscript *j* denotes the corresponding parameter value in volume *j*. Subscripts *k* and *k* + 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 (*C k*+1 *p*,*j* ) can be expressed as:

$$\begin{split} & \left[ \varepsilon + \rho\_{p} K\_{fr} \left( C\_{p,j}^{k+1} \right)^{\frac{1}{h} - 1} + \frac{D\_{t} \Delta t}{r\_{p}^{2}} \frac{1 - \rho}{r\_{p+0.5} - r\_{p-0.5}} \left( \frac{r\_{p+0.5}^{2}}{r\_{p+1} - r\_{p}} + \frac{r\_{p-0.5}^{2}}{r\_{p} - r\_{p-1}} \right) \right] \mathbf{C}\_{p,j}^{k+1} \\ &= \left[ \frac{D\_{c} \Delta t}{r\_{p}^{2}} \frac{1 - \rho}{r\_{p+0.5} - r\_{p-0.5}} \frac{r\_{p+0.5}^{2}}{r\_{p+1} - r\_{p-1}} \right] \mathbf{C}\_{p+1,j}^{k+1} + \left[ \frac{D\_{c} \Delta t}{r\_{p}^{2}} \frac{1 - \rho}{r\_{p+0.5} - r\_{p-0.5}} \frac{r\_{p-0.5}^{2}}{r\_{p} - r\_{p-1}} \right] \mathbf{C}\_{p-1,j}^{k+1} \\ &+ \left[ \frac{D\_{c} \Delta t}{r\_{p}^{2}} \frac{q}{r\_{p+0.5} - r\_{p-0.5}} \frac{r\_{p+0.5}^{2}}{r\_{p+1} - r\_{p-1}} \right] \mathbf{C}\_{p+1,j}^{k} + \left[ \frac{D\_{c} \Delta t}{r\_{p}^{2}} \frac{q}{r\_{p+0.5} - r\_{p-0.5}} \frac{r\_{p-0.5}^{2}}{r\_{p} - r\_{p-1}} \right] \mathbf{C}\_{p-1,j}^{k} \\ &+ \left[ \varepsilon + \rho\_{p} K\_{fr} \left( \mathbf{C}\_{p,j}^{k} \right) \right$$

for shell *p* = 2 to shell *p* = *L* − 1 and

1 ∆*t εC k*+1 1,*<sup>j</sup>* <sup>+</sup> *<sup>ρ</sup>pKf r C k*+1 1,*j* 1 *n* − *εC k* 1,*<sup>j</sup>* <sup>−</sup> *<sup>ρ</sup>pKf r C k* 1,*j* 1 *n* = *De r* 2 1 1−*ϕ <sup>r</sup>*1+0.5−*r*1−0.5 *r* 2 1+0.5 *C k*+1 2,*<sup>j</sup>* −*C k*+1 1,*j r*2−*r*<sup>1</sup> + *De r* 2 1 *ϕ <sup>r</sup>*1+0.5−*r*1−0.5 *r* 2 1+0.5 *C k* 2,*j*−*C k* 1,*j r*2−*r*<sup>1</sup> After transformation : *<sup>ε</sup>* <sup>+</sup> *<sup>ρ</sup>pKf r C k*+1 1,*j* 1 *<sup>n</sup>* −1 + *De*∆*t r* 2 1 1−*ϕ <sup>r</sup>*1+0.5−*r*1−0.5 *r* 2 1+0.5 *r*2−*r*<sup>1</sup> *C k*+1 1,*j* = *De*∆*t r* 2 1 1−*ϕ <sup>r</sup>*1+0.5−*r*1−0.5 *r* 2 1+0.5 *r*2−*r*<sup>1</sup> *C k*+1 2,*<sup>j</sup>* + *De*∆*t r* 2 1 *ϕ <sup>r</sup>*1+0.5−*r*1−0.5 *r* 2 1+0.5 *r*2−*r*<sup>1</sup> *C k* 2,*j* + *<sup>ε</sup>* <sup>+</sup> *<sup>ρ</sup>pKf r C k* 1,*j* 1 *<sup>n</sup>* −1 − *De*∆*t r* 2 1 *ϕ <sup>r</sup>*1+0.5−*r*1−0.5 *r* 2 1+0.5 *r*2−*r*<sup>1</sup> *C k* 1,*j* (A14)

for shell 1 (or the innermost shell, *p* = 1) and:

1 ∆*t εC k*+1 *<sup>L</sup>*,*<sup>j</sup>* <sup>+</sup> *<sup>ρ</sup>pKf r C k*+1 *L*,*j* 1 *n* − *εC k <sup>L</sup>*,*<sup>j</sup>* <sup>−</sup> *<sup>ρ</sup>pKf r C k L*,*j* 1 *n* = *De r* 2 *L* 1−*ϕ <sup>R</sup>*−*rL*−0.5 *R* 2 *C k*+1 *<sup>w</sup>*,*<sup>j</sup>* −*C k*+1 *L*,*j R*−*r<sup>L</sup>* − *r* 2 *L*−0.5 *C k*+1 *<sup>L</sup>*,*<sup>j</sup>* −*C k*+1 *L*−1,*j <sup>r</sup>L*−*rL*−<sup>1</sup> + *De r* 2 *L ϕ <sup>R</sup>*−*rL*−0.5 *R* 2 *C k <sup>w</sup>*,*j*−*C k L*,*j R*−*r<sup>L</sup>* − *r* 2 *L*−0.5 *C k <sup>L</sup>*,*j*−*C k L*−1,*j <sup>r</sup>L*−*rL*−<sup>1</sup> After transformation : *<sup>ε</sup>* <sup>+</sup> *<sup>ρ</sup>pKf r C k*+1 *L*,*j* 1 *<sup>n</sup>* −1 + *De*∆*t r* 2 *L* 1−*ϕ <sup>R</sup>*−*rL*−0.5 *R* 2 *R*−*r<sup>L</sup>* + *r* 2 *L*−0.5 *<sup>r</sup>L*−*rL*−<sup>1</sup> *<sup>C</sup> k*+1 *L*,*j* = *De*∆*t r* 2 *L* 1−*ϕ <sup>R</sup>*−*rL*−0.5 *R* 2 *R*−*r<sup>L</sup> C k*+1 *<sup>w</sup>*,*<sup>j</sup>* + *De*∆*t r* 2 *L* 1−*ϕ <sup>R</sup>*−*rL*−0.5 *r* 2 *L*−0.5 *<sup>r</sup>L*−*rL*−<sup>1</sup> *C k*+1 *L*−1,*j* + *De*∆*t r* 2 *L ϕ <sup>R</sup>*−*rL*−0.5 *R* 2 *R*−*r<sup>L</sup> C k <sup>w</sup>*,*<sup>j</sup>* + *De*∆*t r* 2 *L ϕ <sup>R</sup>*−*rL*−0.5 *r* 2 *L*−0.5 *<sup>r</sup>L*−*rL*−<sup>1</sup> *C k L*−1,*j* + *<sup>ε</sup>* <sup>+</sup> *<sup>ρ</sup>pKf r C k L*,*j* 1 *<sup>n</sup>* −1 − *De*∆*t r* 2 *L ϕ <sup>R</sup>*−*rL*−0.5 *R* 2 *R*−*r<sup>L</sup>* + *r* 2 *L*−0.5 *<sup>r</sup>L*−*rL*−<sup>1</sup> *<sup>C</sup> k L*,*j* (A15)

for shell *L* (or the outermost shell, *p* = *L*).

Based on the mass balance, solute mass change in the external water phase (*Mw*) equals the solute mass change in the spherical particles; for better understanding, the simple case of particles with uniform size is shown:

$$\frac{\partial M\_w}{\partial t} = V\_w \frac{\partial \mathcal{C}\_w}{\partial t} = 4\pi R^2 F N\_p \tag{A16}$$

where *F* [M L−<sup>2</sup> T −1 ] denotes the solute flux density into the external water phase. *R* and *N<sup>p</sup>* denote the radius and the total number of the spherical particles. The latter can be calculated by:

$$N\_p = \frac{m\_d}{\rho\_p \left(\frac{4}{3}\pi R^3\right)}\tag{A17}$$

The solute flux density into the external water phase is given by:

$$F = D\_\varepsilon \frac{\mathbf{C}\_L - \mathbf{C}\_w}{R - r\_L} \tag{A18}$$

Substituting *F* and *N<sup>p</sup>* 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:

*C k*+1 *<sup>w</sup>*,*<sup>j</sup>* −*C k w*,*j* <sup>∆</sup>*<sup>t</sup>* = (1 − *ϕ*) *D<sup>L</sup> C k*+1 *<sup>w</sup>*,*j*−1−2*<sup>C</sup> k*+1 *<sup>w</sup>*,*<sup>j</sup>* +*C k*+1 *w*,*j*+1 ∆*x* <sup>2</sup> − *v C k*+1 *<sup>w</sup>*,*<sup>j</sup>* −*C k*+1 *w*,*j*−1 ∆*x* + 3*Dem<sup>d</sup> ρpVwR C k*+1 *<sup>L</sup>*,*<sup>j</sup>* −*C k*+1 *w*,*j R*−*r<sup>L</sup>* +*ϕ D<sup>L</sup> C k <sup>w</sup>*,*j*−1−2*<sup>C</sup> k <sup>w</sup>*,*j*+*C k w*,*j*+1 ∆*x* <sup>2</sup> − *v C k <sup>w</sup>*,*j*−*C k w*,*j*−1 ∆*x* + 3*Dem<sup>d</sup> ρpVwR C k <sup>L</sup>*,*j*−*C k w*,*j R*−*r<sup>L</sup>* After transformation : [1 + (1 − *ϕ*) 3*Demd*∆*t <sup>ρ</sup>pVwR*(*R*−*rL*) <sup>+</sup> 2*DL*∆*t* ∆*x* <sup>2</sup> + *<sup>v</sup>*∆*<sup>t</sup>* ∆*x* i*C k*+1 *w*,*j* = h (1 − *ϕ*) *DL*∆*t* ∆*x* <sup>2</sup> + *<sup>v</sup>*∆*<sup>t</sup>* ∆*x* i*C k*+1 *<sup>w</sup>*,*j*−<sup>1</sup> <sup>+</sup> h (1−*ϕ*)*DL*∆*<sup>t</sup>* ∆*x* 2 i *C k*+1 *w*,*j*+1 + h (1−*ϕ*)3*Demd*∆*<sup>t</sup> ρpVwR*(*R*−*rL*) i *C k*+1 *<sup>L</sup>*,*<sup>j</sup>* + h *ϕ DL*∆*t* ∆*x* <sup>2</sup> + *<sup>v</sup>*∆*<sup>t</sup>* ∆*x* i*C k w*,*j*−1 + h *ϕDL*∆*t* ∆*x* 2 i *C k <sup>w</sup>*,*j*+<sup>1</sup> + h *ϕ*3*Demd*∆*t ρpVwR*(*R*−*rL*) i *C k L*,*j* + h 1 − *ϕ* 3*Demd*∆*t <sup>ρ</sup>pVwR*(*R*−*rL*) <sup>+</sup> 2*DL*∆*t* ∆*x* <sup>2</sup> + *<sup>v</sup>*∆*<sup>t</sup>* ∆*x* i*C k w*,*j* (A19)

In order to solve this system of equations, we may merge the two concentration vectors to a single one (*C* = [*Cw*; *Cp*]; with the semicolon being a line delimiter):

$$\mathbf{C} = \begin{bmatrix} \mathbf{C}\_{w,1} \\ \mathbf{C}\_{1,1} \\ \vdots \\ \mathbf{C}\_{p,1} \\ \mathbf{C}\_{p,1} \\ \vdots \\ \mathbf{C}\_{w,j} \\ \mathbf{C}\_{1,j} \\ \vdots \\ \mathbf{C}\_{p,j} \\ \mathbf{C}\_{p,j} \\ \vdots \\ \mathbf{C}\_{L,j} \\ \vdots \\ \mathbf{C}\_{L,N} \\ \mathbf{C}\_{p,N} \\ \vdots \\ \mathbf{C}\_{p,N} \\ \vdots \\ \mathbf{C}\_{p,N} \\ \vdots \\ \mathbf{C}\_{p,N} \\ \end{bmatrix} = \begin{bmatrix} \mathbf{C}\_{1} \\ \vdots \\ \mathbf{C}\_{p,1} \\ \vdots \\ \mathbf{C}\_{L,1} \\ \vdots \\ \mathbf{C}\_{(L+1)\*(j-1)+2} \\ \vdots \\ \mathbf{C}\_{(L+1)\*(j-1)+p+1} \\ \vdots \\ \mathbf{C}\_{(L+1)\*j} \\ \vdots \\ \mathbf{C}\_{(L+1)\*N} \\ \mathbf{C}\_{(L+1)\*(N-1)+2} \\ \vdots \\ \mathbf{C}\_{(L+1)\*(N-1)+p+1} \\ \vdots \\ \mathbf{C}\_{(L+1)\*N} \end{bmatrix} \tag{A20}$$

with *pǫ*[1, 2, . . . , *L*] and *jǫ*[1, 2, . . . , *N*].

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]:

*f* (*L*+1)∗(*j*−1)+<sup>1</sup> <sup>=</sup> h 1 + (1 − *ϕ*) 3*Demd*∆*t <sup>ρ</sup>pVwR*(*R*−*rL*) <sup>+</sup> 2*DL*∆*t* ∆*x* <sup>2</sup> + *<sup>v</sup>*∆*<sup>t</sup>* ∆*x* i*C k*+1 *w*,*j* − h (1 − *ϕ*) *DL*∆*t* ∆*x* <sup>2</sup> + *<sup>v</sup>*∆*<sup>t</sup>* ∆*x* i*C k*+1 *<sup>w</sup>*,*j*−<sup>1</sup> <sup>+</sup> h (1−*ϕ*)*DL*∆*<sup>t</sup>* ∆*x* 2 i *C k*+1 *w*,*j*+1 − h (1−*ϕ*)3*Demd*∆*<sup>t</sup> ρpVwR*(*R*−*rL*) i *C k*+1 *L*,*j* − h *ϕ DL*∆*t* ∆*x* <sup>2</sup> + *<sup>v</sup>*∆*<sup>t</sup>* ∆*x* i*C k w*,*j*−1 − h *ϕDL*∆*t* ∆*x* 2 i *C k w*,*j*+1 − h *ϕ*3*Demd*∆*t ρpVwR*(*R*−*rL*) i *C k L*,*j* − h 1 − *ϕ* 3*Demd*∆*t <sup>ρ</sup>pVwR*(*R*−*rL*) <sup>+</sup> 2*DL*∆*t* ∆*x* <sup>2</sup> + *<sup>v</sup>*∆*<sup>t</sup>* ∆*x* i*C k w*,*j f* (*L*+1)∗(*j*−1)+<sup>2</sup> <sup>=</sup> *<sup>ε</sup>* <sup>+</sup> *<sup>ρ</sup>pKf r C k*+1 1,*j* 1 *<sup>n</sup>* −1 + *De*∆*t r* 2 1 1−*ϕ <sup>r</sup>*1+0.5−*r*1−0.5 *r* 2 1+0.5 *r*2−*r*<sup>1</sup> *C k*+1 1,*j* − *De*∆*t r* 2 1 1−*ϕ <sup>r</sup>*1+0.5−*r*1−0.5 *r* 2 1+0.5 *r*2−*r*<sup>1</sup> *C k*+1 2,*j* − *De*∆*t r* 2 1 *ϕ <sup>r</sup>*1+0.5−*r*1−0.5 *r* 2 1+0.5 *r*2−*r*<sup>1</sup> *C k* 2,*j* − *<sup>ε</sup>* <sup>+</sup> *<sup>ρ</sup>pKf r C k* 1,*j* 1 *<sup>n</sup>* −1 − *De*∆*t r* 2 1 *ϕ <sup>r</sup>*1+0.5−*r*1−0.5 *r* 2 1+0.5 *r*2−*r*<sup>1</sup> *C k* 1,*j f* (*L*+1)∗(*j*−1)+*p*+<sup>1</sup> <sup>=</sup> *<sup>ε</sup>* <sup>+</sup> *<sup>ρ</sup>pKf r C k*+1 *p*,*j* 1 *<sup>n</sup>* −1 + *De*∆*t r* 2 *p* 1−*ϕ <sup>r</sup>p*+0.5−*rp*−0.5 *r* 2 *p*+0.5 *rp*+1−*r<sup>p</sup>* + *r* 2 *p*−0.5 *<sup>r</sup>p*−*rp*−<sup>1</sup> *<sup>C</sup> k*+1 *p*,*j* − *De*∆*t r* 2 *p* 1−*ϕ <sup>r</sup>p*+0.5−*rp*−0.5 *r* 2 *p*+0.5 *rp*+1−*r<sup>p</sup> C k*+1 *p*+1,*j* − *De*∆*t r* 2 *p* 1−*ϕ <sup>r</sup>p*+0.5−*rp*−0.5 *r* 2 *p*−0.5 *<sup>r</sup>p*−*rp*−<sup>1</sup> *C k*+1 *p*−1,*j* − *De*∆*t r* 2 *p ϕ <sup>r</sup>p*+0.5−*rp*−0.5 *r* 2 *p*+0.5 *rp*+1−*r<sup>p</sup> C k p*+1,*j* − *De*∆*t r* 2 *p ϕ <sup>r</sup>p*+0.5−*rp*−0.5 *r* 2 *p*−0.5 *<sup>r</sup>p*−*rp*−<sup>1</sup> *C k p*−1,*j* − *<sup>ε</sup>* <sup>+</sup> *<sup>ρ</sup>pKf r C k p*,*j* 1 *<sup>n</sup>* −1 − *De*∆*t r* 2 *p ϕ <sup>r</sup>p*+0.5−*rp*−0.5 *r* 2 *p*+0.5 *rp*+1−*r<sup>p</sup>* + *r* 2 *p*−0.5 *<sup>r</sup>p*−*rp*−<sup>1</sup> *<sup>C</sup> k p*,*j f* (*L*+1)∗*<sup>j</sup>* <sup>=</sup> *<sup>ε</sup>* <sup>+</sup> *<sup>ρ</sup>pKf r C k*+1 *L*,*j* 1 *<sup>n</sup>* −1 + *De*∆*t r* 2 *L* 1−*ϕ <sup>R</sup>*−*rL*−0.5 *R* 2 *R*−*r<sup>L</sup>* + *r* 2 *L*−0.5 *<sup>r</sup>L*−*rL*−<sup>1</sup> *<sup>C</sup> k*+1 *L*,*j* − *De*∆*t r* 2 *L* 1−*ϕ <sup>R</sup>*−*rL*−0.5 *R* 2 *R*−*r<sup>L</sup> C k*+1 *w*,*j* − *De*∆*t r* 2 *L* 1−*ϕ <sup>R</sup>*−*rL*−0.5 *r* 2 *L*−0.5 *<sup>r</sup>L*−*rL*−<sup>1</sup> *C k*+1 *L*−1,*j* − *De*∆*t r* 2 *L ϕ <sup>R</sup>*−*rL*−0.5 *R* 2 *R*−*r<sup>L</sup> C k w*,*j* − *De*∆*t r* 2 *L ϕ <sup>R</sup>*−*rL*−0.5 *r* 2 *L*−0.5 *<sup>r</sup>L*−*rL*−<sup>1</sup> *C k L*−1,*j* − *<sup>ε</sup>* <sup>+</sup> *<sup>ρ</sup>pKf r C k L*,*j* 1 *<sup>n</sup>* −1 − *De*∆*t r* 2 *L ϕ <sup>R</sup>*−*rL*−0.5 *R* 2 *R*−*r<sup>L</sup>* + *r* 2 *L*−0.5 *<sup>r</sup>L*−*rL*−<sup>1</sup> *<sup>C</sup> k L*,*j* (A21)

The residual function vector can be organized as:

$$f\left(\mathsf{C}^{k+1}\right) = \begin{bmatrix} f\_1 \\ f\_2 \\ \vdots \\ f\_{p+1} \\ \vdots \\ f\_{L+1} \\ \vdots \\ f\_{L+1} \\ \vdots \\ f\_{(L+1)\*(j-1)+1} \\ \vdots \\ \vdots \\ (L+1)\*(j-1)+p+1 \\ \vdots \\ (L+1)\*j \\ \vdots \\ (L+1)\*(N-1)+1 \\ \vdots \\ (L+1)\*(N-1)+2 \\ \vdots \\ \vdots \\ (L+1)\*(N-1)+p+1 \\ \vdots \\ \vdots \\ \vdots \\ \dots \\ \end{bmatrix}\_{(L+1)\times N} \begin{bmatrix} \mathsf{R} \\ \vdots \\ \mathsf{R} \end{bmatrix} \tag{A22}$$

Similar to the film diffusion case (see Appendix B), the *C <sup>k</sup>*+<sup>1</sup> vector can be determined by Equation (A9) as well and the Jacobian matrix of intraparticle pore diffusion case can be expressed as:

$$J = \begin{bmatrix} \frac{\partial f\_1}{\partial \mathbf{C}\_1} & \frac{\partial f\_1}{\partial \mathbf{C}\_2} & \cdots & \frac{\partial f\_1}{\partial \mathbf{C}\_{(L+1)\star N-1}} & \frac{\partial f\_1}{\partial \mathbf{C}\_{(L+1)\star N}}\\ \frac{\partial f\_2}{\partial \mathbf{C}\_1} & \frac{\partial f\_2}{\partial \mathbf{C}\_2} & \cdots & \frac{\partial f\_2}{\partial \mathbf{C}\_{(L+1)\star N-1}} & \frac{\partial f\_2}{\partial \mathbf{C}\_{(L+1)\star N}}\\ \vdots & \vdots & \cdots & \vdots & \vdots\\ \frac{\partial f\_{(L+1)\star N-1}}{\partial \mathbf{C}\_1} & \frac{\partial f\_{(L+1)\star N-1}}{\partial \mathbf{C}\_2} & \cdots & \frac{\partial f\_{(L+1)\star N-1}}{\partial \mathbf{C}\_{(L+1)\star N-1}} & \frac{\partial f\_{(L+1)\star N-1}}{\partial \mathbf{C}\_{(L+1)\star N}}\\ \frac{\partial f\_{(L+1)\star N}}{\partial \mathbf{C}\_1} & \frac{\partial f\_{(L+1)\star N}}{\partial \mathbf{C}\_2} & \cdots & \frac{\partial f\_{(L+1)\star N-1}}{\partial \mathbf{C}\_{(L+1)\star N-1}} & \frac{\partial f\_{(L+1)\star N}}{\partial \mathbf{C}\_{(L+1)\star N}} \end{bmatrix}\_{(L+1)\star N\times (L+1)\star N}$$

∆*C <sup>k</sup>*+<sup>1</sup> = *C k*+1 1,*guess*,*new* − *C k*+1 1,*guess*,*old C k*+1 2,*guess*,*new* − *C k*+1 2,*guess*,*old* . . . *C k*+1 *<sup>p</sup>*+1,*guess*,*new* − *C k*+1 *p*+1,*guess*,*old* . . . *C k*+1 *<sup>L</sup>*+1,*guess*,*new* − *C k*+1 *L*+1,*guess*,*old* . . . *C k*+1 (*L*+1)∗(*j*−1)+1,*guess*,*new* − *C k*+1 (*L*+1)∗(*j*−1)+1,*guess*,*old C k*+1 (*L*+1)∗(*j*−1)+2,*guess*,*new* − *C k*+1 (*L*+1)∗(*j*−1)+2,*guess*,*old* . . . *C k*+1 (*L*+1)∗(*j*−1)+*p*+1,*guess*,*new* − *C k*+1 (*L*+1)∗(*j*−1)+*p*+1,*guess*,*old* . . . *C k*+1 (*L*+1)∗*j*,*guess*,*new* − *C k*+1 (*L*+1)∗*j*,*guess*,*old* . . . *C k*+1 (*L*+1)∗(*N*−1)+1,*guess*,*new* − *C k*+1 (*L*+1)∗(*N*−1)+1,*guess*,*old C k*+1 (*L*+1)∗(*N*−1)+2,*guess*,*new* − *C k*+1 (*L*+1)∗(*N*−1)+2,*guess*,*old* . . . *C k*+1 (*L*+1)∗(*N*−1)+*p*+1,*guess*,*new* − *C k*+1 (*L*+1)∗(*N*−1)+*p*+1,*guess*,*old* . . . *C k*+1 (*L*+1)∗*N*,*guess*,*new* − *C k*+1 (*L*+1)∗*N*,*guess*,*old* (*L*+1)∗*N*×1 (A24)

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).

#### **Appendix D. Length of the Mass Transfer Zone (***Xs***) for the First Order Analytical Solution**

Analytical solutions can be derived for the case of the first flooding of the column which are used here for verification of the numerical codes.

#### *Appendix D.1. Analytical Solution Based on the Film Diffusion Model*

During the first flooding of the column, the front water flow is always contacting fresh contaminant material. Therefore, the solute concentration at the particle-water boundary is constant and in equilibrium with the solids:

$$\mathcal{C}\_{w,eq} = \frac{\mathcal{C}\_{s,ini}}{\mathcal{K}\_d} \tag{A25}$$

Inserting Equation (A25) into Equation (7) gives:

$$\frac{\partial \mathbb{C}\_w}{\partial t} = k \frac{\theta (1 - n)}{nd} \left( \mathbb{C}\_{w, \epsilon q} - \mathbb{C}\_w \right) \tag{A26}$$

which upon integration yields the following analytical solution for the initial condition *Cw*(*t*=0) = 0 (desorption):

$$\begin{split} \int\_{0}^{\mathsf{C}\_{\mathsf{w}}} \frac{\partial \mathsf{C}\_{\mathsf{w}}}{\left(\mathsf{C}\_{\mathsf{w},\mathsf{q}} - \mathsf{C}\_{\mathsf{w}}\right)} &= \int\_{0}^{t} k \frac{6(1-n)}{nd} \partial t \\ -\ln\left(\mathsf{C}\_{\mathsf{w},\mathsf{q}} - \mathsf{C}\_{\mathsf{w}}\right) + \ln\left(\mathsf{C}\_{\mathsf{w},\mathsf{q}}\right) &= -\ln\left(1 - \frac{\mathsf{C}\_{\mathsf{w}}}{\mathsf{C}\_{\mathsf{w},\mathsf{q}}}\right) = k \frac{6(1-n)}{nd} t \\ \frac{\mathsf{C}\_{\mathsf{w}}}{\mathsf{C}\_{\mathsf{w},\mathsf{q}}} &= 1 - \exp\left(-k \frac{6(1-n)}{nd} t\right) \end{split} \tag{A27}$$

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.

$$X\_{\text{s,63.2\%}} = \frac{v \ n \ d}{6k(1-n)}\tag{A28}$$

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 (*Xcol*), 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:

$$\begin{array}{rcl} t\_{63.2\%} &=& \frac{X\_{\varepsilon}}{\upsilon} + \frac{(X\_{col} - X\_{\varepsilon})}{\upsilon} R\_d\\ &=& \frac{X\_{col}}{\upsilon} \left( 1 + K\_d \frac{\rho\_b}{n} \left( 1 - \frac{X\_{\varepsilon}}{X\_{col}} \right) \right) \end{array} \tag{A29}$$

Considering fast kinetics ( *X<sup>s</sup>* → 0), *t*63.2%(≈ *Xcol*/(*v*/*R<sup>d</sup>* )) is mainly dominated by the retarded seepage velocity (*v*/*R<sup>d</sup>* ).

#### *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* = *De*/*δp*, where *D<sup>e</sup>* is the effective intraparticle diffusion coefficient (*D<sup>e</sup>* = *Daqε*/*τ* ≈ *Daqε 2* ) and the mean square displacement *δ<sup>p</sup>* (*δ<sup>p</sup>* = √ *πD<sup>a</sup> tc*) representing the diffusion distance, which grows with the square root of contact time between particles and water (*tc*) at early times, leads to:

$$\frac{\partial \mathcal{C}\_w}{\partial t} = k \, A^o \left( \mathcal{C}\_{w, \epsilon q} - \mathcal{C}\_w \right) = \frac{D\_\varepsilon}{\sqrt{\pi D\_a \, t\_\circ}} \frac{6(1 - n)}{nd} \left( \mathcal{C}\_{w, \epsilon q} - \mathcal{C}\_w \right) \tag{A30}$$

The contact time between water and dry particles can be estimated by the ratio of particle size and flow velocity (*t<sup>c</sup>* = *d*/*v*).

For the initial condition *Cw*(*t*=0) = 0, integration of Equation (A30) yields the following analytical solution:

$$\begin{split} \frac{\mathsf{C}\_{w}}{\int\_{0}^{\mathsf{C}\_{w}} \left( \mathsf{C}\_{w,\eta} - \mathsf{C}\_{w} \right)} &= \frac{\mathsf{f}}{\mathsf{0}} \frac{D\_{\mathsf{c}}}{\sqrt{\pi D\_{a}}} \frac{6(1-n)}{nd} \partial t \\ -\ln \left( \mathsf{C}\_{w,\varkappa q} - \mathsf{C}\_{w} \right) + \ln \left( \mathsf{C}\_{w,\varkappa q} \right) &= -\ln \left( 1 - \frac{\mathsf{C}\_{w}}{\mathsf{C}\_{w,\varkappa q}} \right) = \frac{D\_{\mathsf{c}}}{\sqrt{\pi D\_{a} \frac{d}{v}}} \frac{6(1-n)}{nd} t \\ \mathsf{C}\_{w,\varkappa q}^{\mathsf{C}\_{w}} &= 1 - \exp \left( -\frac{D\_{\mathsf{c}}}{\sqrt{\pi D\_{a} \frac{d}{v}}} \frac{6(1-n)}{nd} t \right) \end{split} \tag{A31}$$

ௗ

, μ

The length of the mass transfer zone of intraparticle pore diffusion can be calculated by:

$$\begin{array}{rcl} X\_{\text{s,63.2\%}} &=& \frac{v \text{ n } d}{6 \frac{D\_{\text{f}}}{\sqrt{\pi D\_{\text{u}} \text{f}}} (1 - n)}\\ &=& \sqrt{\frac{\pi d^3 v}{D\_{\text{f}} \left(\varepsilon + K\_d \rho\_p\right)}} \frac{n}{6(1 - n)} \end{array} \tag{A32}$$

*Xs*,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 (√ *Kd* ). 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).

#### *Appendix D.3. Comparison of Analytical and Numerical Solution and Estimation of Mass Transfer Zone Length (Xs)*

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<sup>d</sup>* , and much shorter than for IPD with *X<sup>s</sup>* = 10 cm, 3.5 cm and 1.1 cm for *K<sup>d</sup>* values of 0.1 L kg−<sup>1</sup> , 1 L kg−<sup>1</sup> and 10 L kg−<sup>1</sup> , respectively. The deviations between FD and IPD gradually vanish with increasing *K<sup>d</sup>* values. If the initial concentration in the column leachate is close to equilibrium, it may be used for the determination of *K<sup>d</sup>* (*K<sup>d</sup>* = *Cs*,*ini*/*Cw*,*peak*); *K<sup>d</sup>* is overestimated if the initial effluent concentration does not reach equilibrium (*Cw*,*peak* < *Cw*,*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. ௗ <sup>௦</sup> ௗ − − − ௗ ௗ ௗ = ௦,/௪, ௪, < ௪,)

௪, <sup>−</sup> <sup>−</sup> ௦, μ <sup>−</sup> <sup>−</sup> <sup>−</sup> ,௦ μ **Figure A2.** Concentration increase in a water parcel (*Cw*,*peak*) in the column during the first flooding (up-flow); solid lines: film diffusion; dashed lines: intraparticle diffusion; comparison between analytical (ana.) and numerical (num.) solutions. *<sup>n</sup>* = 0.45, *<sup>v</sup>* = 1.67 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m s−<sup>1</sup> , *α* /*x* = 0 (no dispersion), *Cs*,*ini* = 1000 µg kg−<sup>1</sup> , *<sup>t</sup>* = 5 h, *<sup>D</sup>aq* = 1 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup><sup>2</sup> s −1 , *ε* = 0.05, *dp*,*coarse* = 2000 µm.

ௗ

#### **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 (*dp*, *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<sup>d</sup>* values, are due to kinetics in the numerical solution. Deviations gradually vanish with the increase of *K<sup>d</sup>* .

<sup>−</sup> <sup>−</sup> ௦, μ <sup>−</sup> , μ **Figure A3.** Concentration vs. distance in the up-flow column test after the first flooding of the column (initial condition). Comparison of the analytical solution (Equation (6), dashed lines) and numerical solution (solid lines); *n* = 0.45, *<sup>v</sup>* = 1.67 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m s−<sup>1</sup> , *α* = 0 (no dispersion), *Cs*,*ini* = 1000 µ g kg−<sup>1</sup> , *t<sup>c</sup>* = 5 h, *dp*, *f ine* = 63 µ m. <sup>−</sup> <sup>−</sup> ௦, μ <sup>−</sup> , μ

௪/௪,, <sup>௪</sup> ௨ ௪/௪,, <sup>௪</sup> ௨ **Figure A4.** Normalized and absolute concentration (*Cw*/*Cw*,*eq*, *Cw*) as well as cumulative concentration (*mcum* ) in the column effluent vs. time (expressed as liquid to solid ratio: *LS*) for initial conditions shown in Figure A3; comparison of analytical solution (Equation (6), dashed lines) and numerical solution (solid lines).

#### **References**

