According to the QSAR model theory, activities or properties of pollutants (e.g., sorption coefficients) are explained by their chemical structure (Aranda et al., 2016; Olguin et al., 2019). This implies that the outcome is exclusively dependent on pollutants, or put another way, independent of soil properties.

In the literature, a high correlation between soil organic carbon content (OC) and sorption coefficients has been found for non-ionizable or neutral pollutants sorbed on soils with OC > 0.1% (Olguin et al., 2017; Olguin et al., 2019). This minimizes the soil variability (only OC is important) and allows the extrapolation of trends to other soils having the same sorption mechanism, especially if the mechanism is unspecific (i.e., pollutants may interact with different types of sorption sites) and generalizable (i.e., a specific soil component may represent the whole sorption), such as hydrophobic. Under this scenario, the sorption coefficient is a property of pollutants (Olguin et al., 2017; Cai et al., 2019).

In this sense, most QSAR models predicted K OC (Wang et al., 2015; Aranda et al., 2016; Olguin et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Zhu et al., 2017; Zhang et al., 2018; Cai et al., 2019; Olguin et al., 2019; Kobayashi et al., 2020; Kobayashi and Yoshida, 2021; Muhire et al., 2021; Pandey and Roy, 2021; Cantwell et al., 2022; Jiang et al., 2022) of non-ionizable and neutral compounds, such as diphenyl ethers and biphenyl congeners (Zhu et al., 2017; Zhang et al., 2018; Cantwell et al., 2022), perfluorinated and polyfluoroalkyl substances (Jiang et al., 2022), polycyclic aromatic hydrocarbons and phthalic acid esters (Cai et al., 2019), pharmaceuticals (Rybacka and Andersson, 2016; Berthod et al., 2017), pesticides (Kobayashi et al., 2020), organophosphorus insecticides (Daré et al., 2017; Muhire et al., 2021), or several organic compounds of different chemical classes (Wang et al., 2015; Olguin et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Olguin et al., 2019; Kobayashi and Yoshida, 2021; Pandey and Roy, 2021).

In contrast to QSAR, PTF theory considers that correlational associations among soil properties have explanatory implications for other soil properties (Sidoli et al., 2016; Singh et al., 2016; De Gerónimo et al., 2018; Conde-Cid et al., 2020). Therefore, sorption coefficients of one pollutant among different soils must depend on soil properties.

In the literature, different sorption mechanisms and trends have been found for the same pollutant, e.g., a major and a minor role of OC in the sorption of oxytetracycline and chlortetracycline across different soils (Conde-Cid et al., 2020) or non-hydrophobic sorption mechanisms (Dollinger et al., 2015; Kodešová et al., 2015; Sidoli et al., 2016; De Gerónimo et al., 2018; Klement et al., 2018), in accordance with the PTF theory.

Additionally, the outcome selection process was linked to the evaluation of sorption linearity during the data treatment (Kodešová et al., 2015; Sidoli et al., 2016; Klement et al., 2018) and the analysis of predictor variables related to the environmental and methodological conditions, such as tillage (Singh et al., 2016) or the presence of phosphate (Sidoli et al., 2016) in agricultural soils.

In this sense, PTF models studied one pollutant at a time, including pesticides (Dollinger et al., 2015; Sidoli et al., 2016; Singh et al., 2016; De Gerónimo et al., 2018), pharmaceuticals (Kodešová et al., 2015; Klement et al., 2018), and antibiotics (Conde-Cid et al., 2019; Conde-Cid et al., 2020; Hu et al., 2022), generally ionizable (Dollinger et al., 2015; Kodešová et al., 2015; Sidoli et al., 2016; De Gerónimo et al., 2018; Klement et al., 2018; Conde-Cid et al., 2019; Conde-Cid et al., 2020; Hu et al., 2022), and using non-normalized sorption coefficients ( K d and K F ) (Dollinger et al., 2015; Kodešová et al., 2015; Sidoli et al., 2016; Singh et al., 2016; De Gerónimo et al., 2018; Klement et al., 2018; Conde-Cid et al., 2019; Conde-Cid et al., 2020; Hu et al., 2022) and sorption isotherm linearity (Dollinger et al., 2015) for representing the linear ( 1 / n = 1) and nonlinear ( 1 / n ≠ 1) sorption under specific conditions.

We found different strategies that were employed to simplify the sorption process. QSAR models considered the sorption coefficient as a physicochemical property of pollutants, independent of soils and are therefore valid for all soils, simplifying the outcome selection (P1, Figure 3). PTF models evaluated one pollutant per model, considering that the sorption coefficient is a soil property. This implies that QSAR models, when applicable, are broad and general, while PTF models tend to be specific, acquiring local relevance. Figure 4 shows conditions and assumptions that we propose to understand both kinds of predictive models.

Physicochemical information of sorbents was typically not included in QSAR studies (Wang et al., 2015; Aranda et al., 2016; Rybacka and Andersson, 2016; Berthod et al., 2017; Daré et al., 2017; Olguin et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Zhu et al., 2017; Zhang et al., 2018; Cai et al., 2019; Olguin et al., 2019; Kobayashi et al., 2020; Kobayashi and Yoshida, 2021; Muhire et al., 2021; Pandey and Roy, 2021; Cantwell et al., 2022; Jiang et al., 2022). Additionally, methodological conditions used to obtain the experimental sorption coefficient values were neglected during the development of the datasets. Moreover, some sorption coefficients were theoretically derived from other predictive models using the octanol/water partition coefficient ( K OW ) (Wang et al., 2015; Daré et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Muhire et al., 2021; Cantwell et al., 2022; Jiang et al., 2022). Of course, these theoretical data lack a specific experimental methodology and soil/environmental properties that explain their values and guide their interpretation.

