Development and Validation of a Box and Flux Model to Describe Major, Trace and Potentially Toxic Elements (PTEs) in Scottish Soils

The box and flux model is a mathematical tool used to describe and forecast the major and trace elements perturbations of the Earth biogeochemical cycles. This mathematical tool describes the biogeochemical cycles, using kinetics of first, second and even third order. The theory and history of the box and flux modeling are shortly revised and discussed within the framework of Jim Lovelok’s Gaia theory. The objectives of the investigation were to evaluate the natural versus anthropic load of Potentially Toxic Elements (PTEs) of the Scottish soils, investigate the soil components adsorbing and retaining the PTEs in non-mobile species, evaluate the aging factor of the anthropic PTEs and develop a model which describes the leaching of PTEs in layered soils. In the Scottish land, the soil-to-rock enrichment factorinversely correlates with the boiling point of the PTEs. The same is observed in NW Italy and USA soils, suggesting the common source of the PTEs. The residence time in soils of the measured PTEs linearly correlates with the Soil Organic Matter (SOM). The element property which mostly explains the adsorption capacity for PTEs’ is the ionic potential (IP). The downward migration rates of the PTEs inversely correlate with SOM, and in Scottish soil, they range from 0.5 to 2.0 cm·year−1. Organic Bentoniteis the most important soil phase adsorbing cation bivalent PTEs. The self-remediation time of the polluted soil examined ranged from 50 to 100 years. The aging factor, the adsorption of PTEs’ into non-mobile species, and occlusion into the soil mineral lattice was not effective. The box and flux model developed, tested and validatedhere does not describe the leaching of PTEs following the typical Gaussian shape distribution of the physical diffusion models. Indeed, the mathematical model proposed is sensitive to the inhomogeneity of the layered soils.


Introduction
In the geodynamic earth system, elements cycle between several compartments, most relevant of which are the atmosphere, hydrosphere, pedosphere, biosphere and lithosphere, but also across ecosystems in biogeochemical equilibria with these principal planetary reservoirs. As demonstrated by Jim Lovelock's Gaia theory, many of the important driving forces of the earth's dynamic biogeochemical system, as well as most of the relevant geochemical reaction kinetics, rule the flux of elements between the different earth reservoirs, modulated by the biological activity of living organism [1,2]. Most of the chemical equilibria controlling geodynamics are modulated by simple and multiple retroactive feedback, which ensures some buffering capacity to earth ecosystem and some degree of homeostasis to overall planet earth [3]. Int. J. Environ. Res. Public Health 2021, 18,8930 3 of 17 C, N, S and P biogeochemical cycles, but also to evaluate the short-and long-term effects caused by industrial growth that are kinetically perturbing the biogeochemical cycles of the major elements and trace PTEs present in the environmental ecosystem.
In this study, we used local soil systems to determine the soil-rock enrichment factor of PTEs, the anthropogenic PTE load and the properties of the soil that affect the mobility of that pollution. Furthermore, we evaluate the chemical-physical processes and kinetics that remove low-mobile anthropogenic PTEs adsorbed on soil colloids, in insoluble amorphous precipitates and in non-bioavailable phases occluded in the soil, forming crystalline mineral lattices. The ultimate aim of these measures was to develop and validate a new box and flow model that describes the mobility of PTEs in soils and, finally, to determine the so-called "self-remediation" or flushing rates of anthropogenic soil PTEs.

