Several research studies have proposed RNN models to forecast soil movements and determine various triggering parameters for such movements (Xu and Niu, 2018; Xing et al., 2019; Yang et al., 2019; Jiang et al., 2020; Liu et al., 2020; Meng et al., 2020; Niu et al., 2021). For example, Yang et al. (2019) created an LSTM model for forecasting soil movements in China’s Three Gorges Reservoir area. The model was trained using reservoir water level, rainfall, and soil movement data. In this experiment, the support vector machine (SVM) model was also developed to compare with the LSTM model to forecast soil movements. From the results, it was found that LSTMs could effectively forecast soil movements. Similarly, Xu and Niu (2018) developed LSTM models to forecast the Baijiabao landslide’s displacement in China. The developed model was compared with a SVR and a backpropagation neural network. The LSTM model performed better than the SVR and neural network models. Furthermore, Xing et al. (2020) created LSTM and SVR models to forecast soil movements of the Baishuihe landslide in China. The findings indicated that the LSTM model could be used to forecast soil movements. Next, Liu et al. (2020) created LSTM, gated recurrent unit (GRU), and random forest (RF) models to forecast soil movements. These models were trained on the data recorded from the Three Gorges Reservoir, China. Results showed that GRU and LSTM models performed better compared to the RF model for forecasting soil movements. Furthermore, Meng et al. (2020) created an LSTM model and trained this model on the data collected from the Baishuihe landslide in the China. The recorded data from this landslide included different parameters like weather, rainfall, and soil movements. The univariate and multivariate versions of LSTM models were created on this dataset. The results revealed that the multivariate LSTM model performed better without overfitting. Similarly, Xing et al. (2019) developed a stacked LSTM model, where the sequence of soil movements were split into different subsequences. Next, the model used these subsequences to forecast soil movements. Niu et al. (2021) created an ensemble of EEMD and RNN models (the EEMD-RNN model). The proposed EEMD-RNN model was evaluated and compared to standard RNN, GRU, and simple LSTM models. The results showed that the EEMD-RNN model outperformed the individual RNN, GRU, and LSTM models. Next, Jiang et al. (2020) developed an ensemble of simple LSTM and SVR models to forecast soil movements. The Shengjibao landslide in the Three Gorges Reservoir area in China was taken as a case study. The results showed that the ensemble model outperformed the individual simple LSTM and SVR models.

In addition to the machine-learning literature related to landslides, several ensembles of RNN models have been tried for general sequence forecast and specific landslide susceptibility forecasting problems (Huang et al., 2015; Qiu et al., 2018; Zhang et al., 2019; Barzegar et al., 2020; Cui et al., 2020; Singh et al., 2020; Wang et al., 2020; Behera et al., 2021b). However, ensemble RNN approaches for forecasting soil movements have yet to be developed. For example, Huang et al. (2015) created a Bi-LSTM model for sequence labeling in NLP. According to the study, the Bi-LSTM model could learn past and future input features in a sequence (Huang et al., 2015). Next, Qiu et al. (2018) built the DGeoSegmenter using a Bi-LSTM model, which derived words randomly and combined them into phrases. Similarly, Cui et al. (2020) created a hybrid model of a temporal Bi-LSTM with a semantic gate, namely SG-BiTLSTM. The LSTM model was also developed to compare with the proposed model to identifying the landslide hazard-affected bodies (i.e., roads and buildings) in the images (Cui et al., 2020). Results revealed that the SG-BiTLSTM model was better than the LSTM model to classify the landslide affected bodies and extract features of images. Furthermore, Behera et al. (2021b) developed the Conv-LSTM to analyze consumer reviews posted on social media. An ensemble of CNN networks and a simple LSTM model was created for the sentiment classification of reviews posted across diverse domains. The experimental results showed that the Conv-LSTM outperformed other machine learning approaches in accuracy and other parameters. Furthermore, Singh et al. (2020) also created a Conv-LSTM model for crowd monitoring in large-scale public events. In this research, five different LSTM models were also developed to compare the performance of the Conv-LSTM model (Singh et al., 2020). Results showed that the Conv-LSTM performed best amongst other models. Besides, Barzegar et al. (2020) created an ensemble CNN-LSTM model to forecast the water quality variables. In this experiment, the CNN-LSTM model was developed to compare with other ML models, namely, LSTM, CNN, SVR, and decision trees (DTs). Results revealed that the developed ensemble model performed better than the non-ensemble models (CNN, LSTM, DTs, and SVR). Similarly, Zhang et al. (2019) also developed an ensemble CNN-LSTM model for the zonation of landslide hazards. The developed model was compared to other shallow ML models in this experiment. This experiment showed that CNN-LSTM was better than other shallow ML models.

