Searching for optimum adsorption curve for metal sorption on soils: comparison of various isotherm models fitted by different error functions

Studies comparing numerous sorption curve models and different error functions are lacking completely for soil-metal adsorption systems. We aimed to fill this gap by studying several isotherm models and error functions on soil-metal systems with different sorption curve types. The combination of fifteen sorption curve models and seven error functions were studied for Cd, Cu, Pb, and Zn in competitive systems in four soils with different geochemical properties. Statistical calculations were carried out to compare the results of the minimizing procedures and the fit of the sorption curve models. Although different sorption models and error functions may provide some variation in fitting the models to the experimental data, these differences are mostly not significant statistically. Several sorption models showed very good performances (Brouers-Sotolongo, Sips, Hill, Langmuir-Freundlich) for varying sorption curve types in the studied soil-metal systems, and further models can be suggested for certain sorption curve types. The ERRSQ error function exhibited the lowest error distribution between the experimental data and predicted sorption curves for almost each studied cases. Consequently, their combined use could be suggested for the study of metal sorption in the studied soils. Besides testing more than one sorption isotherm model and error function combination, evaluating the shape of the sorption curve and excluding non-adsorption processes could be advised for reliable data evaluation in soil-metal sorption system.


Introduction
Metal sorption by soils is of particular interest and having importance from both agricultural and environmental viewpoints. Adsorption is a major process responsible for the retention of metals by soils [5]. For the exact characterization of the sorption process, it is crucial to provide the most appropriate information on the adsorption equilibrium. This information is fundamental for the reliable prediction of adsorption parameters and comparing the adsorbents' behavior in varying sorption systems quantitatively [16]. Adsorption isotherms are generally used to describe how adsorbates interact with the adsorbent. Their use is crucial for several reasons, like the study of adsorption mechanism pathways, surface properties, and adsorbents' capacities [12]. In isotherm studies, the experimental adsorption data are inserted into an equilibrium model to achieve the best fit for the sorption system. The better the model's fit to the experimental data, the more accurate the calculation can be [10].
For soil-metal systems, the metal sorption capacity is one of the most critical parameters that have to be determined. Besides sorption capacities, equilibrium models often supply insight into the adsorbent's affinity, its surface properties, and the adsorption mechanism [6]. Besides the well-known Langmuir and Freundlich models, many other isotherm equations can be used to explain soil-metal adsorption systems [3].
Traditionally, isotherm parameters are determined using the linear least-squares method [13]. Nevertheless, a considerable inconvenience related to the linearization of the isotherm models has been identified. According to Hong et al. [22], linearization may produce several different outcomes, and it may alter the error structure strongly. Additionally, it affects both the error variance and the normality assumptions of standard least-squares disadvantageously, leading to the distortion of the data. As several linearization methods exist, their parallel use results in the change of error distribution differently. That is why non-linearized regression analysis became unavoidable as it provides a mathematically rigorous method for determining adsorption parameters using the original form of isotherm equations [18]. Unlike linear regression, it usually involves minimizing the error distribution between the experimental data and the predicted isotherm based on its convergence criteria [44].
Numerous error analysis methods have been used and developed to determine the best-fitting isotherm equation [48]. The most popular tool for model comparison is the use of the coefficient of determination (R 2 ). Despite the apparent simplicity of interpreting R 2 , it is not always suitable for evaluating the goodness of fit and comparing different models. According to El-Khaiary and Malash [13], R 2 is sensitive to extreme data points, influenced by a range of the independent variables, and can be artificially large by adding more parameters. However, in the case of nonlinear methods, previous studies reported that the predicted isotherms were varying with different error functions. Consequently, comparing the usefulness of different error functions should precede selecting the best fitting isotherm [24].
Comparison of a large number of isotherm models with different error functions was carried out on several single-component natural adsorbent-adsorbate systems in the last decade (e.g., [6,17,21,37]. However, such studies are lacking completely for soil-heavy metal adsorption systems. This shortcoming could be primarily due to the soil's substantial heterogeneity that provides several different active mineral and organic surfaces for metal ions. Additionally, precipitation of metal ions may occur during the equilibrium adsorption experiments, making it challenging to use isotherm models based on equilibrium data [38]. This study aimed to fill the gap in comparing adsorption models with different terror functions for the soil-metal system. As the metal sorption capacity is the most important parameter for soils, we studied whether models beyond the widely used Langmuir model support more reliable sorption capacity values. The adsorption equilibrium isotherms of four widely studied heavy metals in four soil samples with different geochemical properties were investigated. A trial-and-error method of fifteen two-, three-, four-, and five-parameter nonlinear models were used altogether to determine the best fitting isotherm. Additionally, seven error functions were also used in the assessment of the adequacy and goodness-of-fitted models. We tried to find the most useful isotherm model and error function combination(s) to predict metal adsorption performance in soils.

Soil samples
As the most important soil components affecting metals' sorption are pH, organic matter, clay minerals, and Fe-oxides [43], these soil properties were considered primarily for the sample selection. Two acidic and two alkaline soil samples were used from a Luvisol and a Phaeozem profile, respectively. The samples were studied for their pH (in 0.1M CaCl 2 ), total organic carbon content (Tekmar-Dohrman Apollo 9000N TOC analyzer), BET surface area (Quantochrome Autosorb-1-MPV N 2 gas sorption system), cation exchange capacity (ISO 23470:007), particle size distribution (Fritsch Analysette Microtech A22 laser diffraction instrument), and dithionite extractable iron content (after Mehra and Jackson [28] using Perkin-Elmer AAnalyst 300 atomic absorption spectrometer). As shown in Table 1, the selected samples show significant differences in their total organic carbon, clay, and dithionite-extractable-iron content.

Adsorption experiments
Competitive batch adsorption experiments were carried out with the metals Cd, Cu, Pb, and Zn. These metals are often related both to environmental and agricultural problems. For example, Cd and Pb are known for their toxic effects on soil biota, whereas Cu and Zn often show a deficiency in soils requiring artificial fertilization [2]. These metals also exhibit contrasting sorption properties in soils since Cu and Pb are retained through inner-sphere complexes, whereas Cd and Zn through outer-sphere complexes predominantly [36].
The experiments were carried out in duplicates. Soil:solution ratio of 1:30 was used, and the initial concentration of metals were set to 0.1, 0.2, 0.5, 1, 2, 5 and 10 mmol/L in 0.01M Ca(NO 3 ) 2 solution. The initial solution's pH was set to 5.5 to avoid metal hydroxide precipitation [47]. Equilibration was carried out for 24 h at 22 °C. Then the soil and the solution were separated by centrifugation at 4000 rpm for 20 min. Solutions were then filtered, and the concentrations of metals were analyzed by atomic absorption spectrometry (Perkin-Elmer AAnalyst 300). The relative standard deviations of duplicate analyses are below 5% for each metal at equilibrium concentrations (C e ) above 100 mg/L and never reached 10% even at lower concentrations.
The concentrations of metals adsorbed by the soil samples were calculated using the following equation: where Q e is the sorbed metal amount per unit weight of the soil (mmol/kg), C e is the equilibrium metal concentration in the solution (mmol/L), C i is the initial metal concentration in the solution (mmol/L), V is the volume of the solution (mL), and W is the weight of the air-dried soil (g).

Determining sorption curve parameters by nonlinear regression
The sorption curve parameter sets were determined by nonlinear regression. Four two-parameter, nine threeparameter, one four-parameter, and one five-parameter model were used to study the metals' sorption performance. Models providing maximum sorption capacity parameters were selected primarily for this study (Table 2). First, seven minimizing methods were used to solve the isotherm equations by minimizing the errors between the models and the experimental data (Table 3). A trial-and-error method was used for the regression analysis to reduce the objective function using the solver add-in function of the MS Excel 2016 Akaike's Information Criterion (AIC) was used to compare the results of minimizing procedures and the fit of sorption curve models. It is a well-established statistical method based on information theory and maximum likelihood theory. It determines which model is more likely to be correct [13]. For a small sample size, AIC is calculated for each model from the equation: where n is the number of data points, SSE is the sumof-squared deviations of the points from the regression curve, and n p is the number of parameters in the model. Freundlich Foo and Hameed [14] Dubinin-Radushkievich Foo and Hameed [14] Jovanovic Q e = Q mJ ⋅ (1 − exp K J ⋅ C e Panahi et al. [32] Three-parameter models Hamdaoui and Naffrechoux [19] Langmuir-Freundlich Hamdaoui and Naffrechoux [19] Radke-Prauschnitz Hamdaoui and Naffrechoux [19] Tóth Hamdaoui and Naffrechoux [19] Vieth-Sladek Hill

Statistical analyses
Statistical calculations were carried out to show whether both different error functions and sorption curve models provide significantly different AIC values. In all cases, the values of asymmetry (skewness) and kurtosis between −2 and 2 were verified to assume normality when the D'Agostino-Pearson test was significant. The results showed that data sets based on the AIC values of the studied error functions do not exhibit a normal distribution. When the data sets were studied separately according to the sorption curve types, data of subtype-1 and subtype-2 showed still non-normal distribution, whereas data of subtype-mx showed normal distribution. On the contrary, the AIC values of the studied sorption curve models exhibited normal distribution mostly, except the data of Dubinin-Radushkievich and Brunauer-Emmett-Teller models for the whole data set, and the data of Sips, Langmuir-Freundlich, and Hill models for the data of subtype-1 curves. Homogeneity of variances was tested through Levene's test and verified for nearly all cases, except when the data sets were studied separately for the sorption curve subtypes. As the ERRSQ error function generally provided the lowest error between the predicted and experimental data, its data set was also studied separately with statistical analyses. This data set showed the normal distribution in each case with homogeneous variance. The one-way analysis of variance (ANOVA) was used for data sets with normal distribution and homogeneous variance. This test shows if there is a significant difference between the AIC values provided by different error functions or models. To discover between which datasets the significance can be found, the Tukeytest was used as a post-hoc analysis. For data sets with non-normal distribution and/or non-homogeneous variance, the Kruskal-Wallis-test was used combined with the Mann-Whitney-test as a post-hoc test. The significance level of each test was set to = 0.05. The statistical analyses were carried out using the StatistiXL add-in of the Microsoft Office Excel 2016 software.

The shape of the sorption curves
The shape of the sorption curves is generally evaluated to predict specific immobilization processes. Different sorption models may be able to describe sorption curves differing in their shapes [4]. Hence, it is essential to specify the sorption curves of the metals before model fitting.
The classification of the curves was carried out according to Giles and Smith [15]. We found that the L1-type curve is the most common in acidic samples (Fig. 1). This curve type was found for Cd in the acidic samples, for Pb in the S1 sample, and for Cu and Zn in the S2 sample. Copper in the S1 and Pb in the S2 samples exhibited L2-type curves, whereas an Lmx-type curve was found for Zn in the S1 sample. Contrarily, an H1-type curve was observed for Pb in both alkaline samples, whereas Cu exhibited H2-type curves in these samples, and this was the case for Cd in the C1 sample. Zinc showed Hmx-type curves in both alkaline samples similarly to Cd in the sample C2. Sorption curves of metals in the acidic samples belong to the main class "L". In contrast, those in the alkaline ones belong to the "H" class, which refers to Langmuirtype and high-affinity classes, respectively. A decreasing slope is characteristic of the L-type isotherm as concentration increases since vacant adsorption sites decrease Ho [20] Sum of the squares of the errors Hadi et al. [18] Hybrid fractional error function (HYBRID) Xue et al. [50] Derivative of Marquardt's Percent McKay et al. [29] The Average Relative Error (ARE) Piccin et al. [33] Sum of the Absolute Errors EABS) Dotto et al. [11] The chi-square statistic ( X 2 ) McKay et al. [29] as the adsorbent becomes covered. The adsorbate's high affinity could explain such adsorption behavior for the adsorbent at low concentrations, which then decreases as concentration increases. For the H-type curves, the adsorbate shows (almost) complete immobilization by the adsorbent suggesting especially high adsorbatesubstrate affinity. Nevertheless, both L-and H-type curves are produced when strong attractive forces exist between the adsorbate and adsorbent, but very weak forces between the adsorbates themselves [42]. Curves in subtype-1 represent systems in which the adsorbate monolayer has not been completed, probably because of experimental difficulties. This rather characteristic in the acidic samples and for Pb in the alkaline samples. The phenomenon can generally be due to the well-known concentration effect [49] on the one hand. It can also suggest the adsorbents' precipitation, primarily when the curve exhibits extremely high slope even at high initial concentrations, respectively [39]. In the subtype-2, a plateau can be identified, representing the completion of the first monolayer. This subtype was characteristic both in the acidic and alkaline samples and generally enabled the correct estimation of the adsorbate's maximum adsorption on the adsorbent [25]. In the subtypemx, the sorption curve has a maximum, which suggests desorption from the surface, and this is also often related to the adsorbent concentration effect.

Comparison of the error functions
The best-fitting isotherm is generally selected based on the error function(s) that produce minimum error distribution between the predicted and experimental isotherms. Several studies reported that the predicted isotherms were varying with the error function when nonlinear methods were used (Kumar and Porkodi 2014). Thus, to analyze the impact of various error functions on the predicted isotherms, seven different error functions were optimized. They are among the most widely used ones, and they have been often developed for adsorption studies [3]. Each error function supports an objective function for the minimum error distribution between experimental and predicted isotherms. First, we have ranked the error functions from the best to worst based on their AIC values. Rank 1 refers to the lowest AIC value (e.g., showing the lowest error distribution between the predicted and experimental curves), whereas rank 7 refers to the highest one (Fig. 2). ERRSQ provided the lowest error distribution for most of the models for each soil-metal pairs, as it was ranked in category 1 at 96% of the cases. If ERRSQ was ranked in another class, it still resulted in similarly low error distribution to the error functions ranked in category 1. We have not found any relationship between the lower ranking of ERRSQ function and the type of the metal, the soil, the sorption curve, or  [15]. C e = equilibrium metal concentration in the solution (mmol/L), Q e = adsorbed metal concentration in the soil (mmol/kg) the sorption model. In the ranking, ERRSQ was followed by R 2 (generally ranked in category 2), MPSD, and HYBRID (ranked mostly in categories 3 and 4, respectively). Contrarily, EABS provided the highest error distribution, as shown by its ranking in category 7 at 52% of the cases. Similarly, high error distribution was shown by the ARE function, which was ranked mostly in categories 6 and 7. The X 2 error function showed a uniform distribution among the ranks between 2 and 6 with frequency ratios from 18 to 22%.
The ranking was carried out separately also for the groups of soil-metal pairs with different sorption curve subtypes. In this case, sorption curves belonging to subtype-1 and -2 showed similar results to those presented above (Fig. 2b, c). Slight differences were found for R 2 , which exhibited higher frequencies in the rank 2 for both subtype-1 and -2 curves. The HYBRID function was ranked in category 5 for the soil-metal pairs with adsorption curves belonging to subtype-2, and the ARE function showed higher frequencies in the ranks 6 and 7 there. On the contrary, soil-metal pairs with sorption curves belonging to subtype-mx exhibited significantly different ranking characteristics for the error functions (Fig. 2d). Although the ERSSQ function could be ranked in category 1 at 98% of the cases, it was followed by the MPSD and X 2 functions, both often ranked in category 2. The HYBRID function still showed similar characteristics to those presented above, and it was ranked mostly in category 4. However, R 2 resulted in much lower fits in this case, and it was ranked in category 6 mostly. The error functions ARE and EABS showed similar characteristics again to those found for the other sorption curve subtypes. The above data showed that ERRSQ provided the lowest error distribution between the experimental and predicted sorption curves for almost every case. Shahmohammadi-Kalalagh and Babazadeh [37] also found that this error function was the best to minimize the error distribution between the experimental data and predicted two-parameter models (Freundlich and Langmuir) for the Zn-kaolinite and Cu-kaolinite sorption systems. However, ranking provides only a relative comparison of the studied error functions, and it does not inform us whether the differences are statistically significant between them. That is why we used statistical analyses (ANOVA and Kruskal-Wallis-tests) to compare the statistical parameters of the AIC values calculated for each error function. These analyses showed that there is a significant difference between the AIC values of the studied error functions. Such analyses were also carried out separately for the different sorption curve types. It was found that there is a significant difference within the AIC values of the error functions for subtype-1 and subtype-2 curves but not for subtype-mx curves. However, ANOVA and Kruskal-Wallis-tests are not able to tell between which datasets the significance can be found. So, Tukey and Mann-Whitney-tests were carried out, respectively, as post-hoc analyses between the AIC values of each two error functions both for the whole data set and for the sorption curve subtypes 1 and 2 separately. For the entire data set, significant differences were observed between the AIC values of the EABS error function and those of all other error functions. Additionally, significant differences were found between the AIC values of ERRSQ, R 2 , and the ARE error functions. For the subtype-1 sorption curves, significant differences were observed between the AIC values of the EABS error function and those of almost all other error functions, except the ARE error function. However, for the subtype-2 sorption curves, significant differences were found only between the AIC values of the EABS error function and those of the ERRSQ, R 2 , and MPSD error functions, and between the AIC values of the ARE and ERRSQ and R 2 functions. Consequently, besides the ERRSQ function, even R 2 , HYBRID, MPSD, and X 2 provided similarly low error distribution between the experimental and predicted sorption curves for all of the studied cases. However, the EABS error function mostly provided a significantly higher error distribution than the other error functions. Additionally, this was the case for the ARE error function compared to the ERRSQ (and R 2 ), which generally provided the lowest error distribution.
Some studies found that different error functions have to be used to minimize the error distribution between the experimental and predicted data for sorption models with a different number of parameters [37]. Typically, if the number of model parameters increases, the fitting error decreases because a better fit can be obtained.
Simultaneously, the generalization bias is expected to increase due to the more considerable model variability [35]. However, smaller models may tend to do better for the same data set. Thus, when models with different parameters are fitted on the same set of data points, the model with a sufficient number of data points in excess will provide the best conformity [30].

Comparison of the sorption curve models
The adsorption curve is an invaluable curve that describes a substance's retention from the aqueous media to a solid phase. Together with the underlying thermodynamic assumptions, its physicochemical parameters provide an insight into the adsorption mechanism, surface properties, and the degree of affinity of the adsorbents [9]. Several equilibrium adsorption models have been developed in the last century regarding the fundamental approaches, such as kinetic considerations, thermodynamics, and potential theory [27]. Recently, the isotherm modeling trend is the derivation in more than one approach, thus directing towards the difference in the physical interpretation of the model parameters [14].
Sorption curve models were also ranked based on the AIC values similar to the error functions. Their ranking was carried out based on the AIC values provided by the ERRSQ function, as this function provided the lowest error function generally. The results of the ranking are shown in Table 4. The sorption curve models could be categorized into six groups based on their ranking characteristics. Fritz-Schlunder, Baudu, and MacMillan-Teller models exhibit very low ranking generally. The Brunauer-Emmett-Teller and Vieth-Sladek sorption models could be characterized by a relatively low ranking in most of the cases. On the contrary, the Hill and Langmuir-Freundlich models showed a rather high ranking in most of the cases, whereas the Brouers-Sotolongo and Sips models exhibited the highest rankings generally. The Jovanovich, Langmuir, and Dubinin-Radushkievich models showed very high rankings for sorption curves belonging to the subtypemx. However, they presented rather low rankings for all other sorption curve types. Finally, the Freundlich, Radke-Prauschnitz, and Tóth models have not shown any specific characteristics in their rankings.
Considerable variation was found among the different sorption curve types fitted by different sorption models. However, their ranking based on their AIC values allowed us to specify which models provide a good fit for most cases and which ones do not. The ANOVA and Kruskal-Wallis-tests were used again to study whether the studied sorption models provide significantly different fits to the experimental data based on their AIC values. They showed that there is a significant difference between the AIC values of the different sorption curve models. The difference was observed both for the entire data set and for the data of varying sorption curve types. Therefore, we used post-hoc tests (Tukey and Mann-Whitney-tests) to find between which two data sets can be the significance found. The post-hoc tests showed that the MacMillan-Teller, Baudu, and Fritz-Schlunder models provided significantly different AIC values than the other sorption models for the whole data set. Additionally, the Dubinin-Radushkievich and the Brunaur-Emmett-Teller models also showed significantly different AIC values to those of Langmuir, Sips, Langmuir-Freundlich, Tóth, Brouers-Sotolongo, and Hill models. For the sorption curves belonging to the subtype-1, the Mac-Millan-Teller and Fritz-Schlunder models exhibited significantly different AIC values to the other models, and so did the Baudu, Dubinin-Radushkievich, and Jovanovich models mostly. On the contrary, the Freundlich, Sips, Langmuir-Freundlich, Tóth, Brouers-Sotolongo, and Hill models have not provided significantly different AIC values to each other. For the sorption curves in the subtype-2 group, the Mac-Millan-Teller, Baudu, and Fritz-Schlunder models provided significantly different AIC values than all other models. On the other hand, the Dubinin-Radushkievich, Brouers-Sotolongo, Langmuir, Sips, Tóth, Hill, and Langmuir-Freundlich models exhibited similar AIC values, which are significantly different from those of the other models in most cases. For the sorption curves belonging to the subtype-mx, besides the MacMillan-Teller, Baudu, and Fritz-Schlunder models, the Langmuir, Dubinin-Radushkievich, and Jovanovich models exhibited significantly different AIC values from all other sorption models. This comparison is in good harmony with the ranking of AIC values of the sorption models. Models with very low ranking (MacMillan-Teller, Baudu, and Fritz-Schlunder) provided significantly different AIC values to other models almost in each case. The Sips, Langmuir-Freundlich, Brouers-Sotolongo, and Hill models exhibited high rankings generally, and they provided similar AIC values. Additionally, they can complement the Tóth and Freundlich models in some cases (e.g., for subtype-1 curves). The Langmuir, Dubinin-Radushkievich, and Jovanovich models exhibited high ranking only for sorption curves Table 4 Sorption curve models ranked from the best to worst based on their AIC values calculated using the ERRSQ error function. Green and red background colors represent the best and worst fitting performances, respectively. The orange background color shows models exhibiting very good performances for type-mx sorption curves primarily.  Consequently, there can be at least 4-7 models providing significantly not different fit to different sorption curve types. The maximum sorption capacity (Q max ) is one of the most important parameters when studying metals' sorption performance in soils. Therefore, we also compared these parameters calculated by the different sorption models. The Kruskal-Wallis-test showed that the Q max values exhibited no significant difference when calculating them based on the different error functions with similar performance in minimizing the error. Therefore, the Q max values calculated using the ERRSQ error functions will be compared below. The averages and ranges of the Q max values are shown in Fig. 3. The coefficient of variation (CV) of the Q max values was generally low (below 0.4) in most cases. Exceptions are the Cd sorption by the samples S1 and S2, where extraordinarily high Q max values were calculated using the Tóth model. If these values are not taken into account, the CV values of the Q max values are below 0.6. Additional exceptions are the sorption of Pb on the C1 and C2 samples. In these cases, several models, like Tóth, Sips, Hill, Langmuir-Freundlich, and Brouers-Sotolongo, supported unrealistic high Qmax values, probably due to the almost complete retention of Pb in these samples. The mineralogical study of these samples showed that Pb precipitated as cerussite (PbCO 3 ) during the experiment. The precipitation of Pb resulted in a unique sorption curve shape [40]. However, the other models provided quite similar Q max values with relatively low variance (with CV < 0.3). Generally, the Tóth model resulted in the highest Q max values, often resulting in outlier values. The Q max values calculated through the Tóth model was usually higher than the upper quartile of all the Q max values. This phenomenon was also observed for the Sips, Langmuir-Freundlich, and Hill models. On the contrary, the lowest Q max values were calculated using the Baudu sorption model. In more than 60% of the cases, the Q max values calculated by this model were lower than the lower quartile of all these values. Similar characteristics were found for the Radke-Prauschnitz model. As expected, the Freundlich model also suggested lower sorption capacities. However, the Freundlich K F value is only equal to the maximum adsorption capacity when n approaches infinity [26]. Although K F is not defined as the maximum adsorption capacity, K F values should be of the same order as Q max values [45]. This phenomenon was also found in this study. On the contrary, the Langmuir, Dubinin-Radushkievich, Jovanovich, Brouers-Sotolongo, Brunauer-Emmett-Teller, and Fritz-Schlunder models provide Q max values falling in the interquartile range mostly. We have not found any relationship between the variation of Q max values and the sorption curve types. Our results showed that different sorption models might provide some variation in their fit to the experimental data. On the one hand, several models generally offer a very similar estimation for the soils' maximum sorption capacity. On the other hand, some models often overestimate or underestimate the maximum sorption capacities systematically, like the Tóth and Baudu models, respectively. Moreover, several sorption models may estimate unrealistic high maximum sorption capacity for metals with sorption curves lacking developed plateau at least partially. If a process other than adsorption results in this sorption curve, sorption models can not be used to calculate Q max values. If such processes can be excluded, experimental conditions have to be re-designed [12]. For the studied soils, the Brouers-Sotolongo model was found to be the most useful one. This model was always among the models showing the best fit to the experimental data for each sorption curve type. Besides this model, other ones like the Sips, Hill, and Langmuir-Freundlich models also exhibited good performances. Other models can also be useful for certain sorption curve types. A further advantage for the Brouers-Sotolongo model is it provided realistic maximum sorption capacities for the studied soil-metal pairs. This model is a combined form of the Freundlich and Langmuir expressions. It is deduced for predicting the heterogeneous and microporous adsorption systems. Additionally, it can measure the width of the adsorption energy distribution and energy heterogeneity of the sorbent surface (Brouers and Sotolongo 2005). This model has not been widely used for the study of soil-metal interaction yet. However, several recent results suggest that the Brouers-Sotolongo model describes the metal sorption on several types of heterogeneous surfaces, like those that can also be found in soils. Examples are the adsorption of Pb on algae [7] that of As on activated carbon and MnFe 2 O 4 composite [34], that of Pb, Cd, Zn, Ni, and Cu on goethite-humic-acid-modified kaolinite [46], and that of Hg and Cr onto Fe 2 O 3 -SiO 2 composites [41].

Conclusions
Our results showed that the ERRSQ error function generally provided the lowest error distribution between the experimental and predicted sorption curves for the studied soil-metal pairs with varying sorption curve types and geochemical characteristics. Additionally, the R 2 , HYBRID, MPSD, and X 2 error functions also provided statistically similar error distribution in most cases. However, the EABS and ARE error functions provided significantly higher error distribution than the other error functions. So these latter error functions could not be suggested for fitting sorption curves in our case.
Using the ERRSQ error function, the Sips, Langmuir-Freundlich, Brouers-Sotolongo, and Hill models showed mostly very good fits for the studied soil-metal pairs with varying sorption curve types and geochemical characteristics. On the other hand, the Fritz-Schlunder, MacMillan-Teller, Brunauer-Emmett-Teller, Baudu, and Vieth-Sladek models generally showed worse fit. That is why they can not be suggested for the evaluation of metal sorption in the studied soils.
Based on the evaluation of the provided Q max values, the combined use of the Brouers-Sotolongo model with the ERRSQ error function could be suggested to characterize metal sorption in the studied soils. This model was always among the models showing the best fit to the experimental data for each sorption curve type. Besides this model, other ones like the Sips, Hill, and Langmuir-Freundlich models can also be suggested with the ERRSQ function, and further models can be useful for certain sorption curve types.
The results showed that several isotherm models and error function combinations could be found, which provide statistically similar estimation for the given soil's maximum sorption capacity. Besides testing more than one sorption isotherm model and error function combination, evaluating the shape of the sorption curve and excluding non-adsorption processes are advised for reliable data evaluation. Additionally, the goodness of fit of the models and the provided physicochemical parameters should be evaluated when selecting the optimum adsorption model.

Data Availability
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Compliance with ethical standards
Conflict of interest The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.