A consistent model for estimating the partitioning of Am, Pu and Se in agricultural soils

The component additive model UNiSeCs II for simulating the physicochemical behaviour of the radionuclides americium, plutonium and selenium in agricultural soils is presented. The model is validated by estimating the distribution coefficients (Kd) of these elements measured in batch experiments from the literature. For all three elements, the resulting average relative deviations from the experimental values are smaller than a factor of 2.5. This indicates that the model has the potential to significantly improve the predictions of radioecological models that normally use tabulated Kd values from the IAEA which are known to have large uncertainties. Using UNiSeCs II, the soil solution parameters most important for the partitioning of Am, Pu and Se are identified by single parameter variations.


Introduction
Contamination of agricultural soils by radionuclides may occur in the course of nuclear accidents as in Chernobyl or Fukushima where the main route is deposition onto plants and into the soil, being transported by rain or irrigation water. Another scenario is the leakage of long-lived radionuclides from a nuclear waste repository (NWR) and their subsequent migration into groundwater and surface soil. The soil can serve as a large reservoir for these nuclides, depending on their chemical properties and the chemical composition of the soil in question. Often a large proportion of the contaminant will be adsorbed to the surface of soil particles, but a certain amount will be present in the liquid phase ("soil solution"). This partitioning is described by the so-called distribution coefficient or K d (in L/kg) which is defined as the ratio of the activity sorbed to the soil solids (in Bq/kg) to the activity in soil solution (in Bq/L). Plant roots generally do not distinguish between contaminants and nutrients and thus, radionuclides may be taken up from the liquid phase and incorporated into various parts of the plant, depending on the chemical nature of the contaminant. If the plant is edible, the radioactive contaminants enter the food chain and are finally ingested by humans.
For the prediction of radiation exposure in a certain scenario an estimate for the K d of radionuclides in soils is required because this largely determines the amount that is available to the plants. Current radioecological models like Ecolego [1] use literature values of K d for exposure calculations. These tabulated values often vary by orders of magnitude even within one soil type [2] which is one of the significant causes for the uncertainty of such calculations. In this situation, geochemical modeling can be helpful to (i) estimate the K d for a given soil/radionuclide combination, (ii) identify the most relevant parameters that determine the partitioning of the nuclide which gives the opportunity for K d uncertainty prediction and (iii) contribute to the understanding of the key processes that lead to plant uptake. To perform the necessary geochemical equilibrium calculations, one can use geochemical codes like MINTEQA2 [3], Chess [4] or PHREEQC [5] that basically solve large mass balance equation systems for thermodynamical equilibrium in solution including the correspondent equation systems describing the interaction of ions with solid surfaces. If suitable thermodynamical data are provided, one can choose a subset of soil components like clay or iron hydroxides that are known to be active sorbents. Thus, in addition to a thermodynamical database for reactions of the relevant nuclides in solution, the inclusion of a database for (surface) complex 1 3 formation is essential. These assemblage models are also called component additive or CA models [6]. Examples for CA models are WHAM [7] which has been developed especially for situations where the chemical speciation is dominated by organic matter or the "generic multisurface sorption model" by Dijkstra et al. [8], which has been devised for the partitioning of heavy metals.
Hormann [9] developed a CA model called UNiSeCs that already had been successfully used for calculating the speciation and partitioning of uranium, nickel, selenium and cesium in agricultural soils [10]. This model has recently been updated [11] and has now been further developed and modified to include the radionuclides plutonium, americium and selenium. These elements all have long-lived isotopes relevant in the NWR scenario but are chemically very different. Pu and Am are present as cations, while Se mostly occurs in the shape of the negative ions selenite or selenate; Pu and Se are redox-sensitive while Am is not. The properties and the behaviour of these elements in the environment have been extensively reviewed in the literature [12][13][14][15][16] and will not be described in detail here.
The extended model is called UNiSeCs II; it uses the well-known geochemical code PHREEQC [5] for speciation calculation and the Thermochimie database [17] for the thermodynamical solution equilibrium data. 1 The model structure is shown in Fig. 1, further details (e.g. calculation of the amounts of binding sites) have already been described elsewhere [9,10]; the modifications and adjustments that have been made for UNiSeCs II are emphasized below.
The model's concept of radionuclide partitioning involves a rather simplified picture of soil composition (Fig. 1). The most important sorbents have been identified as (i) clay minerals, (ii) hydrous oxides of iron, aluminium and manganese and (iii) humic matter [8,18]. This concept is the same for all nuclides, meaning that all components are always present in the modeling assemblage.
Illite is a mineral that has been extensively studied as a sorbent for radionuclides [19,20] and serves as a proxy here because of its intermediate position within the clay minerals regarding number of binding sites [18]. In the calculation of illite binding sites, the clay mass is corrected by subtracting (i) the mass of hydrous oxides of Fe, Al and Mn and (ii) 50% of the mass of organic matter. This has already been justified in [10]. In UNiSeCs II, the thermodynamical data from Bradbury and Baeyens [19,20] are only used for the complexation of Am. For Pu, the thermodynamical constants from Banik et al. [21] are used because they are based on experimental data and not on linear free energy relationships. The complexation of Se (as selenite) on illite is additionally accounted for by the model of Missana et al. [22]. This is not coupled with the Bradbury and Baeyens model because it uses different protonation constants.
The hydrous oxides are being taken into account using the model of Dzombak and Morel [23] for Se. In UNiSeCs II, the model of Zavarin and Bruton [24] is used for Pu and Am because the complexation constants are based on experimental results rather than on the estimation by linear free energy relationships as described in [23] and [25]. This model is non-electrostatic, which implies that the influence of the surface potential is neglected. For both models it is assumed that the complexation constants and the amount of sorption sites for aluminium oxides are similar [26]. In UNiSeCs II, manganese oxides (which are usually present in small amounts) are also included and treated like the other hydrous oxides. In manganese-rich soils and in cases where a radionuclide is specifically bound to manganese oxides, one may have to use a specific model for this component.
Organic matter (OM) and dissolved organic matter (DOM) are the remaining essential soil components that have to be considered. In UNiSeCs II, Model VII by Tipping et al. [27] for the sorption of cations on humic matter replaces the previous Model VI and also the complexation constants from Stockdale et al. [28] for Pu are included. It is assumed that 80% of the immobile OM is humic substance [29,30] and can be attributed to humic acid whereas 38% of dissolved organic matter (DOM) is "active" with respect to cation binding and can be attributed to fulvic acid [31].
All submodels except those of Zavarin and Bruton [24] and Missana et al. [22] allow for the sorption of competing ions such as Ca 2+ and Fe 2+ /Fe 3+ . A synopsis of the submodels is given in Table 6 of the appendix.
For consistent modeling of partitioning, it is desirable that as many modeling assumptions as possible are general and to test these assumptions by validation with chemically different radionuclides.
In this work, UNiSeCs II is validated for each radionuclide by comparison of experimental results from batch experiments in the literature with predictions of model simulations.
To apply the model under field conditions for a certain soil, it is important to know its characteristics not only with respect to its solid composition, but also regarding the soil solution, namely pH, pe, DOC and ions that may act as complexants or as competitors for sorption. This has been evaluated for each radionuclide by varying single parameters within ranges frequently found in agricultural soils and comparing the resulting impacts on the distribution coefficient. Thus, the most important soil solution parameters have been identified for each radionuclide.