In addition, the development of RNN ensemble techniques have shown some scope in the social network analysis, NLP, and similarity measures predictions (Lin et al., 2019; Behera et al., 2020; Behera et al., 2021a; Behera et al., 2019). For example, an ensemble of the RNN models could be used to predict the critical nodes in extensive networks (Behera et al., 2020). Furthermore, Behera et al. (2021a) predicted the missing link using the similarity measures, where similarity measures calculated the similarity between two links. Similarly, Behera et al. (2019) used the similarity measures for community detection in large-scale networks based on similarity measures. An ensemble of RNN models could predict these similarity measures (Lin et al., 2019).

We found that an ensemble of different RNN models has not been explored in the past for soil movement forecasting. However, an ensemble of the RNN model with EEMD or SVR models have been developed (Xing et al., 2019; Jiang et al., 2020; Niu et al., 2021), or individual ML models like LSTM, GRU, SVR, SVM, DT, and RF (Xu and Niu, 2018; Yang et al., 2019; Liu et al., 2020; Meng et al., 2020) have been developed for soil movement forecasting. Furthermore, these developed ensembles and individual RNN models have been compared with individual ML models (LSTM, GRU, SVR, SVM, DT, and RF) to forecast soil movements. However, different variants of an ensemble of RNN models have not to be compared in the past. Moreover, it was found that some RNN models like the CNN-LSTM, Conv-LSTM, stacked LSTM, and the Bi-LSTM performed well in social network analysis and NLP problems. However, these models have not yet been developed for soil movement forecasting.

Overall, this paper’s primary objective is to fill the literature gaps highlighted above by introducing a new RNN ensemble model called “Bidirectional-Stacked-LSTM” or “BS-LSTM” for forecasting soil movements at the Tangni and Kumarhatti sites in the Himalayan region of India. Furthermore, we develop CNN-LSTM, Conv-LSTM, Bi-LSTM, and stacked LSTM models and compare the performance of these models with the BS-LSTM to forecast soil movements. To the best of the authors’ knowledge, this type of study of soil movement forecasting has never been executed in the Chamoli and Solan districts in India. Thus, the main novelty in this work is to develop an ensemble of the RNN models for soil movement forecasting at new sites in the Himalayan region.

The data for training RNN models were collected from the two landslide sites at Tangni and Kumarhatti (see Figure 1A). The Tangni landslide is located in the Chamoli district, India. The landslide is at longitude 79° 27′ 26.3″ E and latitude 30° 27′ 54.3″ N, and at an elevation of 1,450 m (Figures 1A,B). As depicted in Figure 1B, the study area is located on the National Highway 7, which connects Fazilka in Punjab with Mana Pass. This landslide is categorized into a rock-cum-debris slide (THDC, 2009). Furthermore, the Tangni site’s incline is 30° up the road level and 42° beneath the road level. The surrounding area of the landslide is a forest consists oak and pinewood trees. Frequent occurrence of landslides was recorded in this area in 2013, which caused financial losses to the travel industry (IndiaNews, 2013). To monitor soil movements, inclinometer sensors were installed at Tangni site between 2012 and 2014.

The Kumarhatti site is located in the Solan district, India, along the Kalka - Shimla railway track. The site is at longitude 77° 02′ 50.0″ E and latitude 30° 53′ 37.0″ N at an elevation of 1734 m (Figures 1C,D). The landslide debris often damaged the Kalka - Shimla railway line, which was recorded at the Kumarhatti site by the railway department (Surya, 2011; Chand, 2014). A low cost landslide monitoring system was setup in 2020 at the Kumarhatti site to detect soil movements (Pathania et al., 2020; see Figure 1C).

Soil movement data were collected daily between July 1st, 2012 and July 1st, 2014, from the sensors installed at the Tangni site (see Figure 1B). Soil movement data (in meter) was recorded from the accelerometer sensors installed at the Kumarhatti site between September 10th, 2020, and June 17th, 2021. The accelerometer sensors were installed at 1 m depth from the hill surface at the Kumarhatti site.

