A Stacked Machine Learning Algorithm for Multi-Step Ahead Prediction of Soil Moisture

Abstract

A trustworthy assessment of soil moisture content plays a significant role in irrigation planning and in controlling various natural disasters such as floods, landslides, and droughts. Various machine learning models (MLMs) have been used to increase the accuracy of soil moisture content prediction. The present investigation aims to apply MLMs with novel structures for the estimation of daily volumetric soil water content, based on the stacking of the multilayer perceptron (MLP), random forest (RF), and support vector regression (SVR). Two groups of input variables were considered: the first (Model A) consisted of various meteorological variables (i.e., daily precipitation, air temperature, humidity, and wind speed), and the second (Model B) included only daily precipitation. The stacked model (SM) had the best performance (R² = 0.962) in the prediction of daily volumetric soil water content for both categories of input variables when compared with the MLP (R² = 0.957), RF (R² = 0.956) and SVR (R² = 0.951) models. Overall, the SM, which, in general, allows the weaknesses of the individual basic algorithms to be overcome while still maintaining a limited number of parameters and short calculation times, can lead to more accurate predictions of soil water content than those provided by more commonly employed MLMs.

1. Introduction

Soil moisture is a variable that substantially affects the interactions between the Earth’s surface and the atmosphere, both in meteorological and climatic aspects [1]. It plays a fundamental role in rainfall–runoff processes [2], influencing the division of precipitation into surface runoff, subsurface flow, and infiltration. It also affects the transformation of incoming radiation fluxes to the soil into latent and sensible heat fluxes from the soil to the atmosphere. Soil moisture also strongly impacts the interaction between climate and vegetation in its multiple aspects, primarily the phenomenon of evapotranspiration. Moreover, soil moisture is a major discriminating factor in the type and condition of vegetation in a region. Variations in soil moisture can, therefore, have a massive impact on agriculture, forestry, and ecosystems.

Soil moisture measurement can be conducted by using in situ probes [3,4] or by remote sensing methods [5]. The significant impact on infiltration and runoff phenomena gives soil moisture prediction a key role in flood risk management [6] and landslide risk monitoring [7]. Furthermore, predicting soil moisture and its changes is essential for predicting the onset of drought and planning irrigation [8], as soil moisture is a critical limiting factor for crop growth.

Traditional soil moisture prediction techniques include empirical formulas, models based on soil water balance, models based on soil water dynamics, and autoregressive moving average models (ARMA). Empirical formulas are simple but often inaccurate tools, while models based on water balance or soil water dynamics require very detailed knowledge of the investigation site and boundary conditions. ARMA models need large-scale, high-quality datasets for training, which frequently are not available. Compared to these traditional methodologies, higher prediction accuracy and improved reliability in estimating uncertainty can be achieved by models based on artificial intelligence (AI) algorithms, which have found increasingly widespread use in the prediction of hydrological quantities over the past two decades [9,10,11,12,13,14,15,16,17]. However, as with ARMA models, AI models require large, high-quality datasets, as well as the optimization of a number of hyper-parameters.

A large number of studies on soil moisture estimation were carried out using various machine learning algorithms: support vector regression (SVR), artificial neural networks (ANNs), model tree (MT), multivariate adaptive regression spline (MARS), and adaptive neurofuzzy inference system (ANFIS) [18,19,20,21,22,23,24,25,26,27,28].

Elshorbagy and Parasuraman [18] employed two types of ANNs, i.e., multilayer perceptron (MLP) and the higher-order (HO-NN) types, to estimate soil moisture by accumulating field data at three subwatershed soil covers. They considered precipitation, air temperature, net solar radiation, and soil temperature at various depths for feeding MLP and HO-NN models. They found that the HO-NN model had better performance than MLP. Liu et al. [29] proposed a hybrid ANN—SVR architecture to estimate water content at a study site located in Chongqing, China. The authors noted that the hybrid model clearly outperformed the individual models. Additionally, Ahmad et al. [30] used SVR to assess soil moisture at 10 sites in the Lower Colorado River Basin. SVR models were trained using 5 years of data. The best results obtained were characterized by correlation coefficients between 0.34 and 0.77 between the soil moisture predicted with SVR and the value measured with remote sensing techniques, with a root mean square error (RMSE) of less than 2%. Furthermore, the authors made a comparison with the results obtained from models based on ANN and multiple linear regressions (MLR), showing that they were outperformed by SVR.