The previous findings impact the applicability of QSAR models in two ways: (i) empirical values were used together, independent of their quantification method (e.g., one-point vs. isotherm, batch vs. field-based) or sorption trend (e.g., sorption isotherm linearity); and (ii) datasets in the considered literature were built using mean or median sorption coefficient values such that the amount of data is equivalent to the number of pollutants and soil variability was part of the experimental error. Both issues increase the uncertainty when trying to assess the representational value of data for reliable predictions.

Furthermore, findings impact the reproducibility of QSAR models: (i) only three QSAR models provided the units of sorption coefficients (Wang et al., 2015; Rybacka and Andersson, 2016; Berthod et al., 2017), necessary to evaluate the comparability among empirical values obtained from different studies due to possible changes in the isotherm shape (linearity, 1 / n ); and (ii) more recent QSAR studies used data from previous articles that were in turn obtained from even older papers, based on well-known databases and the ability to contrast algorithms among QSAR models if they share the same dataset. As a result, several QSAR articles shared their datasets (Olguin et al., 2017; Olguin et al., 2019; Kobayashi and Yoshida, 2021; Daré et al., 2017; Muhire et al., 2021; Wang et al., 2015; Sabour and Moftakhari Anasori Movahed, 2017) or part of the data (Zhang et al., 2018; Pandey and Roy, 2021; Wang et al., 2015; Aranda et al., 2016; Kobayashi et al., 2020). However, we were unable to determine the empirical origin of data used in those articles. Moreover, the amount of data among QSAR articles sharing datasets differed without an explanation, which affects their interpretation and usability.

In contrast to QSAR, articles concerning PTF included the description of the site and procedure for soil sampling, with most of the sorbents being agricultural soils. Generally speaking, soil samples were superficial, from 0 to 5 (De Gerónimo et al., 2018), 20 (Conde-Cid et al., 2019; Conde-Cid et al., 2020) or 25 cm depth (Kodešová et al., 2015; Klement et al., 2018), with few cases addressing the spatial variability (Kodešová et al., 2015; Singh et al., 2016). In some cases, confounding variables were minimized using soils without application of the pollutants or fertilizers in previous years (De Gerónimo et al., 2018) or by evaluating the presence of pollutants before the sorption study (Conde-Cid et al., 2020).

Additionally, methodological and environmental conditions such as soil:solution ratio, interval of concentrations, background electrolyte, temperature, and contact and equilibrium time impacted the magnitude and interpretation of sorption coefficients and therefore were included during the data collection and treatment. For instance, the selection of the interval of concentrations or the soil:solution ratio impacted the isotherm shape (Dollinger et al., 2015), with several nonlinear sorption isotherms fitting to the Freundlich model (Dollinger et al., 2015; Kodešová et al., 2015; Sidoli et al., 2016; De Gerónimo et al., 2018; Klement et al., 2018; Conde-Cid et al., 2019; Conde-Cid et al., 2020; Hu et al., 2022), including the prediction of the linearity coefficient ( 1 / n ) in addition to the sorption coefficient (Dollinger et al., 2015).

All sorption coefficients used in PTF models were quantified experimentally. Most of the PTF studies produced their own data using the same methodological conditions (Kodešová et al., 2015; Sidoli et al., 2016; Singh et al., 2016; De Gerónimo et al., 2018; Klement et al., 2018; Conde-Cid et al., 2019; Conde-Cid et al., 2020). Only in two cases data were collected from other studies (Dollinger et al., 2015; Hu et al., 2022), where methodological differences among studies were included as predictor variables and analyzed during the mechanistic interpretation of predictive models. Additionally, data treatment included a homogenization step when (i) sorption was approximately linear, so K F and K d were used jointly (Dollinger et al., 2015; Conde-Cid et al., 2019) or (ii) isotherm linearity was highly variable, affecting the comparability and interpretation of K F values, where average linearity coefficients ( 1 / n ) were calculated for each pollutant and their K F values were recalculated (Kodešová et al., 2015; Sidoli et al., 2016; Klement et al., 2018).

Another homogenization technique involved the simplification of the sorption mechanisms by minimizing the variability (i) between pollutants, e.g., using pH-dependent pollutants at similar pH values among soils (Conde-Cid et al., 2019) or avoiding mixtures between ionic and non-ionic forms (Sidoli et al., 2016); or (ii) between sorbents, e.g., using unmodified soils or sediments with OC < 20% and comparable background electrolytes (Dollinger et al., 2015).

QSAR models were focused on a common sorption mechanism applicable through several soils, while PTF models considered local variability and specific sorption mechanisms, in accordance with our proposed assumptions (Figure 4). This affected the data collection and treatment. For QSAR models, data lack methodological and environmental context because they are assumed as independent of soils (e.g., methodological conditions such as soil:solution ratio could be perceived as soil variability and therefore irrelevant). On the other hand, PTF models implemented a highly detailed procedure, including (i) description of soils, (ii) impacts of methodological and environmental conditions, and (iii) classification and treatment of data in agreement with the previous steps.

Only molecular properties were used by QSAR models to predict sorption coefficients. Those predictor variables, called “descriptors” in QSAR models, were generally related to hydrophobicity, such as K OW (Cai et al., 2019; Muhire et al., 2021) or the number of C-F bonds in a molecule (Jiang et al., 2022). In general, the simplified molecular input line entry system (SMILES), international chemical identifier (InChIKey) and structural data file (SDF) were sufficient to represent the 2D and 3D molecules. Then, a molecular structure optimization method was applied (Wang et al., 2015; Zhu et al., 2017; Zhang et al., 2018; Cai et al., 2019; Jiang et al., 2022). Finally, descriptors were quantified implementing several software solutions such as Dragon (Wang et al., 2015; Sabour and Moftakhari Anasori Movahed, 2017), EPI Suite (Aranda et al., 2016; Olguin et al., 2017), MOLE db (Sabour and Moftakhari Anasori Movahed, 2017), MOE (Berthod et al., 2017), Mordred (Kobayashi et al., 2020; Kobayashi and Yoshida, 2021), Multiwfn (Jiang et al., 2022), Open Babel (Kobayashi et al., 2020), OPERA (Kobayashi and Yoshida, 2021), and PaDEL-Descriptor (Aranda et al., 2016; Pandey and Roy, 2021). These software offer thousands of constitutional, topological, geometrical, electronic, thermodynamic, quantum chemical, and other kinds of theoretical and semi-theoretical descriptors. Descriptors quantified by each software were complementary yet also overlapped, making it necessary to eliminate constant or correlated descriptors in a future step.