Methodology
For all simulations, two solutions and the surface assemblage are defined in PHREEQC: • The solution present in the soil before the addition of contaminant (SOLUTION 1). If not analyzed, its composition is approximated by a "standard" soil solution; for details, see Hormann and Fischer [10], Table 2. Fe and Al concentrations are estimated by assuming that the soil solution is in equilibrium with the minerals Goethite and Gibbsite, respectively. To account for the undersaturation of Al in the soil solution with respect to Gibbsite at low pH [32], the saturation index was set to -1 below pH 6.3 as a first approximation. These assumptions yielded the best results in the validation study. • The solution used to contaminate the soil in the experiment (SOLUTION 0). In the case of Se, it is pure water as this was an extraction experiment [33]. • The surface assemblage is defined via the EXCHANGE and SURFACE keywords. The calculation of surface sites has already been explained in [10].
After these definitions, the surface assemblage is brought into equilibrium with SOLUTION 1 (PHREEQC command -equilibrate 1). This is to estimate the initial state of the surface sites regarding occupancy with competing ions. In a final step, the preequilibrated surface assemblage is equilibrated with SOLUTION 0 using a single ADVECTION step.
For the sensitivity study, both solutions are approximated by the "standard" soil solution, with the addition of the contaminant in SOLUTION 0. This is to simulate the arrival of a dissolved "concentration front" in a previously uncontaminated soil volume.
In both sensitivity and validation studies, the K d values are finally calculated in a USER_PUNCH section using PHREEQC basic commands by summing up the radionuclide content of the different surfaces and dividing this by the final solution concentration.