Si et al. [19] employed ANFIS, MLP, and the Bayesian regularization neural network (BRNN) in order to estimate soil moisture content at two different depths: 40 and 60 cm. They applied 900 data sets from field measurement in order to develop the AI models. From their results, it was found that ANFIS provided more accurate prediction of soil moisture than the BRNN and the MLP models. In addition, Zanetti et al. [20] employed the MLP model to assess soil moisture content while considering various properties of five types of soils, such as the apparent dielectric constant, clay and organic matter contents, bulk density and sand, and the silt content. They found that the MLP model with various combinations of input variables, such as organic matter combined with apparent dielectric constant, was particularly effective. Karandish and Simunek [31] evaluated the superiority of ANFIS and SVR with HYDRUS-2D for predicting time-dependent soil moisture content obtained using a physical model under various water stress circumstances over the maize-growing time period of 2010 and 2011. Later, Cui et al. [21] used MLP and MODIS data to estimate soil moisture, achieving an acceptable level of accuracy. In another study, Prasad et al. [23] developed an ensemble committee machine (CoM) learning model based on ANN (ANN-CoM) and utilized it to predict monthly soil moisture in the upper and lower layer of soil. From their study, statistical results indicated outperformance of the ANN-CoM model in comparison with those yielded by extreme learning machine (ELM), RF, and M5Tree.

Moreover, Prasad et al. [24] found superiority of ELM with ensemble empirical mode decomposition and the Boruta wrapper algorithm (EEMD-Boruta-ELM) over standalone MARS, ELM, and the EEMD-Boruta-MARS models for estimating weekly values of soil moisture content. Cai et al. [32] found that the deep learning NN (DLNN) provided a more accurate prediction of daily soil moisture based on various meteorological factors (e.g., daily precipitation, daily mean surface temperature, average wind speed, average relative humidity, average air pressure, and average temperature) than the MLP model at depths of 10 and 20 cm. Achieng [26] successfully used an SVR model by Gaussian kernel to simulate soil moisture content when compared with SVR models developed by polynomial and linear kernels, MLP, and the DLNN models. In recent years, Yuan et al. [27] reported a permissible level of accuracy when the generalized regression NN (GR-NN) was employed in order to estimate the regional surface soil moisture by means of satellite observations as input factors. Maroufpoor et al. [25] proposed a hybrid model based on the adaptive neurofuzzy inference system (ANFIS) and grey wolf optimization (GWO) algorithms, which was then compared with ANN, SVR, and standalone ANFIS. The input parameters of the model were the dielectric constant, bulk soil density, clay content, and organic matter of 1155 soil samples. The ANFIS-GWO model proved to be the most accurate, followed by the standalone ANFIS and SVR models. Adab et al. [33] used RF, SVR, ANN, and elastic network (EN) regression to estimate soil moisture from data obtained from Landsat 8 optical and thermal sensors, and knowledge of land use in a semi-arid region of Iran. The best results, characterized by a Nash–Sutcliffe efficiency value of 0.73, were obtained with the RF algorithm. In Heddam’s [28] study, four machine learning models (MLMs, i.e., MT, RF, MARS, and MLP-NN) were successfully employed to estimate soil moisture content while considering only hourly soil temperature as the input variable (obtained from two USGS stations) and compared with multivariate linear regression (MLR) technique.

Therefore, in the current literature, various MLMs have indicated promising performance in the estimation of soil moisture content for various conditions of soil physical properties. However, there is a shortage of models for predicting future soil water content (SWC), even in the short term, that are both simple, based on a few easily measurable input variables, and highly accurate. The main objective of this study is to propose a novel ensemble daily SWC prediction model obtained by stacking [34] three individual machine learning algorithms: MLP, RF, and SVR. These three standalone algorithms were chosen both because individually they showed good predictive capabilities, and because they have different structures and, thus, their combination can overcome the weaknesses of each algorithm. Furthermore, these three algorithms, compared with more complex algorithms such as deep learning, have the advantage that they depend on few parameters, facilitating training and optimization operations, and are characterized by significantly shorter calculation times. To the best of the authors’ knowledge, there are no applications of stacked algorithms for short-term prediction of SWC in the literature so far. The performance of the stacked model is compared with that of the individual algorithms considering two different scenarios of input variables. The proposed model is trained and tested with data obtained from a measurement site in East Anglia, UK. In addition, changes in model accuracy are statistically analyzed as the prediction horizon increases, while remaining within the scope of short-term forecasts.