Description of Study Area
The soil selected to test and validate the model was developed on the top of a drumlin, several hundred meters long, located in Clelland, South Glasgow, Scotland, UK. The quaternary geology map of the British Geological Survey Service (www.largeimages.bgs.ac. uk/iip/mapsportal.html?id=1003899 (accessed on 21 April 2021)) indicated the soil being studied was derived from a fluvial-glacial parent material that had been deposited during the last glaciations. The alluvial soil was classified as Gleyic Fluvisol [21]. It has been polluted by atmospheric emissions from an adjacent smelter operation, peaking in the 1950s until 1987 when it was decommissioned. Climatic conditions of the area, average of the monthly data, was collected by the Scottish meteorological service (www.metoffice. gov.uk (accessed on 11 June 2000, stored by Gallini L.)), from 1990 to 1995, and are shown in Figure 1.
Therefore, box and flux models are mathematical tools that are widely used to not only model the earth biogeochemical equilibria, but also the response to perturbation of these equilibria. The box and flux model have been applied to describe the major elements C, N, S and P biogeochemical cycles, but also to evaluate the short-and long-term effects caused by industrial growth that are kinetically perturbing the biogeochemical cycles of the major elements and trace PTEs present in the environmental ecosystem.
In this study, we used local soil systems to determine the soil-rock enrichment factor of PTEs, the anthropogenic PTE load and the properties of the soil that affect the mobility of that pollution. Furthermore, we evaluate the chemical-physical processes and kinetics that remove low-mobile anthropogenic PTEs adsorbed on soil colloids, in insoluble amorphous precipitates and in non-bioavailable phases occluded in the soil, forming crystalline mineral lattices. The ultimate aim of these measures was to develop and validate a new box and flow model that describes the mobility of PTEs in soils and, finally, to determine the so-called "self-remediation" or flushing rates of anthropogenic soil PTEs.

Description of Study Area
The soil selected to test and validate the model was developed on the top of a drumlin, several hundred meters long, located in Clelland, South Glasgow, Scotland, UK. The quaternary geology map of the British Geological Survey Service (www.largeimages.bgs.ac.uk/iip/mapsportal.html?id = 1003899 (accessed on 21 April 2021)) indicated the soil being studied was derived from a fluvial-glacial parent material that had been deposited during the last glaciations. The alluvial soil was classified as Gleyic Fluvisol [21]. It has been polluted by atmospheric emissions from an adjacent smelter operation, peaking in the 1950s until 1987 when it was decommissioned. Climatic conditions of the area, average of the monthly data, was collected by the Scottish meteorological service (www.metoffice.gov.uk (accessed on 11 June 2000, stored by Gallini L.)), from 1990 to 1995, and are shown in Figure 1.

Soil Profile
Ap horizon (0-25 cm). This was a typically eluvial ploughed horizon, with a sharp boundary and a strong brown color due to the organic matter content. The porosity was high, and the plasticity moderate. The soil aggregates were irregular. Throughout the horizon, strong biological activity was observed: worms and abundant grassroots were found. A relevant flux of PTEs from soil to land biosphere could be expected from Bt horizon due to the intense plants root activity and plants uptake.

Soil Profile
Ap horizon (0-25 cm). This was a typically eluvial ploughed horizon, with a sharp boundary and a strong brown color due to the organic matter content. The porosity was high, and the plasticity moderate. The soil aggregates were irregular. Throughout the horizon, strong biological activity was observed: worms and abundant grassroots were found. A relevant flux of PTEs from soil to land biosphere could be expected from Bt horizon due to the intense plants root activity and plants uptake.
Ae horizon (25-35 cm): This was an eluvial horizon below the plough layer. It had a more diffuse boundary with respect to the Ap horizon and a yellow-brown color, which indicates the presence of Goethite and may be the rare mineral forming Maghemite, recognized as an indicator of biomineralization. These iron oxides are shown in the literature to be potential biominerals whose crystallization process is mediated by soil bacterial species. The soil-layer colors and presence of the yellow Goethite indicate a poor drainage state and suggest a seasonal redox oscillation within the Ae horizon, which is confirmed by the estimated soil porosity and the hydrologic balance between precipitation and evapotranspiration (EPT). The Ae horizon had low plasticity and moderate porosity. The soil aggregates were polyhedral with irregular fractures. The cracks observed in the Ae horizon suggest periodic changes in water saturation of the soil and the presence in it of swelling clay minerals in some uncertain amount but below 5%. The biological activity was low and reasonably was mainly represented by soil microflora and soil microfauna.
Bt horizon (35-40 cm). It was an illuvial horizon containing yellowish and greenish lenses that suggested patches of strongly reducing conditions. The red color indicated the presence of iron oxides, such as Hematite, and may be inherited Ilmenite. Porosity and plasticity of the horizon were very high. Soil aggregates were polyhedral with regular fractures, a feature that is proving that swelling clay minerals are present in the Bt horizon in more amount than the upper one Ae. Thin clay films were observed in some polyhedral cracks between aggregates, suggesting that translocation of clay minerals from the upper soil horizon is effective, and they accumulate in the Bt horizon in some quantity. In the Bt horizon, we observed cm diameter stones containing strongly weathered lumps ofcoal and mixtures of basaltic and sedimentary rock.