The Tangni and Kumarhatti datasets included four distinct attributes: the timestamp, the borehole number, the sensor depth, and the soil movement. For example, for Tangni dataset, [2, 2, 12, 1.8°] denotes that the sensor in borehole two at a depth of 12 m recorded the soil movement of 1.8° in the upward direction during the second week since deployment. For Kumarhatti dataset, [20, 1, 1, 0.001 m] denoted that the accelerometer at 1-m depth in borehole one showed the downward movement of 0.001 m in the 200th minute since deployment.

Table 1 displays the mean and standard deviation (SD) in the Tangni and Kumarhatti datasets. For example, soil movements in borehole one had an SD of 1.83°. Furthermore, in Table 1, columns three through seven show the correlation value (r) of soil movements between different boreholes at Tangni. The correlation value between boreholes shows that if two boreholes are on the same failure plane of the landslide, both will simultaneously show soil movements with a high correlation. For example, borehole two was highly correlated with borehole four. As shown in Figure 1B, both boreholes are nearby, and it could be possible that both were on the same failure plane. The mean value of the soil movement for the Kumarhatti data was 0.02 m and had a standard deviation of 0.02 m. In Table 1, the borehole is denoted by BH.

The current values within the time series data may correlate to their previous values, which we call lags or look-back periods. The autocorrelation function (ACF) determines the number of lags present in a time series data. The partial autocorrelation function (PACF) determines a direct correlation between the current value and its lag values, removing the correlation of all other intermediate lags. As a result, the ACF value determines how many past forecasted errors are required to forecast the current value, and the PACF value determines how many past forecasted values are required to forecast the current value for the time series data. Figure 4A-4 L show the ACF and PACF plots with a 95% confidence interval for the Tangni and Kumarhatti datasets. As can from these figures, we can see that for the Tangni dataset, the time series one had an ACF value of five and a PACF value of one. The time series two had ACF and PACF values of zero, respectively. The time series three and four had ACF and PACF values of two, respectively. Furthermore, time series five had an ACF value of four and a PACF value of one. Similarly, the time series from the Kumarhatti dataset had an ACF value of forty-two and a PACF value of one. Based upon the range of ACF and PACF values across the two datasets, the look-back period for the developed models was varied from one to five for the Tangni data and one to forty-two for the Kumarhatti data.

Recurrent neural networks (RNNs) are specially designed to discover dependencies between current and previous values in time series data (Medsker and Jain, 1999). RNN networks are composed of a chain of cells linked by a feedback loop, and these cells extract temporal information from time series data. In general, every cell in RNNs has a simple design, such as a tanh function (Medsker and Jain, 1999). In general, the RNN models suffer from the exploding and vanishing gradients problem during the training process (Bengio et al., 1994). The problem arises when the long sequence of small or large values multiplies while calculating the gradient in the backpropagation. The exploding and vanishing gradient problems in the RNN model prevent the model from learning the long-term dependency in the data.

LSTMs solve the exploding and vanishing gradient problems by employing a novel additive gradient structure in the backpropagation. The additive gradient structure includes direct access to the forget gate activations, allowing the network to update its parameters so that the gradient does not explode or vanish (Hochreiter and Schmidhuber, 1997). Thus, LSTMs solve the problem of vanishing gradient and the learning of longer-term dependencies in RNNs models.