2. Materials and Methods

2.1. Standalone Machine Learning Algorithms

In this research, MLP, RF, and SVR algorithms were used both individually and combined through stacking. Only a brief description of the machine learning algorithms used is given here. Readers interested in learning more about the algorithms may consult the relevant references. An MLP is a simple feedforward [35,36] ANN that can approximate any continuous function. An MLP consists of at least three layers of nodes: an input layer, at least one hidden layer, and an output layer (Figure 1a). The input layer includes the nodes that acquire the input data. Each node of the hidden layer processes the values of the previous layer using a weighted linear sum, followed by a non-linear activation function. The output layer receives the processed data from the last hidden layer and transforms them into the resulting values. The training of the algorithm is performed using the backpropagation technique. The neural networks employed in this study had only one hidden layer.

RF [37] is an ensemble prediction algorithm obtained by combining a set of individual regression trees in order to predict a single value of the target variable (Figure 1b). In each individual regression tree [38] it is possible to identify a root node, which comprises the training dataset, a number of internal nodes, which define the conditions on the input variables, and leaves, which represent the actual values assigned to the target variables. A tree regression model is developed by recursively dividing the input dataset into subsets, conducted in such a way as to minimize the internal node variance. A multivariable linear regression model provides predictions for each subset. Each tree grows from a different bootstrap of the training dataset. In addition, at each node, only a portion of the variables are randomly chosen with respect to which to split. The number of these variables is kept constant during the growth of the forest. A pruning process significantly reduces the risk of overfitting.

The idea behind the SVR (Figure 1c) algorithm [39] is to provide an approximation of the true value with a function that is as flat as possible, and that brings the error within a certain threshold, defined by an ε-value. A simple way to understand the SVR algorithm is to imagine a “tube” with an estimated function (hyperplane) as the center line and boundaries on both sides defined by ε. The goal of the algorithm is to minimize the error by identifying a function that places as many points of the training dataset as possible within the tube, while reducing the “slack”. The concept of slack variables is simple: for any value that falls outside ε, its deviation from the margin is denoted as ξ. When these deviations are to be tolerated, the algorithm tends to minimize them as well. Therefore, the deviations ξ are added to the objective function to be minimized in the constrained optimization problem into which the regression problem turns. The need to ensure a balance between the flatness of the regression function and the tolerated slacks is met by tuning a regularization parameter C. In SVR, regression is performed in a higher dimension. For this purpose, a function is required that maps the data points in a higher dimension. This function is defined as the kernel. In this study, the radial basis function (RBF) was chosen as the kernel K(x_i, x_j):

$$K(x_i,x_j)=\exp \left(-\gamma {‖x_i-x_j‖}^2\right),\gamma >0$$

(1)

where x_i, x_j are two input vectors.

2.2. Evaluation Criteria

Four different evaluation criteria were employed to assess the accuracy of the prediction models: coefficient of determination (R²), root mean square error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE). The R² coefficient is an estimation of goodness of fit, taking values in the range [0, 1]. The more accurate a model’s predictions are, the closer its R² will be to 1. It is defined as:

$${\mathrm{R}}^2=\left(1-\frac{{\displaystyle \sum _t{\left(f_t-y_t\right)}^2}}{{\displaystyle \sum _t{\left(y_a-y_t\right)}^2}}\right)$$

(2)

where f_t is the predicted value at time t, y_t is the measured value at time t, and y_a is the averaged value of the measured data.

The RMSE is the standard deviation of the prediction errors, the so-called residuals, which measure the distance of the experimental points from the regression line. In practice, the RMSE quantifies the dispersion of the data around the line of best fit. It is evaluated as:

$$\mathrm{RMSE}=\sqrt{\frac{{\displaystyle \sum _t{\left(f_t-y_t\right)}^2}}{N}}$$

(3)

in which N is the total number of predicted values in the time series.

The MAE estimates the average size of errors in the forecasts as a whole, without taking their direction into account:

$$\mathrm{MAE}=\frac{{\displaystyle \sum _t\left|f_t-y_t\right|}}{N}$$

