Description of the Study Area and Available Datasets

The area of interest—the Zahorska Lowland—is located in central Europe, more specifically in the western Slovak Republic. It is bounded by the river Morava in the west, and the Little Carpathians mountain range forms a natural boundary in the east. Most of the Zahorska Lowland is composed of Neogene sediments of a marine origin and younger Quaternary sediments covering the surface of the plains, which are mainly represented by clayey sands, drift sands, and sandy clays.¹² The main pedogenetic factors of the lowlands are azonal, such as the accumulation activity of streams and soil-disrupting floods along with soil erosion. The most widespread soils in this area are chernozem, arenosol, and fluvial soils on the fluvial plains of river Morava. An intense accumulation process of organic soil matter can be observed in the chernozems; therefore, they are appropriate for a large number of different plant species.

On the contrary, arenosols are soils at an early stage of development, and they contain almost no continuous vegetation on their surface, due to which the organic matter content is very low. They are suitable for growing crops with lower demands such as rye. Fluvisols are periodically disrupted by floods. Their profile is often constantly loaded with new layers of flood sediment (sludge).¹²

Classification Systems used in this Study

Among the various soil classification systems based on the soil texture, the ones most used in the Slovak Republic are as follows:

A comparison of these two classification systems clearly shows various amounts of fractions in each of them and the discordant limits for each fraction. Namely, the fraction of clay particles in the USDA classification system ranges up to a value of 0.002 mm, while the Kopecky classification sets this limit at 0.01 mm. When solving various tasks for which it is necessary to have soil texture data in the USDA classification (eg, the mentioned example with the application of the Rosetta model), datasets using the Kopecky classification system are not compatible, which can be a problem. As we have already mentioned, there are many different texture classification systems in the world. Hence, similar situations could arise more often, and the methodology which the authors would like to propose in this article may be generally useful.

Description of the Datasets Used

Two datasets were used in this article:

PSD Models and Their Fitting

The nine parametric models of PSD involved in this study were evaluated and compared in order to derive the PSD functions. The models were developed by using an optimization procedure in order to choose the most suitable set of parameters of each model. The models chosen are the Fredlund models with three and four parameters (FR3, FR4), the Weibull model with three parameters (WB3), the Andersson model with four parameters (AND4), the van Genuchten model with three parameters (VG3), two Gompertz models (GP2, GP4), and a logarithmical (LG) and an exponential (EXP) model. These models were previously applied in various works, eg, in Zhao et al,¹⁵ for determining the hydrological properties of soils adjacent to dams constructed in China. The parts of the models were chosen according to the results obtained by Botula et al.⁷

In the FR3 Fredlund model, three parameters needed to be optimized¹⁶ ie, a, b, and c, and it uses the predefined parameters d_f = 0.001 mm and d_m (d_f = 0.0001 mm is a parameter related to the amount of fine particles in a soil and d_m is the diameter of the minimum allowable particle size). The FR4 model contains four parameters, ie, a, b, c, and d_f, which should be found by the optimization procedure and the predefined parameter d_m. The WB3 model with three parameters, ie, a, b, and c (and two predefined parameters: d_min = 0.002 mm and d_max= 2.0 mm), was previously used in Refs. 17 and 18 for the creation of PSD curves of several diverse soils. AND4 is a four-parameter model (parameters a, b, c, e) developed by Jauhiainen,¹⁹ based on an original theory of textural and water retention soil properties presented by Andersson.²⁰ The VG2 and VG3 models were proposed by Haverkamp and Parlange² on the basis of van Genuchten's original soil water retention curve model (developed in 1980). The VG3 model was subsequently developed from the first derived VG2 version with two parameters, when the mutual relation m = 1 − 1/n was used. VG3 considers both fitting parameters m and n as different and independent of each other and uses three additional parameters (a, b, and c). The GP2 and GP4 models are two forms of the Gompertz model with two and four parameters (a, b) and (a, b, c and e). Their curves present specific cases of a logistic curve, which is more general than the normal one. The equation of this curve constitutes an asymmetric closed form. Both models were previously used for extracting PSD curves by Silva et al²¹ for soils in Brazil and by Nemes⁸ for soils in Germany and the Netherlands. The soil particles expressed in these models follow the Gompertz distribution. Finally, the logarithmic model and EXP model contain two parameters—LG parameters a and b and EXP parameters c and b.