Despite the previously described procedure, (i) hydrophobicity was in some cases inadequate for representing the sorption of ionic species, affecting the statistical quality of QSAR models (Rybacka and Andersson, 2016; Berthod et al., 2017), and (ii) the use of molecular descriptors (hydrophobic, hydrogen-bonding and charge-related interactions) to address soil variability produced poor predictive ability (Berthod et al., 2017). Both cases support our proposed assumption (Figure 4) as a requisite to apply QSAR models.

Physicochemical properties of soils were used as predictor variables for PTF models, including environmental and methodological conditions that varied within the dataset. The most common properties were pH in different solutions, OC, soil texture, and cation exchange capacity (CEC). More properties were added to represent specific soil orders or land uses, such as variable charge (e.g., exchangeable aluminum, crystallized and amorphous oxy-hydroxides) (Sidoli et al., 2016; De Gerónimo et al., 2018; Conde-Cid et al., 2019), salinity (e.g., CaCO₃ content, hydrolytic and exchangeable acidity, base cation saturation) (Kodešová et al., 2015; Klement et al., 2018), or agricultural context (e.g., total organic N content, available P) (Sidoli et al., 2016; Conde-Cid et al., 2019). Some properties showed statistically significant correlations, such as OC with CEC and clay content with iron and aluminum oxides (De Gerónimo et al., 2018). Since these correlations depend on soils they cannot be generalized. Additionally, the quantification method varied among articles. In these cases, different methodological conditions within a dataset were (i) tested as descriptors, e.g., soil:solution ratio and maximum initial concentration of pollutants (Dollinger et al., 2015; Hu et al., 2022) or (ii) used as part of the data splitting (Dollinger et al., 2015).

According to our proposed assumptions (Figure 4), if the sorption mechanism is independent of soils, then it is not necessary to contextualize predictor variables or include soil descriptors. This explains why QSAR models generally follow an a posteriori approach, where a massive pool of molecular descriptors is used, and the mechanistic explanation is derived from the algorithmically selected descriptors to represent the outcome. On the other hand, PTF models explicitly address the complexity of the sorption process, following an a priori approach, where they explain the conceptual background of the sorption process and consequently propose a few predictor variables to develop predictive models.

Databases were split into a training and a test set, used to develop the QSAR model and assess its predictive performance, respectively. In this sense, data splitting helped to evaluate the reliability of QSAR models when applied to external data, providing an estimate of the performance for new pollutants. Ideally, data are uniformly distributed between training and test set. In the literature, the process was random (Wang et al., 2015; Berthod et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Muhire et al., 2021; Cantwell et al., 2022; Jiang et al., 2022) or rational, depending on the splitting algorithm, where Y-ranking is the most frequent approach (Zhu et al., 2017; Cai et al., 2019; Olguin et al., 2019; Kobayashi and Yoshida, 2021). Considered training:test set ratios include 14:86 (Aranda et al., 2016), 66:33 (Olguin et al., 2017; Olguin et al., 2019; Kobayashi and Yoshida, 2021; Cantwell et al., 2022), 70:30 (Sabour and Moftakhari Anasori Movahed, 2017), 75:25 (Wang et al., 2015; Daré et al., 2017; Muhire et al., 2021; Pandey and Roy, 2021), 80:20 (Berthod et al., 2017; Zhu et al., 2017; Cai et al., 2019; Pandey and Roy, 2021; Jiang et al., 2022), and 88:12 (Kobayashi et al., 2020). In one case, the effect of dataset size was also evaluated, showing that the division of the training set (N = 643) into eight subsets (N = 79 - 81) produced equivalent QSAR models to those for the whole dataset (Olguin et al., 2019).

Some articles applied mechanistic criteria for splitting, such as the development of different QSAR models based on the whole dataset and specific chemical classes (aliphatic groups, monoaromatic hydrocarbons, diphenyl ethers, polyaromatic hydrocarbons and plant protection products), with their respective training and test sets (Pandey and Roy, 2021), the use of different test sets based on intervals of low, medium and high sorption coefficient values (Daré et al., 2017), the use of several training and test sets to evaluate the stability of the mathematical equation and its predictability (Berthod et al., 2017), and the development of QSAR models based on the charge of the pollutants (neutral, positively, and negatively charged) (Rybacka and Andersson, 2016; Berthod et al., 2017).

Different algorithms were used to develop the QSAR models, most of them assuming a linear relationship between descriptors and outcome, such as multiple linear regression (MLR) (Wang et al., 2015; Aranda et al., 2016; Berthod et al., 2017; Kobayashi et al., 2020; Kobayashi and Yoshida, 2021; Pandey and Roy, 2021; Jiang et al., 2022), partial least squares regression (PLSR) (Rybacka and Andersson, 2016; Berthod et al., 2017; Daré et al., 2017; Zhu et al., 2017; Cai et al., 2019; Pandey and Roy, 2021), support vector machines (SVM) (Kobayashi et al., 2020; Kobayashi and Yoshida, 2021), and univariate linear regression (ULR) (Berthod et al., 2017; Olguin et al., 2017; Zhang et al., 2018; Olguin et al., 2019). Nonlinear models included gradient boosting decision tree (Kobayashi et al., 2020; Kobayashi and Yoshida, 2021) and neural network-based models (Berthod et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017) algorithms.