(4)

The mean absolute percentage error (MAPE) evaluates the average of the absolute percentage errors of the prediction model. For the purpose of calculating MAPE, percentage errors are considered without taking the sign into account:

$$\mathrm{MAPE}=\frac{100}{N}{\displaystyle \sum _t\left|\frac{y_t-f_t}{y_t}\right|}$$

(5)

2.3. Stacked Model Development

Stacking [40,41] is an ensemble machine learning procedure that combines a number of classification or regression models through a metaclassifier. Stacking can exploit the capabilities of several well-performing models on a regression task in order to outperform standalone models in achieving predictions. The individual regression models are developed on the basis of the entire training data set, then a metaclassifier is applied on the basis of the outputs (meta-features) of the individual models. The elastic net (EN) algorithm was selected as the meta-classifier to develop the stacked prediction models. EN algorithm [42] is a combination of the two most commonly used regularized variants of linear regression: the least absolute shrinkage and selection operator (LASSO) method and the ridge method. The LASSO method selects the most influential variables by introducing an absolute penalty in the ordinary least squares (OLS) regression. ridge regularization also introduces a penalty in the OLS formulation by penalizing square weights instead of absolute weights. Thus, large weights are penalized significantly, and many small weights are distributed over the feature spectrum.

Two prediction models, differing in input variables, were developed in this study. Each model was developed in four variants, each based on one of the different ML algorithms introduced before, namely MLP, RF, SVR, and the combination by stacking of the previous ones (Figure 2). Model A includes the following exogenous input variables: cumulative daily precipitation (P), average daily air temperature T, average daily wind speed (WS), and average daily relative humidity (RH). On the other hand, Model B only includes cumulative daily precipitation P as an exogenous input. In addition, both models include lagged values of SWC as input variables.

The optimal number of lagged values of SWC, as well as the optimal values of the hyperparameters of the individual ML algorithms, were chosen by means of a grid search optimization procedure aimed at minimizing the RMSE of individual forecasting algorithms. A grid search is simply a heuristic search procedure through a user-specified subset of the hyperparameter space of a learning algorithm. A grid search algorithm must be driven by some performance metric, in this case RMSE.

It was found that, in the case study investigated, the optimal number of lagged values of SWC to be considered as input is 7. In addition, the optimized hyper-parameters of the forecast models are shown in Table 1. Therefore, based on the optimization process, the following input and output values can be indicated for the two forecast models:

Model A—input: SWC_t-6, SWC_t-5, …, SWC_t, P_t, T_t, WS_t, RH_t; output: SWC_t+1, SWC_t+2, SWC_t+3

Model B—input: SWC_t-6, SWC_t-5, …, SWC_t, P_t; output: SWC_t+1, SWC_t+2, SWC_t+3

where subscripts indicate the number of the day. The generic variable was normalized according to the equation:

$$x_{Ni}=\frac{x_i-x_{\min }}{x_{\max }-x_{\min }}$$

(6)

Model training is the core of the process of developing a prediction model, during which the weights and biases of an algorithm are optimized to reduce the loss function in the prediction interval. Applying a supervised learning technique, model training leads to a mathematical representation of the relationship between data features and a target label. The training of each model was carried out using 80% of the available time series, from June 2017 to June 2019. As this is a forecasting problem on a time series, preserving the temporal continuity of the data itself is essential and it is not possible to randomly select the training dataset. The models were then tested on the remaining 20% of the time series. This division allowed the most accurate results to be obtained.

Table 1. Main hyperparameters of the forecasting algorithms.

Algorithm	Hyperparameter	Value
MLP	Number of hidden layers	1
	Number of hidden neurons	5
	Activation function	Sigmoid
RF	Number of trees	100
SVR	Kernel function	RBF
	C	2
	ε	0.01
EN	α	0.3

Table 1. Main hyperparameters of the forecasting algorithms.

Algorithm	Hyperparameter	Value
MLP	Number of hidden layers	1
	Number of hidden neurons	5
	Activation function	Sigmoid
RF	Number of trees	100
SVR	Kernel function	RBF
	C	2
	ε	0.01
EN	α	0.3

2.4. Case Study