The simple LSTM is a type of RNN model that can remember values from previous stages (Medsker and Jain, 1999). The cell state in an LSTM acts as a conveyor belt ( c < t > ) , allowing unaltered information to flow through the units with only a few linear interactions. The internal architecture of each LSTM unit has an input gate ( i < t > ) , forget gate ( f < t > ) , and output gate ( o < t > ) . These three gates control the flow of information and avoid the exploding or vanishing gradient problems during training (Hochreiter and Schmidhuber, 1997) (see Figure 5A). The input gate adds new information from the new input and previous output to the cell state. The forget gate determines what information is retained for a long time and what information is removed from the cell state. The forget gate uses the sigmoid (logistic) function, where the sigmoid function’s output value is between zero and one. The forgot gate’s output zero means to remove the information from the cell state, and output one means to keep the information. The purpose of the output gate is to determine what output value its required from the cell state, and the output gate also updates the previous state of the hidden state ( h < t − 1 > ) . There is also a layer in the LSTM unit that contains the tanh activation function, which is used to update the state of neurons (see Figure 5A). Eqs 1–5 are the fundamental equations of the LSTM cell, with the Hadamard product denoted by the letter ′ ο ′: i < t > = σ ( W x i x t + W h i   h < t − 1 > + W c i ο   c < t − 1 > + b i ) (1) f < t > = σ ( W x f x t + W h f   h < t − 1 > + W c f ο   c < t − 1 > + b f ) (2) c < t > = i < t > ο   c < t − 1 > + i < t > ο   t a n h ( W x c x t + W h c   h < t − 1 > + b c ) (3) o < t > = σ ( W x o x t + W h o   h < t − 1 > + W c o ο     c < t − 1 > + b o ) (4) h < t > = o < t > ο   t a n h ( c < t > ) (5) Where, i < t > is the input gate at timestamp t; f < t > is the forgot gate at timestamp t; c < t > is the cell state at timestamp t; o < t > is the output gate at timestamp t; and h < t > is the hidden state at timestamp t. The variable x _t in the equations represents the input data sequence at timestamp t. The matrices W x i , W c i ,   W x f ,   W h f ,   W c f ,   W x c ,   W h c ,   W x o ,   W h o , and W c o are the weight matrices between two different layers. For example, the W x i is the weight matrix between input data sequence x and input gate i < t > . Similarly, the b i , b f ,   b o ,   and b c are the biases for the input gate, forgot gate, output gate, and cell state, respectively. The t a n h and σ here represent the tan hyperbolic and sigmoid (logistic) activation functions.

Conv-LSTM is a combination of convolution operation of the CNN model and the LSTM model (Shi et al., 2015). As shown in Figure 5B, the convolution operation is applied to the input and the hidden state of the LSTM cells. As a result, at each gate of the LSTM cell, the internal matrix multiplication process is changed by the convolution operation (*). This operation can find the underlying spatial information in high-dimensional data. The Conv-LSTM’s critical equations are provided in Eqs 6–10 below, where the convolution operation is denoted by ′*′ and the Hadamard product is denoted by ′ ο ′: i < t > = σ ( W x i * x t + W h i   * h < t − 1 > + W c i ο   c < t − 1 > + b i ) (6) f < t > = σ ( W x f * x t + W h f   * h < t − 1 > + W c f ο   c < t − 1 > + b f ) (7) c < t > = i < t > ο   c < t − 1 > + i < t > ο   t a n h ( W x c * x t + W h c * h < t − 1 > + b c ) (8) o < t > = σ ( W x o * x t + W h o * h < t − 1 > + W c o ο   c < t − 1 > + b o ) (9) h < t > = o < t > ο   t a n h ( c < t > ) (10) Where, i < t > is the input gate at timestamp t; f < t > is the forgot gate at timestamp t; c < t > is the cell state at timestamp t; o < t > is the output gate at timestamp t; and h < t > is the hidden state at timestamp t. The variable x _t in the equations represents the input data sequence at timestamp t. The matrices W x i , W c i ,   W x f ,   W h f ,   W c f ,   W x c ,   W h c ,   W x o ,   W h o , and W c o are the weight matrices between two different layers. Similarly, the b i , b f ,   b o ,   and b c are the biases for the input gate, forgot gate, output gate, and cell state, respectively. The t a n h and σ here represent the tan hyperbolic and sigmoid (logistic) activation functions.

The CNN-LSTM model is an ensemble of CNN and LSTM models (Wang et al., 2016). As shown in Figure 5C, the CNN model first searches for spatial information in high-dimensional input data and transforms it into one-dimensional data. The one-dimensional data is then fed as an input to the LSTM model. Here the CNN network acts as a spatial feature extractor.

This model is mainly designed to improve the performance of the simple LSTM model. The Bi-LSTM model trains the two parallel LSTM layers simultaneously (Cui et al., 2018). The models train one of the parallel layers in the forward direction of the input data and another layer in the backward direction of the input data (see Figure 6A). The Bi-LSTM model could learn more patterns from the input data than the simple LSTM model in this forward and backward training method. As the input data grows, the Bi-LSTM model identifies the unique pieces of information from the data.