Soil Sampling
Superficial debris and vegetative materials were removed before the collection of samples from the soil surface. The soil was sampled by collecting slices of 5 cm thick at 5 cm increments from a surface of 20 × 20 cm 2 (about 3 kg dry weight for each layer). Soils were collected in plastic bags, homogenized and air-dried to constant weight, at room temperature; they were then softly disaggregated and sieved to remove coarse debris and rubble (>2.0 mm).

Analytical Methods
Physical-chemical parameters of soil were determinate in accord with standard methods [22]. The exchangeable pool was determined by extraction of 2 g soil samples (dry weight equivalent) suspended in 20 mL of BaCl 2 (10% w/v at pH = 5.5), so that the precipitation of metal oxides during extraction and the dissolution of any mineral phases were avoided, and were mechanically shaken for 30 min and sonicated for 15 min. The extraction procedure was applied five times until the residue of the extractable pool was negligible. The supernatant was collected by centrifugation and analyzed by ICP-AES or HGA-AAS.
Concentrations of "pseudo-total" metals were determined after digestion of soil samples with aqua-regia, which does not completely destroy silicates [21,23]. Then 2.0 g of dry soil was placed in a digestion vessel with 6 mL of HCl and 14 mL of HNO 3 and digested in microwave digestion unit (MLS-1200 Mega, Milestone Inc., Shelton, CT, USA), at rising temperature steps, and later filtered. The digested samples were analyzed for Si, Al, Fe, Ca, Mg, K, Na, Ni, Co, Cu, Zn, Cr, Mn, P and S by inductively coupled plasma-atomic emission spectrometry (ICP-OES; model iCAP 6000, Thermo-Scientific, Cambridge, UK). A certified standard soil reference material was also analyzed in order to evaluate the method of digestion and the accuracy of analysis (GBW7404). The external standard solutions were derived from 1000 ppm stock metal solution; multi-element standards for calibration were prepared in the same reagents to minimize interferences and matrix effects. Three replicate analyses for samples were performed for all determinations. The Standard Deviation (SD) was below 10% for pH, texture, Cation Exchange Capacity (CEC), and below 1% for the aqua-regia extractable pool [16,23].
All the chemicals used in this study were analytical grade or better (Sigma-Aldrich, St. Louis, MO, USA) and used without further treatment. Reagent solutions were prepared with deionized water of resistivity not less than 18.2 MΩ×cm -1 . The main soil properties are reported in Table 1. The box and flux model here developed and tested isthe generalization of the model of Bonazzola et al. [12]. The main differences of the new box and flux model are the following: (a) In the top surface reservoirs, a constant atmospheric input I(t) can be considered. (b) In the model here proposed, each soil reservoir (S n ) is characterized by a specific kinetic constant (K n ) which characterizes any different soil layer-indeed, the averaged K kinetic constant fitted across the entire soil profile, such as the Bonazzola's box and flux model calculates. (c) The general solution valid for Sn reservoir is found and proposed-indeed, the only four reservoirs model solution proposed by Bonazzola et al. [12].
As in the Bonazzola's model, the concentration of the PTEs in the soil profile are subsided in reservoir "Sn", and the leaching from a reservoir to the under layer one is assumed to be a first order kinetics. The leaching of the PTEs can be so described by the system of n homogenous differential equations listed below.
The system is resolvable, and the general solutions are as follows: Three geochemical parameters can be derived from the kinetic constant (K n ) and calculated for each soil reservoir (S n ): (a) residence time of the pollutants in the reservoir T r ; (b) half-life time of the pollutant in the reservoir, t 1 2 ; and (c) average downward migration rate V n of the pollutant in the reservoir S n . The relative formulas are listed as follows: where d n is the thickness of the reservoir S n . The residence time (Tr n ) represents the ratio between the amount of the element in the reservoir and the flux of element throughout the reservoir and at the steady state of the system; the residence time represents the mean time of an element that is remaining in the reservoir.
The half-life time (t 1 2 n ) at the steady state of the system is the time required to remove half of the element pool in the reservoir n. The lower the half-life, the faster the turnover of the element in the reservoir.
The downward migration rate (V n ) represents the average downward migration rate in the reservoir (S n ), and it is related to the partition constant for the adsorbing phases following the equation: where V 0 is the downward migration rate in absence of the adsorbing phase, and Kd i and Q i are the partition coefficient and the pseudo-total amount of the adsorbing phases "i" respectively pooled in the reservoir S n . The downward migration rate in absence of adsorbing phases can be assumed to be as follows: where P is the rainfall rate, I c is the infiltration coefficient, ETP is the sum of EPT and θ is the soil porosity.
When the values I (t) , S (n,0) and S (n,t) are known the values of the kinetic constant Kn characteristic of the reservoir Sn can be numerically calculated. To calculate K n from the value of I (t) , S (n,0) and S (n,t) , the solver extension for Excel was used, and the validity of the solutions found were checked by a self-produced software.
Previous studies have discussed whether the box and flux models applied to forecast the leaching of chemicals in soils are mathematical rather than physical models. Some authors agree that the mathematical box and flux models are predominantly equivalent to the advective-diffusion model at an infinite scale, e.g., when the S (n,t) reservoirs are infinite. A truly physical advective-diffusive model applied to soils requires the determination of several soil properties, such as infiltration rates, soil porosity, pore distribution, soil tortuosity, exchange sites selectivity distribution, colloids surface charge and many other soil properties difficult to be measured, detected or estimated. In this study, a box and flux mathematical model was developed and applied instead of a truly advective-diffusive physical model, as it is more straight forward to develop and numerically solve, and it useful when applied to a layered substrate, such as soil profile [24].