Model input and assumptions for the validation study
For validation, experimental studies are needed with a sufficient amount of information on soil and soil solution composition that can be used as input parameters for modeling. Although there are lots of reports with experimental K d values in the literature, in most cases important data about the soil characteristics or the details of the experiments themselves are missing. The choice of studies from the literature was based on the requirement that (i) at least 10 different soils had been analysed, (ii) the soils were covering a large range of properties, e.g. with respect to pH or clay content and (iii) the soils had been analysed for clay, hydrous oxides, organic carbon and soil solution composition including DOC, preferably before and after batch extraction. However, even after an in-depth literature search no study was found that completely fulfilled criterion (iii).
Finally, three studies (for Am: Ramírez-Guinart et al. [34]; for Pu: Miner et al. [35]; for Se: Tolu et al. [33]) were found to be most suitable; each study covers a range of soil properties that is wide enough and provides enough details for modeling. However, in all cases additional assumptions have to be made because data are missing. To exclude any significant effects by sorption to calcite, only soils with a calcite content of ≤5% have been chosen, as there is currently no submodel for sorption to calcite in UNiSeCs II. Highly alkaline soils (pH > 8.5) have also been excluded. All studies had been performed in the form of batch experiments; in the studies for Am and Pu, the soils had been contaminated with a spiked solution, whereas in the Se study, the ambient selenium had been extracted by ultrapure water. The problem with this kind of batch experiments is that the state of the soil/solution system prior to contamination resp. extraction is usually unknown. This leaves an uncertainty concerning the coverage of the soil sorbents with competitive ions like Fe 3+ or PO 4 3− . As a best approximation, this is estimated by equilibrating the soil surface assemblage with the aforementioned "standard" soil solution ("pre-equilibration") with the additional assumption that Fe and Al are in equilibrium with Goethite and Gibbsite, respectively. For pH ≥ 7, Ca and Mg are assumed to be in equilibrium with calcite resp. magnesite if soil solution values are not given. This step has always been done prior to the simulation of the contamination of the soil assemblage.
In the case of Am, the soil solution had been analyzed for Ca, Mg, Na, K, DOC and dissolved inorganic carbon; for modeling, these values have been used both for the pre-equilibration solution and the contaminating solution. The Pu study contains no information about hydrous oxides; thus, they have been estimated using the average of four reference soils (definition of "reference soil" see below). However, as it turns out, this compartment can be neglected in the simulations. Miner et al. [35] specify the pe of the final solutions; in the simulations, these values have also been used for the respective pre-equilibration solutions. The Thermochimie database contains a solid phase "PuO 2 (coll,hyd)" for which the saturation index calculated by PHREEQC is always negative under the conditions assumed here (the presence of DOM strongly reduces the concentration of free Pu ions); thus, the presence of Pu eigencolloids is not accounted for. For Se and Am, the default value of PHREEQC has been applied (pe = 4). The redox sensitivity of Am is expected to be small, because under ambient conditions, it is exclusively present in the trivalent state; for Se, which is redox sensitive, such a low value seems to be justified, because ultrapure water (containing only traces of O 2 ) had been used for extraction.
In all cases, the pre-and final equilibrations have been performed using the proper solid-liquid ratios, assuming room temperature (20 °C). If the content of DOM is not reported in the experimental study, it has been estimated from pH, OM content and solid/solution ratio by the empirical relation found by Römkens et al. [36].
The K d values have been calculated using the sum of the contaminant sorbed to illite, hydrous oxides and humic acid, as well as the respective concentration in solution including the amount sorbed to fulvic acid.