Figure 6A depicts the Bi-LSTM design, which consists of two simple LSTM layers. One layer of the model trains the model forward, while a second layer of the LSTM trains the model backward. The parallel layers of the LSTM model receive the same input data and combine their outputs as one output. Finally, in order to forecast the output, the Bi-LSTM model is linked to a dense layer.

The stacked LSTM could be developed by stacking the two simple LSTM layers. The first layer receives the input from the input layer and provides the input to the next connected LSTM layer (see Figure 6B) (Yu et al., 2019). Stacking multiple LSTM layers on top of each other allows the model to learn different temporal patterns from various timestamps in the input data. This design gives more power to the LSTM models to converge faster.

Figure 6B shows a stacked LSTM with two layers stacked on top of each other. The first layer of the stacked LSTM model is designed to take data from the input layer and process it before passing it to the next layer. The next layer is linked to the dense layer, which processes the output of the first layer before passing it to the dense layer. Finally, the dense layer forecasts the required outputs. The primary equations of the model, which update the model’s initial layer, are as follows: i l < t > = σ ( W l x i h l − 1 < t > + W l h i h l < t − 1 > + W l c i ο   c l < t − 1 > + b l i ) (11) f l < t > =   σ ( W l x f h l − 1 < t > + W l h f h l < t − 1 > + W l c f ο   c l < t − 1 > + b l f ) (12) c l < t > = f l < t > ο   c l < t − 1 >   i l < t > ο   t a n   h   ( W l x c h l − 1 < t > + W l h c h l < t − 1 > + b l c ) (13) o l < t > =   σ ( W l x o h l − 1 < t > + W l h o h l < t − 1 > + W l c o ο   c l < t − 1 > + b l o ) (14) h l < t > = o l < t > ο   t a n   h   ( c l < t > ) (15) Where, the variable x _t in the equations is representing the input data sequence at timestamp t. The matrices W l x i , W l c i ,   W l x f ,   W l h f ,   W l c f ,   W l x c ,   W l h c ,   W l x o ,   W l h o , and   W l c o are the weight matrices between two different layers at the level l. Similarly, the b l i , b l f ,   b l o ,   and b l c are the biases at the level l for the input gate, forgot gate, output gate, and cell state, respectively.

An ensembled version of a bidirectional LSTM and a stacked LSTM (called the BS-LSTM network) is a newly designed model for sequence forecasts. A bidirectional LSTM network is concatenated with a stacked LSTM network (see Figure 6C). First, the bidirectional LSTM network is trained in both directions of the input time series (week one to week 62 and vice versa). Second, the output of the bidirectional LSTM layers is linked to a dense layer. Next, the dense layer’s output is provided as input to the stacked LSTM layers. Finally, the stacked LSTM layers relate to a dense layer. The final dense layer forecast the required outputs.

As seen in Figure 6C, seven layers make up the structure of the BS-LSTM model. The input layer was the first layer, and the following two parallel layers were trained in the forward and backward directions. Next, the output of the Bi-LSTM model was related to the dense layer. The dense layer provided the values to the input layer of the stacked LSTM. In the design of this stacked LSTM, two LSTM layers were stacked, one on top of the other. The output of the last stacked layer was related to the dense layer. Finally, the dense layer of the stacked LSTM forecasted the following week’s soil movements.

The first layer in different models was the input layer. The input layer’s dimension could be determined by the features in the dataset, look-back periods, and batch size. The Tangni dataset had fewer data points than Kumarhatti. Thus, the batch size was selected as 16 for Tangni and 1,024 for Kumarhatti. The range of the look-back periods was estimated from the ACF and PACF values. As per the ACF and PACF values, the look-back period for the Tangni dataset was varied from one to five. Similarly, the look-back period for the Kumarhatti dataset was varied from one to forty-two. The number of the features in both datasets was three (i.e., borehole, depth, and soil movement). The second layer was the hidden layer. The nodes in the hidden layer of the LSTM models are called LSTM units. In this research, the one-step forecasting method has been used to forecast the soil movement at the next timestamp.