Hypothesis Applied in the Box and Flux Model and Validation
The history of pollution inputs into the soil is known; the values I (t) and S (n,0) were estimated by the value S (n,t) of the trace elements. Since the beginning of the 1900s until 1987, the soil was ploughed, and therefore the atmospheric input I (t) mostly accumulated in the plough layer [16]. After 1987, the industrial plant was closed and subsequently the ploughing of the soil ceased and cropping for barley was changed to grassland. Therefore, the trace element atmospheric input I (t) reduced close to zero, and the elements in the ploughed layer was dominated by downward leaching. The input of the trace elements by agricultural fertilizer are assumed to be negligible with respect to the atmospheric input I (t) following Riffaldi and Levi-Minzi [25].
To calculate the kinetic constant (K n ), the values of I (t) , S (n,0) and S (n,t) must be identified and were estimated by the following assumptions: (a) From the beginning of the 1900 until 1987, the soil received an atmospheric input I (t) which has appreciably increased the trace elements concentration in soil; (b) The atmospheric input I (t) can be disregarded after the end of smelter activity; (c) Until 1987 the soil was ploughed, and therefore the I (t) input accumulates in the ploughed layer and the leaching of trace elements is minimal, below 5% of the measured anthropogenic PTEs load; (d) Almost all trace elements were held in the plough layer until the end of ploughing; (e) Trace elements starts to be leached throughout the layers of the soil profile after the end of ploughing; (f) The biological uptake and biological loop of soil-plant-soil transferring fluxes is minimal and can be disregarded at the present time; (g) The lithogenic contribution to the trace element pool extracted by aqua-regia can be disregarded, and the S (n,t) value of PTEs can be estimated by the aqua-regia extraction at the time of sampling; (h) Most of the trace elements pool was bound to the labile fraction and the ageing factor is reduced. Therefore, the value of S (n,t) can be estimated by the aqua-regia extraction at the time of sampling. Anthropogenic PTEs are bound or immobilized in the soil biomass, associated with humic and fulvic compounds, labile soil mineral fraction and biominerals. Anthropogenic PTEs estimated by aqua-regia is almost exhaustive, and only a partially overestimated due to partial dissolution of primary mica by aqua-regia extraction [23,24].

Statistical Analysis
SPSS statistical package (Window version 25, SPSS Inc., Chicago, IL, USA) was used for data analysis. The analysis of the experimental data was carried out by using oneway ANOVA, Duncan multiple comparison test and Pearson correlation analysis. All significance statements reported in this study are at the p < 0.05 level [26,27].