The data used in this study were provided by the COSMOS-UK network of the UK Centre for Ecology and Hydrology. Specifically, data were obtained from the COSMOS-UK site in Fincham (https://cosmos.ceh.ac.uk/data, accessed on 1 September 2022), East Anglia, UK (Figure 3). East Anglia’s climate is generally dry and mild. The region is the driest in the UK. In many areas, rainfall is less than 600 mm per year and is fairly evenly distributed throughout the year (Figure 4). The Fincham site is located in a large flat field planted with winter wheat, oilseed rape, and sugar beet in a 6-year rotation. The soil type is a chalky loam, a calcareous mineral soil. Similar to the other sites in the network, the Fincham site is equipped with an instrument that uses cosmic rays to measure soil moisture. More details on the measurement technique can be found in [43,44,45]. Experimental data are related to volumetric SWC (%) = (volume of water/volume of soil) × 100. The time series of daily hydrological variables of interest analyzed (soil water content, cumulative rainfall, average air temperature, average wind speed, average relative air humidity) include data collected from 22 June 2017 to 31 December 2019. Some essential information regarding the data used in the study is shown in Figure 4 and

Table 2. Essential time series characteristics of measured SWC and other climatic variables.

	SWC [%]	Air Temp. [°C]	Wind Speed [m/s]	Rel. Hum. [%]
Mean	24.18	11.06	3.28	80.12
Median	25.00	11.12	3.03	81.43
Max	34.20	27.36	8.52	99.62
Min	9.40	−4.82	0.67	53.36
St. Deviation	5.16	5.54	1.42	9.50
CV	0.21	0.50	0.43	0.12
1st Quartile	20.55	6.79	2.20	73.00
3rd Quartile	27.90	15.50	4.14	87.82
Skewness	−0.57	0.00	0.88	−0.31

. The SWC data sample has a low coefficient of variation (CV) and negative skewness (Table 2).

3. Results

Table 3. Model A evaluation metrics.

			MLP	RF	SVR	Stacked Model
Model A (Training)	1-day-ahead	R2	0.957	0.992	0.942	0.968
		RMSE	1.092	0.49	1.267	0.937
		MAE	0.816	0.356	0.911	0.694
		MAPE	3.36%	1.49%	3.73%	2.85%
	2-days-ahead	R2	0.940	0.985	0.912	0.953
		RMSE	1.285	0.663	1.569	1.137
		MAE	1.009	0.469	1.139	0.861
		MAPE	4.22%	1.94%	4.68%	3.56%
	3-days-ahead	R2	0.928	0.977	0.891	0.941
		RMSE	1.406	0.829	1.752	1.276
		MAE	1.101	0.571	1.266	0.959
		MAPE	4.66%	2.36%	5.24%	3.99%
Model A (Testing)	1-day-ahead	R2	0.957	0.956	0.951	0.962
		RMSE	0.924	0.985	0.996	0.877
		MAE	0.741	0.787	0.744	0.673
		MAPE	3.41%	3.62%	3.35%	3.05%
	2-days-ahead	R2	0.940	0.938	0.927	0.946
		RMSE	1.146	1.217	1.264	1.053
		MAE	0.942	0.990	0.945	0.821
		MAPE	4.40%	4.59%	4.27%	3.74%
	3-days-ahead	R2	0.921	0.929	0.911	0.935
		RMSE	1.355	1.360	1.442	1.169
		MAE	1.105	1.113	1.069	0.921
		MAPE	5.25%	5.22%	4.83%	4.22%

shows the values of the evaluation metrics for the prediction Model A with reference to the 1-day-ahead, 2-days-ahead and 3-days-ahead SWC. The table shows the metrics for both the training and testing phase, for each of the individual algorithms and for the stacked model (SM).

With reference to the 1-day-ahead forecast, in the testing phase, the three standalone algorithms showed roughly equivalent accuracies, with R² varying between 0.957 (MLP) and 0.951 (SVR), while MAPE varied between 3.35% (SVR) and 3.62% (RF). The SM outperformed all other forecasting algorithms, being characterized by a higher R² of 0.961 and smaller errors, with MAPE of 3.05%. It should be noted that the metrics values for the testing phase were absolutely comparable to those for the training phase. The only algorithm for which there was a perceptible difference between the two phases was RF.

The scatter plots of the predicted SWC values versus the measured values (Figure 5) show the excellent performance of all forecast models, with the points lying along the line of perfect agreement. With reference to the stacked model for the 1-day-ahead forecast, Figure 6a shows the time series of the predicted and measured SWC, while Figure 6b shows the relative error in the same time series. The relative error is defined as the absolute error in the forecast divided by the actual value of the SWC. The SM could accurately reproduce both SWC peak values and value fluctuations. Moreover, the relative error was almost always in the range −5%, +5%, and in a few cases approached ±10%.

Considering the 2-days-ahead forecasts, it can be seen that all variants of Model A underwent a very slight reduction in accuracy, but the forecasts were still very good. With regard to the SM metrics, for example, it can be observed that R² decreased from 0.962 to 0.946, RMSE increased from 0.877 to 1.053, MAE increased from 0.673 to 0.821, and MAPE increased from 3.05% to 3.74%. Again, the stacked model outperformed the standalone models. The decrease in accuracy for the standalone models was of the same order of magnitude as that of the SM, as shown by the metrics (Table 3). Even with regard to 3-days-ahead forecasts, all variants of Model A showed a further slight decrease in accuracy. Again, the three individual algorithms led to comparable results, while the SM outperformed them all, as proved by the higher R² value and lower RMSE, MAE, and MAPE values. This reduction in accuracy is also found when examining the time series of the relative error, with reference to 2- and 3-days-ahead forecasts. In the former case, the relative error most frequently falls in the ±5–10% range, in the latter case it often exceeds 10%, even reaching 15% or more.

Focusing on the evaluation metrics of Model B (

Table 4. Model B evaluation metrics for MLP, RF, SVR and stacked model.

			MLP	RF	SVR	Stacked Model
Model B (Training)	1-day-ahead	R2	0.946	0.990	0.934	0.965
		RMSE	1.222	0.533	1.365	0.989
		MAE	0.914	0.394	0.979	0.737
		MAPE	3.72%	1.62%	3.94%	3.02%
	2-days-ahead	R2	0.919	0.976	0.892	0.943
		RMSE	1.495	0.835	1.749	1.258
		MAE	1.161	0.586	1.274	0.964
		MAPE	4.77%	2.38%	5.13%	3.98%
	3-days-ahead	R2	0.900	0.960	0.863	0.925
		RMSE	1.658	1.073	1.989	1.441
		MAE	1.286	0.745	1.479	1.109
		MAPE	5.32%	3.01%	5.98%	4.62%
Model B (Testing)	1-day-ahead	R2	0.951	0.943	0.941	0.949
		RMSE	0.982	1.145	1.069	0.976
		MAE	0.745	0.937	0.810	0.751
		MAPE	3.42%	4.28%	3.64%	3.39%
	2-days-ahead	R2	0.928	0.916	0.907	0.924
		RMSE	1.249	1.456	1.381	1.224
		MAE	0.964	1.198	1.028	0.973
		MAPE	4.48%	5.53%	4.59%	4.45%
	3-days-ahead	R2	0.903	0.896	0.880	0.902
		RMSE	1.513	1.667	1.606	1.411
		MAE	1.185	1.381	1.193	1.144
		MAPE	5.56%	6.43%	5.32%	5.29%

), it can be seen that, with regard to 1-day-ahead forecasts, MLP (R² = 0.951, RMSE = 0.982, MAE = 0.745, and MAPE = 3.42%) led to better results in the testing phase than RF and SVR. The SM (R² = 0.949, RMSE = 0.976, MAE = 0.751, MAPE = 3.39%) led to results practically equivalent to those obtained with MLP. The ensemble model, in this case, did not lead to better results than the most accurate standalone algorithm. Furthermore, the predictions provided by Model B were slightly less accurate than the corresponding ones provided by Model A, with the exception of the MLP algorithm, for which negligible differences were observed.

Figure 7 shows the scatter plots of the predicted SWC values compared to the measured values for Model B. Again, the regular arrangement of the points along the line of perfect agreement can be seen, with small deviations.

Referring to the SM for the 1-day-ahead prediction, Figure 8a shows the time series of the predicted and measured SWC, while Figure 8b shows the relative error in the same time series, in the case of Model B. Again, the SM was able to accurately reproduce both the peak values of the SWC and the value fluctuations. Moreover, the relative error, although again almost always in the range of −5%, +5%, in some cases exceeded ±10%, even approaching 15%.

Focusing on the 2-days-ahead forecasts, it can be seen that, even for Model B, all variants suffered a reduction in accuracy. Furthermore, all variants underperformed compared to the corresponding variants of Model A. However, the forecasts were still satisfactory. MLP (R² = 0.928, RMSE = 1.249, MAE = 0.964, MAPE = 4.48%) and the SM (R² = 0.924, RMSE = 1.224, MAE = 0.973, MAPE = 4.45%) again led to the best results. Finally, 3-days-ahead forecasts showed a further reduction in accuracy. The SM provided the best results, and its metrics took the following values: R² = 0.902, RMSE = 1.411, MAE = 1.144, MAPE = 5.29%. The forecasts were still very good, even though all Model B variants slightly underperformed compared to the corresponding Model A variants.

The results obtained were compared with those provided by a more traditional autoregressive integrated moving average with exogenous variables (ARIMAX) model [46,47], which is effective in forecasting non-stationary time series and was adopted as a benchmark. The details of such a model are not described here in order not to overburden the discussion. For these details, please refer to the essential literature. Considering again 80% of the training data, including the same exogenous variables as in Model A, and using the two optimal values of 6 and 2 for the two parameters p (autoregressive order) and q (moving average order), the results shown in Figure 9 and

Table 5. ARIMAX model evaluation metrics.

ARIMAX
1-day-ahead	R2	0.944
	RMSE	1.046
	MAE	0.831
	MAPE	3.88%
2-days-ahead	R2	0.924
	RMSE	1.284
	MAE	0.997
	MAPE	5.23%
3-days-ahead	R2	0.910
	RMSE	1.542
	MAE	1.124
	MAPE	5.63%

are obtained.

It is evident that the ARIMAX model slightly underperformed in all MLMs considered. In particular, it appears to be less effective in accurately reproducing all SWC fluctuations. However, the results can still be considered valuable in this case as well.

4. Discussion

The results shown above demonstrate that both Model A and Model B are able to provide satisfactory predictions of short-term SWC. Model A proved to be slightly more accurate. The presence of air temperature, relative humidity and wind speed among the input data allows for the consideration of evapotranspiration, which depends on the aforementioned climatic variables and, in most cases, is the main outflow of moisture from the soil. However, even the availability of daily cumulative rainfall data as the only exogenous variable allowed for accurate short-term SWC forecasts in the case study. The local climatic conditions, characterized by low evapotranspiration values, probably had a significant influence on the results, giving Model B even greater practical utility than Model A, due to the smaller number of exogenous variables needed. Obviously, this result needs to be verified under different climatic conditions.

The SM generally outperformed the standalone models. In some cases, for Model B, it provided comparable performance to the most accurate individual algorithm. It seems that the SM performs significantly better than the individual models from which it is combined if the number of input variables is increased. This statement, however, needs further investigation.

Further insight into the accuracy of the different prediction models can be pursued by analyzing the violin plots in Figure 10, which show the relative error distributions of all variants of Model A and Model B, for the three forecast horizons considered. The same violin plots also include the corresponding box plots. The following insights can be deduced from these plots:

• In the case of Model A, only the SVR-based variant was characterized by an appreciable bias, whereas in the case of Model B, an appreciable bias could be found in both the MLP- and SVR-based variants.
• The distribution of the relative error in both models was asymmetrical in many cases.
• The error distribution tended to become flatter as the forecast horizon increased, and the IQR of the relative error expanded as the forecast horizon increased. The standard deviation of the residuals increased as the forecasting horizon increased. For example, with reference to Model A based on the SM, it was 0.874, 1.049, and 1.164% for the 1-day-ahead, 2-days-ahead, and 3-days-ahead forecasts, respectively. With reference to the SM-based Model B, the standard deviation of the residuals was 0.978, 1.227, and 1.414 for the 1-day-ahead, 2-days-ahead, and 3-days-ahead forecasts, respectively.
• The number of outliers resulting from forecasting models was very low.

This additional information provided by the violin plots enhanced the understanding of the results described above in terms of metrics. The improvement achieved through the stacked model is not only recognized in terms of metrics, but also in terms of error distribution, as can be seen from the violin graphs. The stacked model has a less asymmetric error distribution than the standalone models, in fact it is almost symmetrical. Furthermore, the relative error range is smaller, with a higher frequency of lower relative errors. This, for example, can certainly be an advantage in applications related to the optimization of irrigation planning under conditions of severe water scarcity, or in the assessment of landslide triggers.

The lack of benchmark datasets [48] and closely comparable studies prevents direct comparisons of the results. There are also very few studies focused on soil moisture that use stacking algorithms for purposes other than forecasting. A recent study by Das et al. [49] aimed to map soil surface moisture with a spatial resolution of 30 m in a semi-arid region using optical, thermal, and microwave remote sensing data, and applying machine learning techniques such as bagging, boosting, and stacking. The authors found that the stacking of the cubist, gradient boosting machine (GBM), and RF algorithms led to better results than the individual algorithms, in agreement with the findings of this study.

Other recent studies on SWC forecasting are based on the use of hybrid models. In terms of quantitative comparisons, the statistical measures provided by Models A and B showed that an improved MLM (i.e., the stacked model) outperformed MLP. This finding was evident in the comparison with Ahmad et al. [30] (R² = 0.2601 and 0.1764 for SVM and ANN, respectively). In the investigation by Ahmad et al. [30], the main limitations were that the input variables were obtained through satellite images, producing a high degree of uncertainty in the angle of incidence from the tropical rainfall measuring mission (TRMM), and in the normalized difference vegetation index (NDVI) from the advanced very high-resolution radiometer (AVHRR).

The performance of the ML models considered in this study is slightly better than that seen by Si et al. [19], who used ANN-Bayesian regularization (R² = 0.929) and ANN-Levenberg–Marquardt (ANN-LM) (R² = 0.932). It can be noted that the general structures of some ML models used here (i.e., RF, SVM, and stacked model) are more complex than those applied in the research of Si et al. [19].

Moreover, Cai et al. [32] provided soil moisture content predictions using deep neural network regression (DNNR) with as satisfying a degree of accuracy (R² = 0.98) as in the present research. Their success in the evaluation of soil moisture was due to considering a variety of input variables, such as average temperature, average pressure, relative humidity, wind speed, land temperature, daily precipitation, and initial soil moisture.

Furthermore, the performance of the present ML models was slightly better than that obtained by Heddam’s [28] investigation (R² = 0.925, 0.929, and 0.931 for M5MTree, MARS, and RF, respectively). In addition, the MLP-based model by Heddam [28] had rather lower accuracy results (R² = 0.885) than those reported in the present research for both Model A and Model B. Heddam [28] did not refer to the climatic variables that were considered in the present research. In fact, he used the soil temperature, the year number, the month number, and the day number in order to estimate the soil moisture content. His study indicated that climatic variables play a key role in improving the accuracy levels of ML models.

The main limitation of this study is that it considers only one case study. Therefore, the possible influence of different climatic conditions on the forecast models is not taken into account here. It will be interesting, in future developments of this study, to address the prediction problem under climatic conditions characterized by intense evapotranspiration and periods of widely varying rainfall (e.g., tropical climates). It will also be interesting to compare the results provided by the stacked model with those provided by models based on deep learning algorithms, which are known to perform very well in predicting time series [50,51]. Finally, the most ambitious goals will be pursued, such as developing models with a more distant forecasting horizon and models dependent only on exogenous climate variables.

5. Conclusions

This study introduced a novel forecast algorithm of daily volumetric soil water content, based on the stacking of the multilayer perceptron, random forest, and support vector algorithms. Two different input variable scenarios were considered, in order to develop two forecast models: Model A, which included daily precipitation, air temperature and humidity, and wind speed as exogenous variables, and Model B, which instead included only daily precipitation as an exogenous variable.

Both models provided very accurate predictions, with the coefficient of determination R² greater than 0.9 and MAPE not exceeding 5% in almost all cases, and with Model A generally outperforming Model B. In addition, for both models, the stacked algorithm-based variant generally outperformed the standalone algorithms. Both models experienced a modest reduction in accuracy as the forecast horizon increased, remaining within the range of short-term forecasts. In any case, even a model that only requires precipitation as an exogenous input variable is capable of providing adequate predictions for practical applications.

The proposed stacked model is simple, based on few parameters, very accurate, and has a very limited computational time. In the context of current research, which shows a marked tendency towards increasingly complex models, the proposed model can be considered an effective tool for facilitating the planning of irrigation activities and supporting flood risk management [52].