MLR was the most frequently used algorithm due to the simple mechanistic explanation derived from those QSAR models. However, previous treatments were required due to the large number of descriptors, e.g., stepwise selection (Wang et al., 2015; Berthod et al., 2017; Pandey and Roy, 2021; Jiang et al., 2022), best subset selection (Pandey and Roy, 2021) and replacement method (Aranda et al., 2016). Additionally, PLSR was useful for dimensionality reduction and to avoid multicollinearity (Zhu et al., 2017; Cai et al., 2019; Pandey and Roy, 2021). On the other hand, nonlinear regressions were able to represent complex relationships between sorption coefficients and descriptors. Finally, ULR was used to assess simple relationships between sorption coefficients and hydrophobicity.

Data splitting was always related to mechanistic issues, such as the use of different outcomes (Hu et al., 2022), the use of one pollutant per model (Kodešová et al., 2015; Sidoli et al., 2016; Singh et al., 2016; Klement et al., 2018; Conde-Cid et al., 2020) and/or per site or plot from which the data were obtained (Singh et al., 2016), and the methodological conditions and available information about specific physicochemical properties (Dollinger et al., 2015). However, only three articles included training and test sets (Singh et al., 2016; Conde-Cid et al., 2019; Hu et al., 2022), negatively affecting the validation of most of the PTF models.

Algorithms for developing PTF models included MLR (Dollinger et al., 2015; Kodešová et al., 2015; Sidoli et al., 2016; De Gerónimo et al., 2018; Klement et al., 2018; Conde-Cid et al., 2019; Conde-Cid et al., 2020; Hu et al., 2022) with stepwise selection (Dollinger et al., 2015; Sidoli et al., 2016; De Gerónimo et al., 2018), PLSR (Singh et al., 2016) and nonlinear algorithms (Sidoli et al., 2016). One study used principal component analysis (PCA) to guide the interpretation of descriptors for the MLR (De Gerónimo et al., 2018). Some MLR algorithms included the transformation of physicochemical soil properties that did not follow a normal distribution (Conde-Cid et al., 2020), or the use of exponential relationships between sorption coefficients and predictor variables to improve the statistical quality (from R² < 0.75 to > 0.92) but changing the mechanistic interpretation (Sidoli et al., 2016).

Data splitting preceded the generation of a mathematical equation. This had a statistical explanation in QSAR models (i.e., reliability by using training and test sets), with a few articles including characteristics of pollutants. The procedure and the use of several explicit algorithms for developing the QSAR models agrees with the OECD principle P2 (Figure 3). Inversely, PTF models lacked a statistically validated procedure but evaluated diversity and complexity of data and used the data splitting to better represent sorption mechanism, according to pollutant-soil interactions. This procedure is not included in QSAR models because the sorption mechanism was assumed as generalizable and independent of soils (Figure 4).

The validation of QSAR models was based on the OECD principle P4 (Wang et al., 2015; Sabour and Moftakhari Anasori Movahed, 2017; Kobayashi et al., 2020; Kobayashi and Yoshida, 2021), considering several statistical parameters and procedures classified in two steps: internal validation (training set) related to the goodness-of-fit and robustness, and external validation (test set) to evaluate the predictive performance.

Common statistical parameters for goodness-of-fit included R² (Berthod et al., 2017; Daré et al., 2017; Olguin et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Zhu et al., 2017; Zhang et al., 2018; Cai et al., 2019; Olguin et al., 2019; Kobayashi et al., 2020; Kobayashi and Yoshida, 2021; Muhire et al., 2021; Pandey and Roy, 2021; Cantwell et al., 2022) and adjusted R² (Wang et al., 2015; Muhire et al., 2021; Jiang et al., 2022), F-Test value (Zhang et al., 2018; Muhire et al., 2021; Cantwell et al., 2022), standard error of the estimate (Zhu et al., 2017; Zhang et al., 2018; Cai et al., 2019; Kobayashi et al., 2020; Cantwell et al., 2022), residual sum of squares (Olguin et al., 2017; Kobayashi et al., 2020), root mean square error (RMSE) (Wang et al., 2015; Daré et al., 2017; Olguin et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Olguin et al., 2019; Kobayashi and Yoshida, 2021; Muhire et al., 2021; Jiang et al., 2022), p-value of descriptors contained in the QSAR model (Zhu et al., 2017; Cai et al., 2019; Muhire et al., 2021; Jiang et al., 2022), variance inflation coefficient to evaluate multicollinearity (Jiang et al., 2022), mean absolute error (Pandey and Roy, 2021), and concordance correlation coefficient (CCC) (Olguin et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Olguin et al., 2019; Kobayashi and Yoshida, 2021).

Robustness was typically addressed through one of the following procedures: “leave-one-out” cross-validation (Wang et al., 2015; Aranda et al., 2016; Berthod et al., 2017; Daré et al., 2017; Olguin et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Zhang et al., 2018; Cai et al., 2019; Olguin et al., 2019; Kobayashi et al., 2020; Kobayashi and Yoshida, 2021; Muhire et al., 2021; Pandey and Roy, 2021; Cantwell et al., 2022; Jiang et al., 2022), “leave-many-out” cross-validation (Olguin et al., 2017; Olguin et al., 2019), and bootstrapping (Wang et al., 2015). Additionally, Y-scrambling was used to detect or discard chance correlations (Aranda et al., 2016; Daré et al., 2017; Olguin et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Zhang et al., 2018; Olguin et al., 2019; Kobayashi and Yoshida, 2021; Muhire et al., 2021; Pandey and Roy, 2021). Other statistical parameters such as RMSE and CCC for the internal validation were also used (Daré et al., 2017; Olguin et al., 2017; Zhang et al., 2018; Olguin et al., 2019; Kobayashi and Yoshida, 2021; Muhire et al., 2021).