Soil Profile Properties
The very coarse fraction reflects the sedimentary (primary) features of the parent materials because they are unaffected by pedological processes [28]. The stones are scattered in the Ap horizon, caused by agricultural disturbance, and progressively decrease to zero towards the bottom layer. It is suspected that the parental material was originally sediment in shallow and slow-moving water in fluvial-glacial environment, affected by current activity. A portion of the clay in the Bt horizon is likely to be primary and due to the Würm post glacial warm period. The pedological horizons reflect the soil clay content, which is homogeneous in throughout the three detected horizons. Only a slight increment in clay content is observed in the top layer of the Ap horizon and in the bottom layer of the Ae horizon. These variations are common in soils with high vertical flushing by soil water and assisted by organic colloids from the surface layers. This is a typical condition in a boreal Atlantic strong-weathered soil environment due to the high rainfall rates. The Soil Organic Matter (SOM) content is highest in the Ap horizon and exponentially decreases sharply in the Ae horizon. There is an increase in the organic matter towards the top of the Ap horizon because of the root zone and degrading organic matter accumulating in the soil canopy. The CEC increases both toward the top and the bottom of the profile. Both organic and inorganic colloids appear to affect the CEC and strong multiple correlations are observed between CEC, SOM and clay content: CEC = (0.144 ± 0.031) × Clay + (0.829 ± 0.185) × SOM; n = 9, r 2 = 0.948; F = 9.1 × 10 −5 The clay and organic matter contribution to CEC can be estimated respectively as 14 and 80 Eq·kg −1 . The clay minerals have a low CEC and poorly crystalline illitic clay minerals are expected to dominatein the soil profile. The amount of exchangeable Na, K, Mg and Ca (Eq·kg −1 ) exceeds the CEC; therefore, the soil has an 100% base saturation and is oversaturated with respect to carbonates, condition which theoretically stabilize the formation of pedogenic swelling clay minerals.
The pH ranges from sub-acid to acid and is quite constant along the soil profile. A little increase of the pH is observed in the organic top layer, and a subtle decrement is observed in the clay minerals rich bottom layer of the Bt horizon. The observed pH range is above the pK both of the SOM carboxyl and the clay minerals basal surfaces but is below the pK of Mn, Fe, and Al oxides surfaces. The SOM and the clay minerals are expected to be mostly negatively charged and the oxides are positively charged in the pH range observed in the soil profile. The exchangeable H + was estimated by pH (KCl) and was constant through the Ap horizon, and increases exponentially below it, in agreement with the observations on the soil colloids, and CEC.
Reducing conditions are evident in the Bt horizon where gray-greenish lens is observed. The relative humidity (RH) can be considered as a broad indication of reducing conditions and correlate with the clay and SOM: RH% = (6.311 ± 2.943) + (0.277 ± 0.060) × Clay + (1.993 ± 0.331) × SOM; n = 9, r 2 = 0.884, F = 1.58 × 10 −3 A reducing environmental low oxidative potential is expected to prevail in the bottom of Bt horizon due the high quantity of clay, high relative humidity and low oxygen activity. The pH value is sub-acid for most of the Ap and Ae horizons, but it results in being acid in the Bt, where it is probably buffered by the dissolution and precipitation reactions of oxides.