Model input and assumptions for the sensitivity study
To identify the most significant solution parameters, K d values have been calculated varying single soil solution parameters in the range that is frequently found in agricultural soils according to Table 5.2 in [18]. The calculations have been performed using the soil composition data of a so-called Refesol (reference soil), taken from a set of soils being used as a reference for chemical soil testing in Germany [37]. The soil used here is Refesol 01A, which is a common soil type in northern Germany; its characteristics relevant for modeling are shown in Table 1. In the simulations presented here, the soil assemblage is assumed to be at field capacity and already in equilibrium with the soil solution when a trace amount (in the µg range, to exclude saturation effects) of contaminant is introduced. All calculations except for the Table 1 Relevant characteristics of Refesol 01A, as analyzed by the Fraunhofer Institute for Molecular Biology and Applied Ecology (IME); clay determined after DIN ISO 11,277, C org determined after DIN ISO 10,694. The subscript "ox" means "oxalate extractable".
Analytical uncertainties are ≤3% [38] Clay (%) Field capacity (g kg −1 ) pH (CaCl 2 ) C org (%) CEC eff (mmol kg −1 ) Al ox (g kg −1 ) Fe ox (g kg −1 ) Mn ox (g kg −1 ) pe resp. pH variations themselves have been performed at pe = 8 and pH = 5.6. First of all, the effect of DOC, pH (in the range 4-8.5) and pe (in the range 0-12) on K d has been tested. Additionally, the concentrations of the most abundant elements/ions in the soil solution that may interact with the surfaces present in the model (Ca, Mg, inorganic Fe, inorganic Al, phosphate, sulfate, carbonate and nitrate) have been varied, because they may compete with radionuclides for sorption or act as complexants in solution. However, it has to be mentioned that in the case of pH and pe, these are not strictly single parameter variations because in the model, Fe and Al concentrations are determined by solid phases and the solubility of those phases depends on acidity and redox status. The Gibbsite saturation index correction for Al mentioned above has been omitted for simplicity and is only shown in Fig. 5.
Moreover, in the case of Fe, the solution concentration is determined by DOM-bound Fe and also colloid Fe nanoparticles may occur in solution [39][40][41]. The omission of Fe colloids in solution is probably the reason why the calculated values for Fe in solution are underestimated by an order of magnitude [42]. However, as the calculated amount of free and DOM-bound Fe will not depend on whether the Fe mineral controlling these values is suspended in solution or not, the modeled distribution coefficients will be unaffected. Effects of radionuclide sorption to suspended mineral particles have not been included in this analysis because there are neither sufficient data on representative ranges of particle concentrations nor respective data on specific surface site concentrations (the particles may have a higher surface to volume ratio, but may be also coated by humic substances).
In principle, temperature can also be varied in PHREEQC, but as the the standard enthalpy changes for the complexation constants are not given in the literature, the temperature dependence of these constants cannot be calculated by the van't Hoff equation and therefore, the results would probably be misleading.

Validation
In this section, modeling results are compared to the experimental data from the three sources cited above. In the case of Am and Se, the experimental uncertainties are so small that they do not exceed the dimensions of the diamonds given in the respective figures. In the bar diagram for Pu, the error bars are much smaller than the differences between the measurements at different concentrations; thus, they have been left out for convenience.
Considering the uncertainty of all parameters used in this model (e.g. soil composition, number of binding sites, state of the system prior to contamination and the general assumptions that have been described above) it is clear that one cannot expect a precise prediction of the K d for any particular soil.