External validation was evaluated through different indicators, such as variance explained in external prediction (Wang et al., 2015; Daré et al., 2017; Olguin et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Zhu et al., 2017; Zhang et al., 2018; Cai et al., 2019; Olguin et al., 2019; Kobayashi and Yoshida, 2021; Muhire et al., 2021; Pandey and Roy, 2021; Cantwell et al., 2022; Jiang et al., 2022), modified coefficient of determination of the external validation (Daré et al., 2017; Olguin et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Zhang et al., 2018; Olguin et al., 2019; Kobayashi et al., 2020; Muhire et al., 2021; Pandey and Roy, 2021), standard error of prediction (Zhu et al., 2017; Cai et al., 2019; Cantwell et al., 2022), RMSE of validation (Wang et al., 2015; Sabour and Moftakhari Anasori Movahed, 2017; Zhang et al., 2018; Olguin et al., 2019; Kobayashi and Yoshida, 2021; Muhire et al., 2021), external CCC (Zhang et al., 2018; Olguin et al., 2019; Kobayashi and Yoshida, 2021), and mean absolute error (Pandey and Roy, 2021).

Validation of PTF models was scarce and focused on goodness-of-fit, using R² or adjusted R² (Kodešová et al., 2015; Sidoli et al., 2016; Singh et al., 2016; Klement et al., 2018; Conde-Cid et al., 2019; Conde-Cid et al., 2020), p-value (Kodešová et al., 2015; Klement et al., 2018; Conde-Cid et al., 2019) and standard error (Singh et al., 2016), but not including an analysis of overfitting (e.g., observation:descriptor ratio). Two articles assessed the model robustness through cross-validation (Dollinger et al., 2015) and bootstrap methods (De Gerónimo et al., 2018). Articles that considered validation sets followed different approaches to validate their models: (i) developing the PTF model using the training set but quantifying R², Nash-Sutcliffe efficiency, RMSE, and absolute error only for the validation set (Hu et al., 2022), thereby causing the model to lack goodness-of-fit and robustness, (ii) plotting the measured vs. estimated sorption coefficient values for the training and test set and checking how many data fell within the interval of the measured value ±   2 units (Conde-Cid et al., 2019), or (iii) applying the same goodness-of-fit parameters for the validation set (R², standard error) (Singh et al., 2016).

The applicability domain (AD) is a theoretical chemical space where QSAR models make reliable predictions (Kobayashi et al., 2020; Kobayashi and Yoshida, 2021). This region is defined by the diversity of pollutants (molecular structures) in the training set and the descriptors that are used to predict their endpoints (Aranda et al., 2016; Zhu et al., 2017). The AD is specific to each QSAR model, and reliability is only possible to assess for molecules and properties that fall within this chemical space, based on similarity (Daré et al., 2017; Zhu et al., 2017). Otherwise, a prediction is an unreliable model extrapolation (Aranda et al., 2016; Zhu et al., 2017). Without an explicit definition of the AD, a predictive model does not meet the OECD principle P3 (Zhang et al., 2018).

AD was quantified by standardization (Wang et al., 2015; Daré et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Zhu et al., 2017; Zhang et al., 2018; Cai et al., 2019; Olguin et al., 2019; Kobayashi et al., 2020; Kobayashi and Yoshida, 2021; Pandey and Roy, 2021), leverage (Wang et al., 2015; Aranda et al., 2016; Daré et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Zhu et al., 2017; Zhang et al., 2018; Cai et al., 2019; Olguin et al., 2019; Muhire et al., 2021), Euclidean distance (Kobayashi et al., 2020; Kobayashi and Yoshida, 2021), and one-class support vector machines (Kobayashi and Yoshida, 2021). These methods were applied to detect outliers (usually defined as > 3.0σ from the mean) and/or influential points (leverage > threshold). AD was commonly visualized through a Williams plot considering the standardized cross-validated residuals versus leverage values of pollutants (Wang et al., 2015; Sabour and Moftakhari Anasori Movahed, 2017; Zhu et al., 2017; Zhang et al., 2018; Cai et al., 2019; Olguin et al., 2019; Muhire et al., 2021; Jiang et al., 2022). Here, outliers and influential points were molecular structures, and their analysis was based on molecular descriptors, delimiting the scope and interpretability of QSAR models (Cai et al., 2019; Jiang et al., 2022).

The considered PTF models did not include an AD. Therefore, PTF model applicability was not statistically validated.

QSAR models presented a clearly defined AD, guided by the OECD (P3, Figure 3). AD was independent of methodological or environmental conditions from which the empirical values are obtained, which fits with the QSAR assumption that we proposed (Figure 4) and impacts the OECD principle P5 (Figure 3). On the other hand, PTF models lack a defined AD but used soil physicochemical properties and methodological conditions as predictor variables, implying that their AD would be defined by the diversity of pollutant-soil interactions, whose interpretability depends on soils and local conditions.

In general, QSAR models were considered as successful for estimating sorption coefficients and were applied to environmental risk assessment (Berthod et al., 2017; Sabour and Moftakhari Anasori Movahed, 2017; Zhang et al., 2018; Cai et al., 2019; Kobayashi et al., 2020; Kobayashi and Yoshida, 2021; Muhire et al., 2021; Jiang et al., 2022), helping the design of new molecules with less environmental impact (Sabour and Moftakhari Anasori Movahed, 2017; Kobayashi et al., 2020; Muhire et al., 2021; Jiang et al., 2022), providing objective decisions (Kobayashi et al., 2020; Kobayashi and Yoshida, 2021) and assessing soil remediation and potential leaching (Sabour and Moftakhari Anasori Movahed, 2017).