The number of LSTM units in the hidden layer was changed between 1 and 400 with a step size of 50. The hidden layer’s output vector size was the same as the number of LSTM units in this layer. The dimension of the model’s output could be changed by a dense layer. The dense layer consists of several neurons with a linear activation function. The dense layer’s output vector size was the same as the number of neurons in this layer. Thus, the dense layer with one neuron was connected to the last hidden layer. Table 2 shows the various parameters used by this model to forecast the soil movement in the Tangni and Kumarhatti datasets.

In this model, the output dimension of the convolution layer was determined by the batch size, new number of rows, new number of columns, and the number of filters. In the convolution layer, the batch size was selected as 16 for Tangni and 1,024 for Kumarhatti. The new number of rows and columns were set as one and two, respectively. The number of filters was selected as 62. The kernel size of the convolution layer was set to 1 × 2 . The output vector dimension of the convolution layer was multidimensional, so a flatten layer was added after the convolution layer. The flatten layer converted the multidimensional data into one dimension. The dimension of the model’s output could be changed by a dense layer. The dense layer consists of several neurons with a rectified linear unit (ReLU) activation function. The dense layer’s output vector size was the same as the number of neurons in this layer. Thus, the dense layer with one neuron was connected to the last hidden layer. Table 2 shows the various parameters used by this model to forecast the soil movements in the Tangni and Kumarhatti datasets.

In this model, the time distributed convolution layer executed the same convolution operation to each timestamp as the LSTMs unrolled. The output dimension of the convolution layer was determined by the batch size, new number of rows, new number of columns, and the number of filters. In the convolution layer, both new rows and columns were set as one. The number of filters was selected as 62. The kernel size of the convolution layer was set to 1 × 1 . A max-pooling layer was used to reduce the spatial dimensions. The pool size of the max-pooling layer was set at 2. The max-pooling layers’ output was multidimensional, so a flatten layer was added after the max-pooling layer. The flattening layer converted the multidimensional data into one dimension. The flatten layer provided the one-dimensional data to the next hidden layer. The next layer was the hidden layer. The nodes in the hidden layer are called LSTM units. The number of LSTM units in the hidden layer was changed between 1 and 400 with a step size of 50. The hidden layer’s output vector size was the same as the number of LSTM units in this layer. The dimension of the model’s output could be changed by a dense layer. The dense layer consists of several neurons with a rectified linear unit (ReLU) activation function. The dense layer’s output vector size was the same as the number of neurons in this layer. Thus, the dense layer with one neuron was connected to the last hidden layer. Table 2 shows the various parameters used by this model to forecast the soil movements in the Tangni and Kumarhatti datasets.

In this model, the nodes in the hidden layer was called LSTM units. The number of LSTM units in both hidden layers was changed simultaneously between 1 and 400 with a step size of 50. The hidden layer’s output vector size was the same as the number of LSTM units in this layer. The dimension of the model’s output could be changed by a dense layer. The dense layer consists of several neurons with a linear activation function. The dense layer’s output vector size was the same as the number of neurons in this layer. In this research, the one-look-ahead forecasting method has been used. Thus, the dense layer with one neuron was connected to the last hidden layer. Table 2 shows the various parameters used by this model to forecast the soil movements in the Tangni and Kumarhatti datasets.

In this model, the number of LSTM units in the hidden layer was changed between 1 and 400 with a step size of 50. The hidden layer’s output vector size was the same as the number of LSTM units in this layer. The dimension of the model’s output could be changed by a dense layer. The dense layer consists of several neurons with a linear activation function. The dense layer’s output vector size was the same as the number of neurons in this layer. Thus, the dense layer with one neuron was connected to the last hidden layer. Table 2 shows the various parameters used by this model to forecast the soil movements in the Tangni and Kumarhatti datasets.

This model has seven layers, three in the Bi-LSTM model (Input: one; Hidden: one; Output:1) and four in the stacked LSTM model (Input: one; Stacked: two; Output:1). The architecture of the Bi-LSTM and stacked LSTM layers was the same as Bi-LSTM and stacked LSTM model developed in this paper. Table 2 shows the various parameter values used by this model to forecast the soil movement in the Tangni and Kumarhatti datasets. For example, the BS-LSTM parameters for the Tangni dataset were varied as per the following: the batch size was fixed at 16; the LSTM units in the hidden layers were selected between 1 and 400 with a step size of 50 for the Bi-LSTM model and changed between 1 and 400 with a step size of 50 for both stacked layers in stacked LSTM; the number of epochs was changed as 10 or 50; the look-back period was varied from one to five; and, the inputs were passed with shuffling turned on or turned off. The BS-LSTM parameters for the Kumarhatti dataset were varied as per the following: the batch size was fixed at 1,024; the hidden layer’s LSTM units in Bi-LSTM and stacked LSTM were varied between 1 and 500 with a step size of 10; the number of epochs was changed as 10 or 50; the look-back period was varied between 1 and 42 (as estimated from ACF and PACF); and, the inputs were passed with shuffling and without shuffling.