Americium
The results for americium are shown in Fig. 2 where the calculated K d is plotted versus the experimental values from [34]. In 8 of the 11 cases, the values calculated by the model deviate by a factor of less than 2.5 from the "ideal line" depicted in Fig. 2. This is remarkable considering the conceptual uncertainties of the model. Looking closer, this result may be partly due to the fact that the simulations show that generally the partitioning of Am is entirely governed by the interplay between DOC and immobile OC. Experiments showing a strong association between Am and organic matter have already been described in the literature, see e.g. [43][44][45]. In effect, this simplifies the model structure because (with the exception of very acid soils, see sensitivity section below) the contributions of hydrous oxides and illite can be neglected, thus removing the uncertainties associated with the assumptions for calculating the respective number of binding sites. A major reason for the remaining discrepancies is probably due to the complete identification of DOM with fulvic acid and OM with humic acid, as well as the assumptions for the percentage of "active" binding sites. It is not yet clear whether attempts to optimize the model in that respect would lead to significant improvements, because the composition of DOM and OM varies from soil to soil [46]. For the remaining three soils, the K d has been underestimated by a factor of 3.8, 5.4 and 9.9, respectively. The reason for these larger deviations cannot be identified, they may be due to the uncertainty concerning Al-and Fe content of the pristine soil solution, but may also be caused by an unknown soil component that acts as an additional sorbent for Am.

Plutonium
The comparison of experimental and simulated distribution coefficients for plutonium is presented in a bar diagram (Fig. 3a) because in this representation, it can be clearly seen that the model overestimates K d for pH ≈ 8, while at acid to neutral values, the simulations are in the correct order of magnitude. These overestimations are probably mainly caused by two reasons: as we will see later, the K d decreases with DOC in solution and increases with Fe in solution.
The formula of Römkens et al. [36] by which DOC is calculated has been derived from soils in the pH range of 1.8 to 7.9 with a median value of 4.8, meaning that it may not be valid at the upper end of the pH range, as well as for the four slightly alkaline soils from Fig. 3. If we assume a DOC content of 5.5 mg/L (median value from Table 5.3 in Blume et al. [18], divided by the dilution factor of 5 from the experimental study [35]) for all soils with pH > 6, the overestimations are strongly reduced. Moreover, it is well known that Fe concentrations in solution are very low already under slightly alkaline conditions [47,48]. If we tentatively take this into account by setting the saturation index of Goethite to -1 in the preequilibration solution, we get even more satisfactory results. The effect of these corrections is shown in the bottom of Fig. 3. The remaining K d overestimations have most probably been caused by non-equilibrium in the respective batch solution experiments. Only in these cases, the deviations from the arithmetic mean of the experimental values are larger than a factor of 5.

Selenium
In the selenium batch experiment of Tolu et al. [32], the range of experimental K d lies within one order of magnitude (Fig. 4). The calculated K d was initially at least one order of magnitude too low. This was most probably due to the high phosphorus concentration (3.2 mg/L) assumed for the pre-equilibration solution. This number, derived from Tab. 5.2 given in [18] is the median value for soils that are being fertilized. Gradually decreasing this value gave a minimum number of residuals at a value of 0.06 mg/L which is a value for inorganic phosphorus that is still in the range that can be found in soils [48,49]. After this adjustment, most of the calculated values deviate by a factor of less than 2, the largest variation being a factor of 2.9. In all cases, the model predicts that ≥90% of the sorbed  [50,51]; the interaction of Se with organic acids is reported to be "not well understood and remains controversial" [52]. In this case, the partitioning of Se can obviously be described reasonably well without a model for reactions with OM. However, under the actual conditions in the soil and at larger time scales this may be different, as the rate of reaction of OM seems to be a slow process [53,54]. Moreover, reduction processes (e.g. SeO 3 2− to Se 0 ) take place that are mediated by microbial action [16]. The extent of these effects regarding Se partitioning is not yet clear.