PTF models were also considered successful and useable for environmental fate and risk assessment (Dollinger et al., 2015; Kodešová et al., 2015; Singh et al., 2016; Hu et al., 2022), and identification of vulnerable soils (Conde-Cid et al., 2020; Hu et al., 2022), helping the developing of mitigation strategies and management practices when necessary (Conde-Cid et al., 2020; Hu et al., 2022). However, the use of PTF models was generally suggested more cautiously than for QSARs. For instance, PTF models were proposed as screening methods for regional approaches if local soils were included (Singh et al., 2016). So, their validity depends on non-soil variables such as the pollutants and environmental conditions (Hu et al., 2022). Furthermore, their descriptors and the interval of initial concentration for calculating sorption coefficients should be verified with new datasets if used on different soils (Dollinger et al., 2015; Kodešová et al., 2015).

QSAR models can be used for regulatory purposes for any soil at any condition. However, QSAR reliability is restricted to the validity of their assumption (Figure 4), which in turn depends on the main sorption mechanism within the dataset. On the other hand, PTF models represent diverse sorption mechanisms at local scale. In this sense, QSAR and PTF models are not exclusive but complementary.

We identified three pesticide types in the literature: herbicides (N = 24, 63% of studied pesticides), insecticides (N = 10, 26%), and fungicides (N = 4, 11%), which reflects the trend in global usage (De et al., 2014). Among all articles, N = 24 (77%) included herbicides.

A large percentage of pesticides had no or low pH-dependence: N = 14 (37%) non-ionizable, N = 6 (16%) (very) weak acid, and N = 3 (8%) (very) weak base. Furthermore, acids were more common than bases, with only one article studying a strong base.

Heterocyclic compounds and amide derivatives were frequent among pesticides N = 10 (26%) and articles N = 13 (42%).

A heterogeneous distribution of targets, acid-base activities and chemical classes was found in the literature for articles and pesticides. We related trends to the occurrence of the most frequent pesticides: glyphosate, 2,4-D and diuron (Table 3; Supplementary Figure S2), which are strong acid herbicides and are present in 10 out of 31 articles (32%). These pesticides have been found frequently in water bodies in North America (glyphosate) and Australia (2,4-D and diuron) (Sharma et al., 2019). Furthermore, glyphosate and 2,4-D are commonly used in Argentina, North/Central America, and Africa (Sharma et al., 2019), and studied using PTF models (Dollinger et al., 2015; Sidoli et al., 2016; Singh et al., 2016; De Gerónimo et al., 2018).

Sources of variability contrasted in complexity: (i) soil variability involved several sorbents, treatments and/or topographic conditions. Samples had different physicochemical properties, so the interpretation of trends may involve multiple potentially valid explanations, including the presence of confounding variables. Then, these studies are useful to represent real scenarios of environmental concern, but their evaluation of sorption mechanisms is limited due to the inherent uncertainty of the experimental design. On the other hand, (ii) other sources of variability generally involved one sorbent, modifying one specific variable at controlled values (e.g., temperature, pH), giving simpler interpretations that are only locally valid but have potential to be extrapolated in future research.

Three different approaches were used to interpret results, based on their complexity and limitations: (i) conceptual comparisons without statistical tools when only two contrasting cases were studied, e.g., two soils, treatments, or pollutants; (ii) correlations between selected physicochemical properties of sorbents or pollutants (Khorram et al., 2018; Pose-Juan et al., 2018; Silva et al., 2018; Ben Salem et al., 2019; Góngora-Echeverría et al., 2019; Wang et al., 2020; Xu et al., 2020); and (iii) statistical tools applied to all possible variables with the purpose of reducing uncertainty when interpreting experimental studies. Common statistical tools were similar to those used in predictive models: PCA (Caceres-Jensen et al., 2020; García-Delgado et al., 2020; Meftaul et al., 2021; Pavão et al., 2022), regression models (García-Delgado et al., 2020; Cáceres-Jensen et al., 2021; Siek et al., 2021), cluster analysis (Caceres-Jensen et al., 2020), multivariate ANOVA analysis (Alfonso et al., 2017), and correlation matrices between (i) sorption coefficients and physicochemical properties (Mosquera-Vivas et al., 2018; Loffredo et al., 2019; Agbaogun and Fischer, 2020; Cáceres-Jensen et al., 2021; Chen et al., 2021; Siek et al., 2021), (ii) sorption coefficients and pairwise interactions to represent the effect of the interaction between physicochemical properties on sorption (Cáceres-Jensen et al., 2021), and (iii) physicochemical properties to explore multicollinearity and avoid biased interpretations (Agbaogun and Fischer, 2020; Cáceres-Jensen et al., 2021; Siek et al., 2021).

Soil variability, especially SV [soil], was more frequent in the literature than other sources of variability (Table 2; Figure 5). Additionally, differences in complexity produced different approaches to interpret results, where the use of statistical tools within studies and comparison among studies affect their reliability and extrapolation.

The sorption of neutral and non-ionizable pollutants occurs preferably on OC and could be represented by K OC , especially the hydrophobic sorption mechanism. However, other sorption mechanisms are also possible and K OC of non-ionizable and neutral pollutants may vary among soils, principally for aromatic and heterocyclic compounds, affecting the validity of QSAR assumption (Figure 4).

Sources of variability can improve the connection between representational value of data and predictive models when used to (i) explore the correlations C1 to C3 (SV [soil]), (ii) detect properties of pollutants and soils affecting sorption (OV [pollut], SV [soil]), and (iii) describe the impact of methodological and environmental conditions in the sorption process in specific contexts (e.g., climate, agricultural practices; OV [exp], SV [treat]). As a result, empirical findings may guide the selection of suitable approaches (QSAR or PTF), outcomes and predictor variables to avoid misrepresentations (e.g., unexplored or underrepresented but relevant sorption mechanisms).