Pseudo-Total Elemental Pool
The aqua-regia extracts all the PTEs bound to the soil humic and fulvic compounds, all the PTEs incorporate into soil microfauna and microflora pool, almost all the PTEs bind to soil minerals labile pool, most of the PTEs bound to superior vegetal residues and only a little fraction of the PTEs sinks into primary soil forming minerals. Most of the elements extracted by aqua-regia are therefore the anthropogenic soil PTEs concentration, here mainly due to the smelter plant emission in the atmosphere and its deposition on the top investigated soil layers. The concentration of pseudo-total element pool extracted in the investigated soil reservoirs is reported in Table 2. It is evident from Table 2 that the PTEs load throughout the examined soil profile is inhomogeneous among the different soil layers, and specific for any element. No element studied shows a Gaussian distribution across the profile, as usually observed in most of the soil analysis published in the literature and observed in the field [23]. It is hypothesized that the PTEs load in any soil layers reflects the adsorbed component in the soil, as well as specific chemical properties of individual PTEs.
When comparing the aqua-regia extracted pool of the trace elements with the averaged soil composition reported in the literature, it is apparent that most of the trace elements in the investigated soil have a concentration closed to the global soil average [29,30]. Only in the case of Cd is the soil averaged value exceeded by the aqua-regia' extract soil-to-crust enrichment factor, which is close to two. It should be noted that the trace metal concentration in the aqua-regia extractable pools and the soil average concentration reported in the literature are not directly comparable, because the global average soil composition quoted in the literature results from a total (HF) digestion and not by aquaregia extraction. The aqua-regia soil extract and does not allow us to estimate the PTEs occluded in most of the silica primary minerals, and it is believed to be almost everywhere strongly sensitive to anthropogenic soil PTEs load [23,25].
The correlation matrix among the elements of the aqua-regia pool is shown in Table 3, and it allowsus to broadly distinguish four groups of the elements with similar distribution in the soil profile. The groups are (a) Na; (b) Al, Fe, Mg, K and Si; (c) Ni, Co, Cu, Zn, Cd, Pb, Mn and Ca; and (d) Cr, P and S.  The concentration of Na has a complex distribution throughout the soil profile and its pseudo-total concentration does not correlate with any other elements. It is reasonable to suppose that the concentration of Na in the soil is affected by the sea aerosol input for the sample site from on-land marine transfer. The concentration of K, Fe, Mg, Al and Si was almost constant within the soil profile, and a concentration peak of these elements are observed in the Bt horizon. Elements K, Mg, Al, Fe and Si results strongly correlate to the clay fraction content, and therefore is reasonable believe that the soil clay minerals are the main source of the aqua-regia extracted pool of these soil chemical element. The group of the bivalent cation elements Ni, Co, Cu, Zn, Cd, Pb, Mn and Ca has complex patterns in the distribution of the concentration along soil profile. Correlations among the concentration of these elements in the soil profile are observed, and correlation varies from weak to strong. The concentrations of the elements of this group are correlated with SOM and for some elements even with the soil clay fraction, such as Ca, Co and Ni. The group of elements stable as anions (Cr, P and S) has a pseudo-total concentration almost monotonically decreasing across the soil profile. Positive strong correlations between pseudo-total pool and organic matter are observed, and it is reasonable to assume that these anionic species are strongly bounds to the SOM colloids.

Soil Enrichment Factor
The soil to crust enrichment factor is the ratio between the element concentration in the soil and the element concentration in the crust. The soil enrichment factor (SEF) in the investigated soil was calculated by assuming that the pseudo-total pool of the measured elements is due to the anthropicfluxes and to the perturbation of the natural biogeochemical cycles, dividing the aqua-regia extracted pools values for the average crust concentration of the PTEs. The SEF calculated for the soil is reported in Table 4. The enrichment factor for Ca, Mn, Cd, Co, Cr, Cu, Ni, Pb, Zn, P and S in the studied soil result correlates with the USA soil-to-crust enrichment factor.
If the SEF of the S is disregarded, a strong correlation is observed between the pseudototal SEF in Clelland and the soil enrichment in the USA: This result is in agreement with previous investigation on anthropic Zn, Pb and Cd in NWItaly [30]. The SEF of Cd, Co, Cr, Cu, Mn, Ni, Pb and Zn observed in the soil does not correlate with the ionic potential (IP) of the element; instead, it negatively correlates with the boiling temperature of the element (Figure 2). The correlation observed between SEF and element boiling temperature suggests that the SEF of the PTEs is due to the smelter plant activity rather than soil chemical pedological processes both in Clelland, NW Italy, the USA and in Canadian soils [13,15,16,18,31]. with the boiling temperature of the element (Figure 2). The correlation observed between SEF and element boiling temperature suggests that the SEF of the PTEs is due to the smelter plant activity rather than soil chemical pedological processes both in Clelland, NW Italy, the USA and in Canadian soils [13,15,16,18,31].