Quantitative analysis
To evaluate the results from the previous sections, the two quantities F and F' are defined:  Table 2. In the case of Pu, two batch experiments have been discarded where the resulting solutions were most probably not in equilibrium (see Table 3 in [35], soils CO-B and WA-A).
On average, the model tends to underestimate the results for Am and to overestimate those for Pu by a factor of less than 2, while for Se, it matches the experimental average quite closely. The IAEA values on the other hand  is no discernible difference, which is probably due to the fact that important parameters like Fe and Al in solution as well as DOC had to be estimated. The importance of these parameters is evaluated in the following section. Even so, for Pu, the predictions of the model are on average not worse than the use of the IAEA best estimate and for Am and Se they are apparently superior. The average deviation of the simulations relative to the experimental values is ten to the power of F' sim ; for all three elements, these values are lower than 2.5.
For K d estimations in the field, the predictive power of the model is probably better, because in the experiment, perturbing effects (e.g. drying, storing and subsequent dilution) may occur in a way that is not the case in the field and that cannot be quantified by the model. For testing this hypothesis, appropriate data would have to be available. On the other hand, the best K d estimates from the IAEA have mostly been derived from batch or column experiments that do not necessarily represent the conditions in the field.
Of course, a reliable simulation requires that all important soil parameters have to be known. This has been one of the reasons for performing the sensitivity study.

Sensitivity
The next three tables show the minimum and maximum values of the distribution coefficients of Am, Pu and Se calculated in the given parameter ranges. Only the most significant results are included where the ratio defined by R = K d,min /K d,max is higher than 2. The parameter ranges are those given by Blume et al. [18] for frequently occurring values in agricultural soils. Figures are only included if the parameter in question is included in the tables and if the K d dependence on this parameter is strongly nonlinear. In the case of Se and inorganic P, the results for very small P concentrations are also shown because the K d is highly sensitive to variations of this parameter. For pH, the range is extended to include all the soils used for validation; pe values were between 0 and 12 to cover the range between anoxic and moderately oxic conditions frequently found in soils. It has to be emphasized that in many cases, soil parameters are correlated (e.g. pH and pe); this has not been taken into account in the simulation. The calculations have also been carried out for Refesol 03G, which is a soil that has significantly higher amounts of clay, hydrous oxides and OM. While the R values may be slightly different in some cases, the conclusions regarding parameter importance are the same. Thus, these results are not shown here.

Americium
The most important parameter determining the partitioning of americium is pH, where the R value is 3.8. Varying over an order of magnitude, K d (pH) shows a maximum at a value of about 7.1 (Fig. 5); this coincides with the findings of Bruggeman et al. [45], who experimentally analyzed Am sorption to organic matter extracted from Boom clay. Above pH 5, the americium in the solid phase is predominantly bound to OM (>98%). At lower pH values, illite becomes increasingly important; at the same time, K d decreases with DOC content, because at higher concentrations more Am is kept in solution being bound to dissolved organic matter. The distribution coefficient increases by a factor of two if the Al content is increased from 1 to 5 mg/L, which is due to the competitive effect of Al sorbed to dissolved organic matter. However, as these high Al contents are usually only found in very acid soils [18], this effect is likely to be negligible in the range considered here (see Fig. 5). In the range below an Al concentration of ~1 mg/L, which is predicted by the simulations, the K d variability is low; this is why the effect of the Gibbsite saturation index correction is small.
Overall we can conclude that in most agricultural soils, the dominant parameters for Am partitioning are pH and DOC (Table 3).