Before entering the time series data into models, we divided the input data into several packets based on borehole and depth. Packets from the first 80% of data were utilized to training the models, and packets from the last 20% of data were used for testing the models. Each packet was a combination of X, Y, where X was the predictor, and Y was the forecasted soil movement value (see Figure 7). As shown in Figure 7, the X predictor formed a movement vector consisting of soil movements recorded by a sensor in a borehole in a timestamp (M), borehole number of the sensor (Borehole), and the depth of the sensor (Depth). The M was the timestamp value of the soil movement recorded by a sensor in a borehole in a particular time. The Borehole value was one of the five boreholes at the site. The Depth value was the depth of the sensor in a specific borehole. For example, for a look-back period of three in the Tangni dataset, packet X₁ may contain the movement recorded by the sensor in the first 3 weeks in borehole one at a depth of 3 m. The corresponding Y₁ contained the actual value of the movement recorded by the same sensor in borehole one in week four, forecasted by a model. Thus, there was a one-look-ahead soil movement forecast (Y) for a certain look-back period and for a particular sensor at a specific depth (X) (the look-back was passed to models as a parameter). These X and Y packets could be shuffled before input to a model, where the shuffle operation would shuffle the two lists (i.e., X and Y) in the same order (see Figure 7).

The Tangni dataset has a limited amount of data, where the training dataset has only 62 data points to train LSTM models. When the training dataset is limited, the system can overfit the model’s parameters during training. The dropouts can be applied to the layers of the LSTM to prevent overfitting (Pham et al., 2014; Gal and Ghahramani, 2016). Different combinations with the probability (p) of dropout were applied on the input-output and recurrent connections of the LSTM layer. The probability value (p) was varied between 0.0 and 0.8 with step size 0.2, where a p-value of 0.0 means no dropout applied. For example, a combination (0.2, 0.8) represents 20% dropout applied to the LSTM unit’s output and 80% dropout applied on the recurrent input of the LSTM unit (Gal and Ghahramani, 2016).

The soil movement forecasting is a regression problem, where the resulting soil movement is assumed to have a floating value. Thus, an error can be calculated between the true value and the forecasted value of the soil movements. The different performance measures have been used to evaluate the performance measure of the various models (Behera et al., 2018). In this paper, four performance measures are used: mean relative error (MRE), root mean square error (RMSE), normalized root mean squared error (NRMSE), and mean absolute error (MAE). Eqs 16–19 were used to calculate the values of these measures to find the difference between actual points of datasets from the forecasted points of datasets. R M S E = 1 n ∑ i = 1 n ( T r u e   a n g l e − F o r e c a s t e d   a n g l e ) 2 (16) N R M S E = Σ i = 1 n ( T r u e   a n g l e − F o r e c a s t e d   a n g l e ) 2 Σ i = 1 n ( F o r e c a s t e d   a n g l e ) 2 (17) M A E = 1 n ∑ i = 1 n ( | T r u e   a n g l e − F o r e c a s t e d   a n g l e | ) (18) M R E = 1 n ∑ i = 1 n ( | T r u e   a n g l e − F o r e c a s t e d   a n g l e F o r e c a s t e d   a n g l e | ) * 100 (19) where, n denotes the number of data points in the Tangni or Kumarhatti dataset, the true angle denotes the actual observed value of the soil movements in the dataset, and the forecasted angle denotes the forecasted value of the soil movements by the model.

We created a grid search procedure to calibrate the parameters in various models. In this process, we changed the different set of parameters (described in Table 2) in a LSTM model. After feeding a combination into the LSTM models, we recorded the MAE, RMSE, NRMSE, and MRE, and we chose the parameters with the lowest error in the model.