All the models listed were used for determining the PSD functions by the optimization method. The L-BFGS-B method used box constraints, which means that for each variable (the model's parameter), a lower bound and an upper bound are given. This method is a limited-memory modification of the quasi-Newton method.²² It was implemented in R language.²³ The purpose of this optimization was to predict the points of the grain curve by each model as closely as possible to the observed data by searching for the proper parameters of the model. The problem to be solved should be defined by the objective function, which in this article is proposed to have the following form:

Table 1.

PSD models used in this study.

Description and Tuning of the RF Model

The usual process for finding the best model for obtaining a proper theoretical parametric PSD function means applying more methods, eg, the models from Table 1, comparing their predictive ability with the help of the observed data and some statistical goodness-of-fit indicator, and then finally choosing the best performing model. However, there is usually no best parametric model that is superior under all circumstances. Various parametric models frequently show different degrees of precision in different soils and different environments, so the application of a single parametric model often leads to a functional relationship that could be more precise in one part of a textural domain but less suitable in other parts. One of the possible solutions of this problem is the application of the ensemble methodology, which uses the best features of various parametric models for achieving more general results from fitting the data to the actual values measured. Moreover, as was proved in this study, such a meta-model usually has the capability to fix systematic errors, if they are produced by the individual models (underestimation, overestimation, multiplicative error, etc.).

The goal of the ensemble methodology is to combine the outputs of several models in order to improve the generalizability/robustness that could be obtained from any of the constituent models. The nine parametric models described hereinbefore were used for obtaining the parametric PSD curves selected for the present study. The outputs of these parametric models are inputs for the ensemble model, namely, a data-driven model based on the RF algorithm, which means that a stacking type of ensemble was used. Stacking introduces the concept of (1) base models (PSD models) and (2) a meta-model, which computes the final results and replaces the averaging procedure used, eg, in bagging. In such a way, stacking tries to learn which base models are more reliable than others, using the mentioned meta-model (RFs in the solved task) to discover how best to combine the output of the base models to achieve the final results. The results of the base models are de facto new data for the learning problem, and an RF algorithm is used to solve this problem.

The RF algorithm consists of a set of regression trees (in this study, the authors are addressing a regression problem; it could also be populated with classification trees in the case of studying a classification problem). The resulting RF prediction is an average of the values of these many tree outputs, each one of which is grown on a bootstrap sample of the training data. The user chooses the number of trees that will be in the RF ensemble. A bootstrap sample means that each tree is trained using a sample obtained by randomly drawing N cases with replacements from the original dataset, where N is the number of variables in that dataset. With each of these bootstrapped training sets, a different tree is obtained. For the regression, the values predicted by each tree are averaged to obtain the RF prediction. More details and a more mathematically based explanation of the RF algorithm could be found in Breiman.¹¹

Optimization of the RF Model

Data-driven models must be optimized to obtain reliable results that are as precise as possible. The optimization of the model mainly means finding the optimal model parameters. An RF has three tunable parameters: ntree (the number of trees to grow), mtry (the number of variables randomly sampled as candidates at each tree split), and nodesize (the minimum size of the terminal nodes), which has the main effect on the final precision of the model. Two concepts are applied in this work for tuning an RF: grid search and repeated cross-validation.

The grid search is designed in the optimization process to choose the values for each parameter of the model by consecutively picking them from a grid of predefined values and then calculating with these parameters. The best combination of the parameters is chosen from that iteration in which the highest degree of precision of the model was achieved.

This precision is evaluated as an average value of more runs of the so-called cross-validation process. The so-called repeated cross-validation is used in the present article.¹⁸ In this process, the initial step consists of randomly dividing the training data into several approximately equal-sized datasets called folds. The training process uses all the folds except one as the inputs to the model, and the one unused fold is used as the validating data. This process runs as many times as the number of folds created with each combination of parameters. Repeated cross-validation means that the initial random splitting of the training data into folds is repeated more than once. This repetition is applied to obtain a more generalized evaluation of the model. The precision of the model in each iteration of the grid search is in fact the average value of the assessed statistic (eg, root mean square error (RMSE)) from all runs of the model with one set of parameters, eg, if there are two repetitions and five folds, the resulting statistic is the average value from 10 values.

This tuning concept involves two purposes in this study: (1) finding the best parameters of the RF algorithm and (2) estimating the precision of the proposed model, which could be expected for future data.

Fitting of the Parametric PSD Functions

When dealing with the transformation of a soil texture classification from the Kopecky system to the USDA classification, it is only necessary to extrapolate one point of the PSD, ie, to find the grain curve value for a soil particle diameter of 0.002 mm. Other USDA fractions could be derived from the Kopecky classification system by basic arithmetic operations. Nevertheless, from a more general point of view, when it is necessary to deal with other classification systems, and more points or different points of a grain curve are needed to be modeled, the methodology is the same as proposed in this article for this one point. Moreover, in previous works, when finding the percentages of unknown fractions in a soil sample, interpolation problems were mainly solved.^7,15 In this work, it is necessary to solve the extrapolation problem on the left tail of a PSD distribution, which is more complicated, so the methodology will surely be suitable while solving relatively easier interpolation problems.

Figure 1 shows the results of fitting various PSD functions to the soil texture data measured by the approach described in the “Materials and methods” section. As can be seen from the graphic representation on the left side of the correlation matrix and the correlation coefficients on its right side, the majority of models work quite well. This is especially apparent in the last column of the correlation matrix, where the correlation of the model with the measured data is evaluated. The results are similar to those that have been obtained in other works.^7,15 Additional statistical coefficients, which serve for an evaluation of the model's fitting, are available in Table 2. These coefficients assess the correspondence between measured and computed values. The statistics used are mean error (ME), mean absolute error (MAE), mean squared error (MSE), RMSE, the percent bias between simulated and observed values (PBIAS%), and the Pearson correlation coefficient (r).

Figure 1.

Correlation matrix of PSD models and measured data (in the clay column).

Table 2.

Goodness-of-fit measures for the PSD models.

Figure 2 shows an evaluation of the best models for various soils by their absolute errors. In every sample of the data (dataset A), the model with the lowest error, while predicting the grain curve point of 0.002 mm needed for translation to the USDA classification, is indicated by the color key. As can be seen, there is no single best model that could be preferred, either for a whole set of data samples or for samples of various soil types. This means that the proposed ensemble methodology, which is a combination of all the successful models, could be useful in this task.

Figure 2.

Evaluation of the PSD fittings for various USDA soil types.

As can be seen in Table 2, some models are not suitable for extrapolation problems regarding PSD fitting, so they were excluded from the final ensemble modeling. This inappropriateness for extrapolation is especially clear for the Weibull model. This is the case not only in the study presented but also could be considered as a general result. It can be explained by the following: Besides the empirical results, which were obtained by the computational experiment accomplished in this work and evaluated in Table 2, an important issue for this model is that parameter d_min should be set to represent the minimal diameters of the soil particles, which are assumed to be present in the sample. However, in the case of the left tail extrapolation, this is the unknown information for which one is searching, so it cannot be correctly set in advance. The second model that was excluded from the final modeling was the LG model, because negative values of the clay content were computed by this model for some samples. This means that the results of the seven models finally served as inputs to the ensemble modeling.

According to the statistical values in Table 2, evaluations of the results by various statistical coefficients differ, eg, the EXP model is evaluated as the best model by the correlation coefficient (r), but the FR4 model is evaluated as the best model by the RMSE. Various goodness-of-fit statistics evaluate different aspects of fitting, eg, the emphasis of one is more from a perspective of the variances and other statistics capture the bias better. For example, although it is possible to see in Figure 1 that the best correlation coefficient is for the EXP model, its prediction has multiplicative errors, which are not evident if one is only using an evaluation by the correlation coefficient (see Fig. 1, left side).

In this article, the authors propose a methodology inspired by the idea of ensemble learning, in which the RF algorithm is built on the top of the predictions computed by various parametric PSD models (they are inputs to the RF) and the optimal final result is obtained with this ensemble.

Fitting of the RF algorithm

While producing the final model for fitting the PSD, there are two basic tasks that are necessary to deal with: (1) to find the optimal model (eg, the optimal parameters of the RF algorithm that are suitable for the task to be solved) and (2) to evaluate the model's expected performance. The predicted values of the USDA clay fractions in the seven PSD models derived from the dataset A were used for the model's calibration or so-called training. This dataset is used because in the training phase of data-driven modeling, it is necessary to know not only the input data (Kopecky grain fractions) but also the measured outputs (USDA clay fractions), which are known in this case as mentioned in the descriptions of dataset.

The basic problem with the training dataset in this task is that it is relatively small (43 samples). The usual, so-called validation set approach, which involves randomly dividing the available set of samples into two parts, ie, a training set and a validation or hold-out set,¹⁸ is not appropriate to apply here. Instead of this method, the authors used the repeated cross-validation approach described in the “Materials and methods” section of this article. Through the know-how of the data mining community as expressed in various books and papers,¹⁸ fivefold cross-validation was principally used with two repetitions. The resulting RF model is based on the best parameters obtained from the cross-validation evaluations. The expected precision of the model is computed by using the computed and observed data from the folds held out in each iteration of the cross-validation.

The fitted RF model has the following parameters: 500 trees, four variables randomly sampled as candidates at each tree split, and terminal nodes with a minimum size of 5. As has already been stated in the “Materials and methods” section, the purpose of the cross-validation was not only to find these optimal parameters but also to evaluate the precision of the proposed model, which could be expected for future data. The precision using regression coefficient r was evaluated for the ensemble model of the soil texture transformation from the Kopecky to the USDA classification as 0.971 and using RMSE as 0.0343. When these values are compared with the results of the individual models of the ensemble model (see Table 2), it can be seen that no model achieved such a degree of precision; hence, the usefulness of the ensemble approach for this study was verified.

Some data-driven models, such as artificial neural networks, have often been criticized because of their black box character. It is true that these models are generally not meant to be descriptive and are usually not well-suited for inferences. Because of this, the authors used RF model in this work, as it not only generates very accurate estimates and is considered to be one of the most effective data-driven algorithms but also offers some information that helps understanding the modeled task. The RF model also included the ability to measure the importance or influence of each of its variable.¹¹ For each tree, the MSE on the out-of-bag portion of the data is recorded. Then the same is done after permuting each predictor variable. The differences between the two are then averaged over all the trees, and in such a way, the importance of the variable is obtained (eg, a decrease in accuracy after the permuting of the variables over all the trees is measured). The importance of each variable of the proposed ensemble model (eg, the predicted values by the 10 parametric PSD models) is scaled from 0 to 100 and displayed in Figure 3. With this evaluation, it is possible to see which model is more important and useful for the final prediction. The most important ones are the EXP model and the Gompertz model with two parameters. Although the authors would like to underscore the general usefulness of the proposed methodology, it has to be said that in the case of the other tasks, especially in the case of interpolation problems, other models may have a greater impact on the ensemble model. This result allows us to assume that in the context of ensemble modeling, the correlation coefficient is more important for the model selection than the other statistical measures that evaluate errors (such as MSE, RMSE, etc.), because these error coefficients are better, for instance, for both Fredlund models (see Table 2). This is true only in the context of ensemble modeling; otherwise, the Fredlund models should be chosen (when one is deciding only between the individual models for the final modeling). This could be explained by the better ability of the ensemble to repair systematic errors than the individual inaccuracies. Surprisingly, the Andersson model, which is not the most accurate one, plays quite an important role in the final ensemble. This is due to the fact that the Andersson model has a small degree of correlation with both of the best models (Fig. 2), so this means it is different. An efficient ensemble should be composed of predictors that are not only sufficiently accurate but also dissimilar, in the sense that the predicted errors occur in different regions of the input space.²⁴ Obviously, combining several identical models results in no gain in precision. From the evaluation in Figure 3, it can be seen that the ensemble mechanism is also capable of excluding a model if it is redundant (the Van Genuchten model in our case).

Figure 3.

Importance of individual PSD models in the final RF ensemble.

The authors would like to emphasize the following practical aspect about ensemble modeling. According to the so-called no free lunch theorem, it is never clear in advance which PSD model best suits a particular task. For this reason, it is usually necessary to try more models. On the basis of the results of this article, it could be said that when more PSD models are already fitted, instead of selecting and using only the best one, it is better to compose an ensemble prediction based on all of these already fitted PSD models (or on the basis of a subset of these models). Forming an ensemble usually brings an improvement in precision as has also been confirmed in the case study in this article, and ensemble prediction is relatively easy to accomplish when the fitted models for a particular task are already available.