Plutonium
The K d for Pu shows a strong dependence on pH above values of 6.5 and a maximum at pH ≈ 7.7 (Fig. 6). This is the region where 99% of the Pu is bound to Illite and only one third of the Pu in solution is bound to DOM. At pH = 5.6, 75% of the Pu is bound to Illite, 25% is bound to OM. At this pH, ca. 98% of the Pu in solution is bound to DOM; this percentage strongly decreases with increasing pH to values below 5% for pH > 8. The model also predicts for RefeSol 01A that in the region of neutral to alkaline pH, Pu sorption to hydrous oxides also plays a rôle when pe is above ~9, but for acid to neutral conditions below pe 9, less than 1% of the sorbed Pu is present on these minerals. As the pe approaches oxic values, the K d maximum disappears (Fig. 6).
If the redox voltage is varied, the distribution coefficient exhibits a plateau between pe 4 and 11 with a value of 430 L/ kg (see Fig. 7). This is the region where Pu(IV) is the dominant sorbing species; this range is extended and shifted to a lower pe range at pH 8 while it is shortened and shifted to a higher pe range at pH 4.5, according to the stability range of Pu(IV) in the Pourbaix diagram; see e.g. [55,56]. Pu(III) and Pu(VI) dominate at low resp. high pe values and are considered to be less relevant in most natural systems [57]. Even so, these species are included in the sorption model. Their weaker binding to the soil components is reflected by the smooth decrease of K d on both ends of the pe range. According to the model, the K d is quite insensitive to pe changes in the region that is characteristic for the soil subsurface. If the soil becomes waterlogged, pe decreases and the model predicts that eventually Pu will be mobilized. The complexity of the dependence on pH and pe of the plutonium K d is also shown in two three-dimensional plots in the appendix (Fig. 12).
The soil solution speciation of iron is determined by DOM-bound Fe while also colloid Fe nanoparticles may occur [39][40][41]. The omission of suspended Fe colloids in this model is probably the reason why the calculated values for Fe in solution are strongly underestimated [48]. However, as the amount of free and DOM-bound Fe calculated in the model will not depend on whether the Fe mineral controlling these values is suspended in solution, the resulting distribution coefficients will be unaffected. There is a strong correlation between K d and the sum of inorganic Fe and Fe-DOM in solution. The K d shows a large maximum around Fe concentrations between 0.01 and 0.1 mg/L (Fig. 8). In this region, 90-96% of the sorbed Fe is predicted to be present on clay mineral surfaces. At decreasing Fe concentrations, the percentage of DOM-bound Pu increases, thus keeping the Pu in solution and decreasing the K d . A similar effect is found when the DOM concentration is increased (this has also been suggested by Kersting [56]), as was the case for americium. At high Fe concentrations in solution, the competition of Fe for binding sites for Pu will be high, again leading to lower K d values. The change of partitioning with respect to the different components of the surface assemblage is shown in Fig. 11 in the annex. Similar considerations apply to Al, where these effects are only relevant at high Al concentrations (>0.5 mg/L), because sorption is only taken into account for OM and DOM, where Al is less tightly bound in comparison to Fe. The strong sensitivity of K d to Fe concentration in solution is obviously a major source of uncertainty in the estimation of Pu partitioning if Fe is estimated by equilibrating the solution with a single mineral.
Consequently we find that in the case of plutonium, its partitioning depends on the complex interplay of the soil parameters pH, pe, Fe, Al and DOM (Table 4).