Residence Time of Element in the Reservoirs
The residence times of the pollutants in the reservoirs of the soil profile are shown in Table 5.
The soil properties which control the residence time have been investigated by forward multiple regression analysis. The main soil property which controls the residence time of Na is the clay content. The main soil property which controls the residence time of the bivalent cations, i.e., Ca, Co, Ni, Cu, Zn, Cd, Pb and Mn, is the SOM, while for the anions (Si, P and S), it is both the organic matter and clay fraction: Assuming that the SOM is the dominant soil property controllingresidence time (RTE) of the elements in the topsoil layer, the RTE can be described by the following simple linear model: The parameter α has been calculated by multiple regression for Ca, Co, Ni, Cu, Zn, Cd, Pb, Mn, Si, P and S; the results are shown in Table 6, and the slopes of the linear regression results α has been correlated withthe IP of the element. The results differ between divalent cations and anions. For the divalent cations, the constant α is inversely proportional to the IP of the metal, while, for the anions, it is directly proportional to the IP (Figures 3 and 4).   The regressions observed between parameter α and the IP of the elements are as follows:  The regressions observed between parameter α and the IP of the elements are as follows: The regressions observed between parameter α and the IP of the elements are as follows: Bivalent cations: α = (7.85 ± 1.95) − (2.50 ± 0.80) × IP − 1; n = 8, r 2 = 0.622, F = 2.00 × 10 −2 Anions: α = (1.88 ± 0.58) + 0.306 (±0.05) × IP; n = 3, r 2 = 0.973, F = 1.05 × 10 −1

The Downward Migration Rate of the Pollutants
The calculated downward migration rates are listed in Table 7. The downward migration rates increase in the following order:  Table 7. Downward migration rates of the pollutants, cm·year −1 (*** = not determinable). *** *** *** *** *** *** *** *** *** *** *** *** As a general rule in the profile, the downward migration rates have a minimum in the top layer, increase to a maximum in the Ae horizon and decrease again approaching the Bt horizon.
It has been assumed that the abundance of the adsorbing phases should be linearly correlated with the logarithm of the downward migration rates, V (E) , of the element. This has been tested by multiple regression analysis. In Table 8, the correlation matrix between ln(V E ) and the soil properties is shown. Most of the PTEs correlate both with SOM and CEC, but the correlation with clay content alone is poor.   (Table 9). Forward multiple regression analysis demonstrates that the bulk clay fraction is not meaningfully correlated to ln(V E ) when SOM and clay fraction are considered together. It is reasonable to assume that the downward migration rate of the considered PTEs is sensitive of the mineralogical composition to the clay fraction, as previously demonstrated for 137 Cs [15,31,32]. Indeed, the Gibbs free energy of the adsorption of cations on clay minerals is very different among the several mineralogical species [24,33].
It has been suggested that the above-reported correlation could be expressed as a linear function of the SOM following the model: Following the proposed linear model, ln(V 0 ) is the logarithm of the downward migration rate in absence of organic matter, and Kd SOM is the fraction of the mobile element retained in the layer because of the adsorption of the element into the organic matter. The values of ln(V 0 ) and Kd SOM have been estimated by linear regression above.
It is reasonable to assume both the parameters ln(V 0 ) and Kd SOM could be related to chemical properties of the elements.
With respect to anions, no meaningful correlation was found between element properties and ln(V 0 ) or Kd SOM . With respect to bivalent cations Co, Ni, Cu, Zn, Mn, Cd, Ca and Pb, the parameter ln(V 0 ) is linearly correlated to the ∆G* of the adsorption of the cation on bentonite for the divalent cations, as shown in the following equation: ln(V 0 ) = (0.689 ± 0.128) + (0.0123 ± 0.0047) × ∆G*; n = 6, r 2 = 0.634, F = 5.80 × 10 −2 The downward migration rate (V 0 ) of the element therefore decreases exponentially as the adsorption constant on bentonite increases. The parameter Kd SOM for the bivalent cations Co, Ni, Cu, Zn, Mn, Cd, Ca and Pb is resulted inversely proportional to the ionic radius of the element following the linear model: Kd SOM = (0.026 ± 0.049) + (0.471 ± 0.121) × IP − 1; n = 8, r 2 = 0.671, F = 7.93 × 10 −3 Assuming that the soil properties that control the downward migration rates of the elements do not change with time and are in the steady state, the proposed box and flux model allowsus to predict the evolution with time of the element concentration throughout the soil profile. This assumption of the steady state of the soil properties controlling the PTEs leaching is reasonable when the concentration of the SOM pool in the soil profile has reached the steady state. It is not true if the climate change impact on land biomass and plant productivity are considered.
It is evident from the applied model that the time needed to decreases at half the concentration of the pollutants in the ploughed layer ranges from 25 years for the slowest mobile element (S) to 10 years for the faster moving (Na). The distribution of the pollutant throughout the soil profile calculated by the proposed box and flux model does not follow a Gaussian shape distribution, unlike in the case of the usual convective-advective physical models. It is because some soil reservoirs selectively accumulate some PTEs with respect to others PTEs. The reason for this phenomenon is that the proposed box and flux model is sensitive to the inhomogeneity commonly observed in the layered soil profile.