We suggest having a clear objective according to the expected performance and applicability of predictive models, especially the scale and degree of specificity. For instance, variable charge soils (Figures 4A, D) are not common in the literature and their sorption behavior may be hidden among the most common trends (Figures 4B, C, E, F). However, variable charge soils possess a high agricultural productivity, which makes them relevant from a regulatory perspective, specifically for agriculture-based economies from emerging and developing countries (Caceres-Jensen et al., 2019). Then, a predictive model created from large datasets without a mechanistic treatment of data will have general applicability, making the conventional sorption mechanisms and predictive models less useful to predict and interpret specific and important scenarios (e.g., Figure 4D). In this sense, models that have different scope or address infrequent cases provide information that complements our understanding of the sorption process.

Predictive models should represent the variability among data in the simplest way possible, considering the available information under comparable conditions. As an example, changes in sorption coefficients of soil samples at different depths (SV [spatial]) or influenced by different treatments (SV [treat]) were generally explained by changes in the physicochemical properties of soils (e.g., OC content (endogenous + exogenous), kind of OC, pH, salinity) (Khorram et al., 2018; Pose-Juan et al., 2018; Silva et al., 2018; García-Delgado et al., 2020; Siek et al., 2021; Mosquera-Vivas et al., 2018; Sousa et al., 2018; Dos Santos et al., 2019; das Chagas et al., 2020; Marín-Benito et al., 2017) (see Table 2). Thus, they behaved as different (but conceptually connected) soils and could be included in the same dataset without requiring predictor variables accounting for depth or treatment information. Nevertheless, the extent of change in the soil properties varied among pesticides (Marín-Benito et al., 2017; Pose-Juan et al., 2018) and soils (Pavão et al., 2022; das Chagas et al., 2020), indicating that specific pesticide-sorbent and soil-amendment interactions were relevant when SV [treat]+OV [pollut] and SV [soil, treat] were addressed, respectively.

The data treatment should consider sources of variability together with statistical tools to decide whether the findings are reliable, especially if the same pollutant was addressed at different levels of complexity and uncertainty (e.g., SV [soil] versus OV [exp]). Those cases require a well-defined method to compare studies with different experimental designs to analyze the information and improve our comprehension of the environmental fate.

Sources of variability are helpful for exploring new predictor variables and improving the AD of models in specific scenarios. The high presence of SV [soil] in the literature may be used to improve the predictability of QSAR models by incorporating soil descriptors based on empirical trends when hydrophobic sorption is dominant but not unique, while SV [treat] and SV [spatial] are helpful for exploring the impact of physicochemical properties of soils on pesticide sorption in agricultural contexts.

We recommend explicitly assessing sources of variability within the data and doing so from the simplest to the most complex (e.g., starting with OV [exp], finishing with soil variability). Thus, the first (simple) findings act as a conceptual basis to contextualize and simplify the exploration of later (more complex) trends.

Scientific interpretation should guide the decision-making process (Figure 1) when selecting scientific information for regulatory purposes. For example, if the agricultural impact of sorption studies is required to propose, apply or evaluate an environmental policy for soil productivity, then it is necessary to consider SV [spatial] and SV [treat], i.e., models that explicitly included those sources of variability and their interpretation, both related to management practices and application of amendments in agricultural soils. As a result, QSAR and PTF models act as a bridge among scientific and regulatory dimensions, involving complex and diverse decisions to help environmental entities to promote and select strategies for applying adequate models in proper scenarios.

Correlations C1 to C3 represent the most probable scenarios we found in the literature, but other simplifications could be applied to non-organic soil components if they are relevant in a subgroup of pollutant-soil systems. To this end, we propose to assess (i) pollutant properties that seem relevant (acid-base activity, chemical class, etc.), (ii) a set of relevant soil components, and (iii) relevant functional groups in the selected soil component. For instance, an alternative sorption mechanism could consider correlations of sorption coefficients and anionic pollutants (alternative C1) when sorption occurs mainly in oxide minerals (alternative C2), specifically in aluminum oxides (alternative C3), which fits with Figure 4D if pollutant-soil interactions involving OC and non-oxide minerals are negligible or less relevant.

Let us suppose that non-hydrophobic interactions are relevant (e.g., sorption of glyphosate). In that case, OC and pH do not necessarily correlate with sorption coefficients. Moreover, correlations may be positive or negative depending on the pollutant-soil interaction. However, correlations C1 to C3 might be applicable in a subset of data, given specific experimental conditions (e.g., inclusion of ionizable compounds in their neutral form, quantification of average K OC values based in soils where sorption coefficients and OC correlates positively). Then, pollutants with non-hydrophobic interactions may be included within a diverse set of molecules for developing QSAR models.

From a QSAR perspective, sorption is independent of soils and therefore, QSAR models are assumed to be generalizable (see their assumption, Figure 4). In this case, the extrapolation of the subset of data used for developing the model as if they represent the entire dataset (or available information from the literature) produces a hasty generalization fallacy, with the consequent risk of bias when making decisions in soils whose properties were not considered in the subset of data.

The data selection process may help us to address this issue. For example, a wide interval of K OC values among soils has been observed for the same pollutant when sorption involves different mechanisms (Skeff et al., 2018; Caceres-Jensen et al., 2019; Cáceres-Jensen et al., 2021). If this occurs, data will be chosen from a subgroup that minimizes variability, reducing the diversity of soils where the predictive model is applicable.

Data selection should be described in terms of soil physicochemical properties instead of minimization of the standard deviation only. Otherwise, this practice may support the belief that sorption is strictly hydrophobic and unexplained variability is attributable to incorrect experimental values instead of other sorption mechanisms (Olguin et al., 2019). In this scenario, the explicit description of mechanistic limitations should help environmental agencies to use predictive models for regulatory purposes in a narrow but valid group of pollutant-soil systems, based on the methodological and environmental conditions used to develop the model.