Selenium
The distribution coefficient for selenium strongly depends on the amount of phosphate in the soil solution (Fig. 9). This is caused by the strong complexation of phosphate in comparison to Se by hydrous ferric oxide surfaces [58], which is reflected by the respective complexation constants in the data base. This results in a very low selenium K d in fertilized soils which also has been reported by Eich-Greatorex et al. [59]. The model predicts that at high phosphate concentrations, illite is the dominant sorbent, while at P < 0.1 mg/L, most of the selenium (>95%) is bound to hydrous oxide surfaces.
The pH dependence of K d exhibits a plateau from acid to nearly neutral values (Fig. 10); at low phosphate concentrations, the curve has a moderate positive slope in that region. As the pH approaches the point of zero charge (~pH 8 [60]) of ferrihydrite, the K d decreases drastically. The hydrous oxide/illite sorption ratio generally decreases with increasing pH.
Similar to Pu, the K d shows a plateau between pe 3 and 9, which is the region where selenite is most stable [61]. At lower pe, Se (0) and selenide dominate; in our calculations we assume that elemental selenium stays suspended in solution, while for selenide, complexation constants have been derived by the linear free energy relationships for anions given in [25], taking the proton dissociation constant of HSe − given in the Thermochimie database. This explains that the K d does not reach zero in the calculations. However, predictions for Se partitioning in this region should be treated with particular caution because it has been reported that selenide ions can be sorbed to iron selenide surfaces, such as Pyrite [62] or Mackinawite [63], which may be present in a soil under anoxic conditions. They are not yet included in the model because data on the amount of these minerals as well as corresponding suitable sorption models are not available to date. At higher pe, the K d decreases because the dominating species selenate is less tightly bound to solid surfaces than selenite [22,64,65]. The plateau region is shifted towards lower pe when pH is increased and vice versa (see Fig. 13 in the appendix).
Thus it can be concluded that in most agricultural soils, the partitioning of Se is governed by a combination of the parameters phosphate in solution, pH and pe (Table 5).

Conclusions
In this work, the model UNiSeCs II is described and it is applied for the estimation of speciation and partitioning of three long-lived radionuclides (Am, Pu and Se) whose physicochemical properties in agricultural soils are considerably different. The model is consistent in the sense that for all three nuclides, the same set of sorbents and general assumptions is applied. It is shown that experimental distribution coefficients from laboratory experiments can be approximated reasonably well, even if the relevant soil solution parameters are not known and have to be estimated. A comparison with the best estimates from the IAEA indicates that UNiSeCs II has the potential to significantly improve K d estimations for agricultural soils. The average deviation of the simulations relative to the experimental values are lower than 2.5 for all three elements. However, to quantify this reliably, more experimental studies are required that include the necessary amount of data for modeling.
For each radionuclide, the calculated distribution coefficient in a single representative soil (Refesol 1A) has been analyzed for its dependence on soil solution parameters. The most important parameters are pH in all cases, DOC for the two actinides, pe for Pu and Se, Fe and Al in solution for Pu and phosphate in solution for Se. In some cases, the ratios between the calculated maximum and minimum K d values are larger than 100 (e.g. pe in the range 0-12 for Se), but almost constant if the parameter range is reduced (e.g. to pe 3-9 for Se). On the other hand, the K d variation may be large in a comparatively small parameter range (e.g. a factor of 6 at pH 6.6-7.2 for Se). Thus, it is important to realize that even for one particular soil, the K d may vary considerably if different sets of soil solution parameters are assumed. Moreover, as the model is multidimensional and parameters may not always be independent, it cannot be excluded that other parameters become relevant in certain regions of parameter space. This is a subject for further study.
It has also to be pointed out that for now, the model is limited to soils where the components clay, hydrous oxides and organic matter are the dominant sorbents. Nonetheless, if thermodynamical sorption data are provided, extensions for UNiSeCs II are possible that take into account additional sorbents.
Currently, UNiSeCs II is being applied for calculating look-up tables of distribution coefficients similar to the "smart K d " concept [66], containing values for possible combinations of the relevant soil parameters. Within the ongoing joint project Trans-LARA, these tables are intended to be used for hydrological simulations of radionuclide migration from groundwater to the upper soil horizon and also for biosphere simulations with Ecolego [1] to estimate the radiation burden of the human population after the contamination of agricultural soils. In the next step, the model will be re-evaluated for uranium, nickel and cesium and it is planned to include other radionuclides as well. According to the parsimony principle, the model will be kept as simple as possible but as complex as necessary.

Appendix
Complexation and exchange submodels for Am, Pu, Se and competing ions in UNiSeCs II (Table 6).
Additional figures from the sensitivity analysis (Figs. 11, 12, 13).  Availability of data and material All data used in this article were taken from the cited literature.