Conclusions
The main pedological processes that are active in the soil profile are ruled by SOM accumulation; weathering of the primary minerals, particularly the mica family; the formation of swelling clay minerals, which are stabilized by the very high soil base saturation; and the translocation and precipitation of soil colloids and oxide phases. The soil profile studied is characterized by almost constant pH values across the profile, which range from acid to slightly-acid. A sharp pH decrease is observed at the top of the Bt horizon, where it is reasonable to suppose it is buffered by iron oxide precipitation. As we observed soil structures and the related colors, wide Eh values are expected because of seasonal weather fluctuations and local soil structures at microscales. Among the elements studied K, Mg, Al and Si extracted by aqua-regia are predominantly derive from illitic minerals, while the aqua-regia extractable pool of Co, Ni, Cu, Zn, Cd, Pb, Mn and Cr is soil contamination derivedfrom the local smelter plant emissions tothe atmosphere and deposition on the topsoils. The soil-to-crust enrichment factor of the PTEs is explained by the element boiling point. As the boiling point decreases, the soil-to-rock enrichment factor increases, and with increasing boiling point, the SEF decreases. The same trends were observed in Northern Italy, the USA and in Canadian soils [13,15,16,18,31]. This demonstrates that extensive impact industrial emissions can have on food production. The average downward migration rates range from 0.89 to 1.75 cm x y −1 , in the following order: S (0.89) < Zn (0.92) = P (0.92) ≈ Pb (0.93) < Mn (1.06) < Cu (1.14) < Ca (1.21) ≈ Ni (1.23) ≈ Co (1.,24) ≈ Si (1.26) < Cd (1.42) < Na (1.75) The main soil component retarding the leaching of the monovalent Na is the clay fraction. The SOM is not effective in Na adsorption in a detectable way. The main soil component adsorbing the divalent cations and anionic PTEs is the SOM. The SOM and, to the lesser effect, the soil organic Bentonite phase are the main soil components affecting the PTEs mobility. The affinity of the PTEs with the SOM is inversely correlated to the IP. The residence times of the individual PTEs in soil are in the following order: Na (3.3) < Ca (4.8) < Co (5.5) ≈ Si (5.7) < Ni (6.3) ≈ Cu (6.2) < Zn (7.5) < P (7.8) < S (10.1) < Cd (11.7) < Pb (23.9) < Mn (28.9) The residence time of the Na is ruled by adsorption on clay minerals. The residence time in soils of divalent cations Ca, Co, Ni, Cu, Zn, Cd, Pb and Mn, and the anions Cr, S, Si and P, are ruled by adsorption on SOM. The main PTE property which explains the PTE affinity to the SOM is the IP. The correlations found among the measured PTEs and affinity to soil matter should allow us to estimate the affinity to SOM of PTEsnot-measured in our study, such as W, Se, As or the Lanthanides elements for which few data are reported in the literature. The biogeochemical kinetics of the trace PTEs obey first-order kinetics in the box and flux model. The box and flux models proposed and tested here allow us to forecast the downward migration rates and the theoretical time to self-clean polluted soil. The box and flux models describing the leaching in soils of PTEs do not follow the Gaussian distribution shape typical of the physical advective-diffuse models because it is sensitive to layered soil horizon in homogeneity. The result obtained in the Scottish soil agree with the data obtained investigating Alpine soils developed on crystalline basement [16,23,34]. The box and flux model has provided a useful geochemical tool to investigate, describe and forecast the fate of trace PTEs contaminating the environmental ecosystems, as well as the potential human health impact of soil contamination by smelter plant activity [13].