In the literature, both soil and pollutant properties produced changes in sorption coefficients. Furthermore, most of the empirical studies were focused on variability among soils, with only a few addressing different pesticides (SV [soil] versus OV [pollut]) (see Table 2). However, no study contrasted both sources of variability.

We analyzed two studies with non-normalized sorption coefficient values at comparable conditions (R² ≥ 0.95, same units) addressing SV [soil] and OV [pollut] for >2 soils and >2 pesticides (Agbaogun and Fischer, 2020; García-Delgado et al., 2020).

We calculated five K F pest (phenylurea herbicides) and eighteen K F soil (alfisols, inceptisols and entisols) values based on findings in the literature (Agbaogun and Fischer, 2020). The obtained COV were 46 ± 6 % for K F pest and 52 ± 4 % for K F soil , with OV [pollut] slightly higher than SV [soil]. Additionally, this study found a correlation between sorption coefficients and OC (Agbaogun and Fischer, 2020), suggesting that normalize K F to OC (i.e., K OC ) reduces the soil variability. When we normalized K F (Agbaogun and Fischer, 2020), the new COV values were 32 ± 6 % for K OC pest and 53 ± 7 % for K OC soil , increasing the relevance of OV [pollut] versus SV [soil].

Another study quantified the sorption of four herbicides in two soils, four amendment materials, and the amended soils (García-Delgado et al., 2020). We quantified two K F pest per pesticide considering (i) the untreated and treated soils, and (ii) the isolated amendments. We also quantified fourteen K F soil (2 control soils + 4 × 2 amended soils +4 isolated amendments).

Considering soils (control + amended), the obtained COV were 47 ± 24 % and 87 ± 14 % for K F pest and K F soil , respectively. Additionally, COV values for isolated amendments (excluding sewage sludge, that showed a different behavior) were 22 ± 9 % and 56 ± 9 % for K F pest and K F soil , respectively. The relevance of OV [pollut] is explained by the different chemical classes of herbicides studied, while the low variability in SV [treat] is related to a common sorption mechanism on amendments (hydrophobic sorption on aliphatic and aromatic carbon) (García-Delgado et al., 2020).

The applicability and performance of QSAR models increase when (i) SV [soil] or COV value of K F pest are minimized (e.g., normalizing sorption coefficients to soil properties), and (ii) OV [pollut] or COV value of K F soil are maximized. Therefore, the variability among sorption coefficients is better explained by changes in the molecular structure of pollutants than by soil properties (OV [pollut] > SV [soil]), which supports the assumption of QSAR models (Figure 4). Otherwise, PTF (OV [pollut] < SV [soil]) or hybrid models (OV [pollut] ∼ SV [soil]) become relevant. From studies one and 2, OV [pollut] > SV [soil] when (i) sorption coefficients partially or completely follow the correlation C2 and are then normalized to the corresponding soil properties, or (ii) the sorbents share similar properties due to the treatments, making their sorption mechanisms more similar among sorbents.

A third option is that QSAR models represent unrealistic but useful overall trends using average values as outcomes. Consider the following two studies.

A study quantified the sorption of four herbicides in four Mexican soils (Alfonso et al., 2017). We quantified K F pest = 2 , 0 ± 0 , 3   µ g 1 − 1 / n m L 1 / n g − 1 for sulfotep considering all the soils, and K F soil = 2 , 8 ± 0 , 5   µ g 1 − 1 / n m L 1 / n g − 1 for a soil from Chablekal, using two structurally similar pesticides (sulfotep and dimethoate) (Alfonso et al., 2017). The low COV (13% and 16%, respectively) suggests that the use of average values based on similarities to reduce or even neglect the variability is valid for both pesticides and soils (all soils behave as one unique sorbent, while both pesticides behave as one unique pollutant). In this case, the structural differences between both pesticides were not enough to produce an important change on sorption. Thus, both pesticides present the same sorption coefficients when used in a QSAR model.

The opposite situation occurs when the structural variability produces important changes on sorption. From another study involving five pesticides and 18 soils (Agbaogun and Fischer, 2020), we obtained similar K F pest values for pesticide with similar structure, such as diuron and linuron ( 9 ± 3 and 11 ± 5 mg  kg − 1 mg   L − 1 − 1 / n , respectively), or monuron and isoproturon ( 4 ± 2 and 2 ± 1 mg  kg − 1 mg   L − 1 − 1 / n , respectively). Additionally, we calculated K F soil values and observed a variability from 4 ± 2 (Ibd soil, alfisol) to 16 ± 8 mg  kg − 1 mg   L − 1 − 1 / n (Uib soil, inceptisol). Interestingly, K F pest and K F soil values correlated positively with lipophilicity of pesticides (e.g., log ⁡ K OW ) and the OC content, respectively, probably due to hydrophobic sorption. This supports the findings from QSAR models, where the normalization of K F soil to OC makes sorption independent of soils, while K F pest is described exclusively by hydrophobic molecular descriptors.

Similarity among pollutants and soils within the dataset affects their variability and mechanistic interpretation. If sorption is generally represented by one overall sorption mechanism, then the use of average values reduces the variability. Moreover, the normalization of sorption coefficients to relevant soil properties (i.e., correlation C2) makes their variability dependent on pollutant properties (case 2). In this scenario, QSAR models have an interpretation and physicochemical meaning, despite their conceptual issues: they represent a sum of unknown trends with compensatory effects that hide specific sorption mechanisms and methodological differences, showing general trends that do not represent local pollutant-soil interactions but allow to approximate them around a probable value. This might explain the unclear impact of some molecular descriptors when complex sorption mechanisms appear (Rybacka and Andersson, 2016). Therefore, these QSAR models cannot predict sorption coefficients in local contexts and should be useful as exploratory analysis prior to the screening step in risk assessment.