Inﬂuence of Meteorological Parameters on Explosive Charge and Stemming Length Predictions in Clay Soil during Blasting Using Artiﬁcial Neural Networks

: The inﬂuence of the meteorological parameters (precipitation and air temperature) during blasting in clay has a direct impact on the success of blasting. In the case of large amounts of precipitation (rain and snow) recorded in the subject area, blasting in clays cannot be carried out due to the grain of the clay and the inability to access the subject area. Moreover, the air temperature in the subject area affects the blasting performance. The most ideal temperature for blasting in clays is between 15 and 25 ◦ C because then the clay has the best geotechnical characteristics. The research was conducted on the exploitation ﬁeld Cukavec II, which is located near the city of Varaždin in the Republic of Croatia. Amount of precipitation and air temperature were considered to obtain the best blasting effect. Inﬂuence of meteorological parameters on the amount of the explosive charge and stemming length when blasting in clays was demonstrated via models based on Artiﬁcial Neural Networks (ANN). The ANN model network consists of a Long Short-term Memory (LSTM) part to process time dependent meteorological data, and fully connected layers to process blasting input data. Two types of explosive charges were compared, Pakaex and Permonex V19.


Introduction
Research conducted in the exploitation field Cukavec II ( Figure 1) in the period from 2014 to 2016, proved the dependence of the amount of explosive charge and stem length on the formation of the volume of the spherical expansion, resulting from the detonation of the explosive charge in a borehole of 131 mm in diameter. The main goal was to expand the knowledge about the possibilities of using explosive charges in geotechnical practice. This especially refers to blasting in clay soil by which spheres or cavities of other shapes are formed at different depths below the soil surface by activating a certain type and mass of explosives [1,2].
The structure of the deposit on the exploitation field Cukavec II is dominated by kaolinite, quartz, feldspar, and chlorite, and the granulometric composition shows that the mineral raw material contains about 75% of the clay component, while the rest is silt. Previous research was conducted with the goal of examining geotechnical properties of the clay soil, where cohesion c, internal friction angle ϕ and volume weight γ were determined. Cohesion is c = 23.4 KN/m 2 , internal friction angle ϕ = 19.8 • , and volume weight γ = 18.7 KN/m 3 . Due to the characteristics of hydroalumosilicate clays, which means that it absorbs water and therefore becomes brittle, blasting is directly affected [2,3].
Mass ranges of two different types of explosive charges (Permonex V19 and Pakaex) which achieve spherical expansions in clay soil were measured and determined. Experimental blasting was performed on a total of 48 boreholes, out of which 6 were test boreholes ( Figure 2).  The expansions created by the blasting in the boreholes were different in shape and size, but they were all basically in the shape of a sphere ( Figure 3). Experimental blasting has determined the effective ranges of explosive charge masses of Permonex V19 and Pakaex which achieved spherical expansion in clay soil. It has been found that spherical expansion for a given borehole diameter of 131 mm and a borehole depth of 2.00-3.00 m can be achieved by explosive charge of Permonex V19 and Pakaex, ranging from 0.2 to 1.6 kg [2,5]. The most significant result was certainly the determination of the volume of the resulting expansions in the clay soil during the blasting of a series of boreholes placed on a certain profile ( Figure 2). To calculate the resulting expansion volume, the application "Bušotine" was designed [2,5]. In addition to the application calculating the resulting expansion, it draws 2D and 3D models of it ( Figure 3). Further, the compatibility of the application with other CAD tools greatly contributes to the verification and more detailed 3D model of the resulting expansion. Regarding the need to use the data processing method for the obtained data, Volume of the resulting cavity, V rc [m 3 ], Resulting expansion of the borehole, L re [m], and Deepening of the resulting expansion, D re [m], during the preparation of the doctoral dissertation, Težak, D. 2018 [2] confirms and additionally proves this by introducing an additional supplemented model [2,6]. The introduction of a supplemented model in defining the resulting volume of spherical expansion, expansion and deepening of a borehole improved the original method presented in Težak, D. 2018 [2], and further facilitated the use of the application "Bušotine" when blasting in clays in geotechnical practice.
Since the clay is hydroalumosilicate, which means that it absorbs water and therefore becomes brittle, the blasting in the clay can only be carried out in dry months and when the air temperature is between 15 and 25 • C. The data obtained also showed that there is a dependence of air temperature and precipitation on the time period in the year in which clay exploitation is possible, and guidelines for successful clay exploitation are presented [7,8]. The amount of precipitation has a direct impact on the success of blasting in clay, because when the precipitation in the clay soil during detonation of explosive charge increases, the pressure wave of the explosion has a stronger impact and larger spherical expansions volume occur [9][10][11].
In the field of geotechnics and blasting in clay soils, artificial intelligence (AI) was mainly investigated in terms of soil properties prediction models, i.e., the study [12] used genetic programming (GP), support vector regression and ANN-s to model the relationship between volumetric water content in soil and input parameters in the form of solid density, initial water content, and soil suction. Similarly, in study [13], generalized regression neural networks (GRNN) and backpropagation neural networks (BPNN) were used to model soil temperature based on half-hourly gauged air temperature, relative humidity, wind speed, vapor pressure deficit, and solar radiation, while in the study [14], an ANN-based predictor was developed to model dependence between current hourly soil temperature and hourly precipitation rates, previous average soil temperature and average air temperature. Another study [15] used adaptive network-based fuzzy inference system (ANFIS) to predict and evaluate ground vibrations and frequencies caused by blasting in a tunnel tube. Different kernel based Support Vector Machine (SVM) models were used in study [16] to predict peak particle velocity induced by mine blasting, while the model inputs depend on results of boosted-Chi-Squared Automatic Interaction Detection (CHAID).
Scarcity of research has been identified that investigates reliable predictions of the influence of air temperature and precipitation on blasting performance in clay soil. By using Artificial Neural Networks (ANN) with meteorological input data in the form of precipitation and air temperature in addition to blasting related input data, it is possible to determine the optimal amount of explosive charge and the length of the borehole stem for successful blasting in clay soil.
However, when dealing with sequential data, as i.e., time series data like precipitation and/or air temperature, often a time dependence is present (data at time t, depends on the data preceding it), so the application of a traditional feedforward neural network is limited [17]. Therefore, when processing timeseries data (like i.e., meteorological measurements) Recurrent Neural Networks [18] (RNN), and especially the more advanced version of RNNs, namely, Long Short-term Memory (LSTM) [19] networks, are often chosen. Several applications of LSTM networks in hydrology [17,[20][21][22], remote sensing [23,24] and air quality modelling [25] tasks have been found. However, to our knowledge no published research investigates the implementation of LSTM neural networks and meteorological data influence on blasting in soft clay.
Therefore, the integration of ANNs and LSTMs in the field of blasting to predict the amount of the explosive charge and stemming length is presented for the first time. The predictions are done by processing meteorological parameters such as daily precipitation rates, daily air temperature and/or relative humidity and blasting related parameters, i.e., borehole depth, volume of resulting borehole cavity, resulting expansion of the borehole etc., for successfully conducted blasting in soft clay. By processing these parameters, it is desired to prove that the use of Artificial Neural Networks can reliably determine the amount of explosive charge and stemming length that will be used for blasting in clays in geotechnical practice.

Study Area and Data
The research was divided into three phases [1]. The first phase included a study conducted in 2014 where the blasting training area consisted of the two test boreholes PMB1 and PMB2, and eight boreholes, MB1-MB8 ( Figure 2). The second phase was carried out in 2015 and consisted of the two test boreholes, PMB3 and PMB4, and 24 boreholes, MB13-MB36 ( Figure 2). The third phase was conducted in 2016 and consisted of the two test boreholes, PMB5 and PMB6, and nine boreholes, MB37-MB45 ( Figure 2).
The test boreholes were used to gain a better insight into the behaviour of the explosive charges of Pakaex and Permonex V19 used for this purpose and to define the length of the stemming. Based on the two test boreholes and eight boreholes from 2014, it was concluded that the best sand stemming is from 0.5 m to 1.00 m long.
The spherical expansions of boreholes were created by the detonation of the Ammonia Nitrogen Powder Explosive Permonex V19 and ANFO explosives commercially known as Pakaex. The velocity of detonation (VOD) is a detonation parameter that can be used to determine the performance of used explosive in different rock types also in a coherent clay. The basic characteristics of Pakaex explosive charge are density 0.87 g/cm 3 , VOD 2950 m/s, and gas volume 984 L/kg, and of explosive charges Permonex V19 are density 0.95 g/cm 3 , VOD 4500 m/s, gas volume 900 L/kg. Furthermore, it should be noted that during the research, the boreholes PMB4, MB33, and MB30 were destroyed, by the work of mechanization in the exploitation field Cukavec II [3].
As  Table 1. Based on the data displayed in Table 1, the volume of the resulting cavity, V rc [m 3 ], and resulting expansion of the borehole, L re [m], at the same explosive charge mass, Q [kg], were higher after the activation of Permonex V19 than Pakaex because Permonex V19 has a higher density and VOD. Moreover, deepening of the resulting expansion, D re [m], at the same explosive charge mass, Q [kg], is larger with Pakaex than with Permonex V19.
Obtained data was used to define resulting spherical cavity in clay after blasting, during which guidelines were given to geotechnicians and builders on how the spherical cavity could be used in practice. The main effects of blasting in soils, particularly clay, are shown in the installation of structural elements for anchoring of foundations and retaining walls, permanent stabilization of clay slope and stabilization of various commercial structures such as overhead transmission towers, tunnels, etc. [2,3].
The meteorological input data in the form of daily precipitation, minimum and maximum daily air temperature, and relative humidity was acquired from Croatian Meteorological and Hydrological Service [26] in the form of tabular data, which was extracted and formatted as timeseries data suitable for input to ANN-s. The data was acquired for two gauges nearest to the exploitation site, main automatic station Varaždin, and climatological station Novi Marof.

Artificial Neural Networks (ANN)
An artificial neural network is a data processing system built of simple processing units (neurons) with the tendency of storing experience for later use [27]. Each ANN consists of neurons interconnected by directed communication links (synaptic weights of neurons). The weights are calibrated when solving tasks in which the network acquires "knowledge" by learning from available data [27]. The most common networks are socalled Feedforward (fully connected) or shorter just ANN, which consist of one input layer, one or more hidden layers, and one output layer which are fully connected (each neuron of the previous layer is connected with each neuron in the next layer).
The mathematical process in each neuron k of the network can be described by the following two equations [27]: and where: x 1 , x 2 , . . . , x m are the inputs; w 1 , w 2 , . . . , w km are the synaptic weights of the neuron k; u k is the linear combiner output due to the inputs; b k is the bias term; ϕ is the activation function; and y k is the output of the neuron k. A visualization of the processes (Equations (1) and (2)) in the neuron is given in Figure 4. When data is delivered to an untrained network, the output differs from the expected value because the synaptic weights (w) and bias (b) are not yet adjusted to the data. The process of adjusting the network weights and biases is called training. During the training phase, input data is delivered to the network (forward propagation process), the outputs are calculated and compared with the observed (expected) output, the loss is calculated using the target (loss) function. The type of function depends on the type of the problem, for regression tasks a popular choice is the mean squared error (MSE) function, while for classification tasks the cross-entropy loss function is often chosen. The calculated error propagated backwards (backpropagation process) through the network by calculating the gradients of the loss function (gradient descent). Depending on the loss value during gradient descent, the weights and biases get updated to minimize the error. The procedure is repeated either a fixed number of times (fixed number of epochs) or until certain conditions are met, e.g., the network has reached a predefined error value.

Recurrent Neural Networks (RNN)
Recurrent Neural Networks (RNNs), however, are especially successful when dealing with sequential data. In opposite to a fully connected layer from an ANN, the recurrent layer in an RNN, consist of a loop (a chainlike structure) which enables information flow "through time". A neuron handling data from timestep t-1, provides an output while also passing information ("Hidden State") to the neuron handling timestep t, and so on. The result is better handling of data where the current timestep depends on data from several timesteps before.
The calculations at each timestep of the input sequence can be represented with the two following equations [28,29]: and where: W h is the weight matrix between the input and hidden layers and U h is the recurrent layer weight matrix at adjacent time steps. The vectors b h and b y are the bias terms which give an offset to each node (neuron). However, in real world data, a situation where the data multiple time steps before influences the data at the current moment often occurs. No matter the effectivity of RNNs when dealing with timeseries data, in such cases, problems occur in the form of vanishing or exploding gradients. With the increase of the time dependency gap, the performance of traditional RNNs drops significantly [18,30].

Long Short-Term Memory Neural Networks (LSTM)
During the late 1990s, Hochreiter & Schmidhuber [19] provided the solution for exploding and vanishing gradients in the form of an upgraded variant of the RNN, namely, Long Short-term Memory (LSTM) which is capable of dealing with long-term dependencies. To overcome the problems associated with long-term dependencies, the LSTM cell ( Figure 5) contains a so-called "Cell State" which enables almost unobstructed information flow through time. Further, the LSTM cell has three building blocks, called "gates", which regulate the data flow through the cell. The "Forget Gate" (Equation (5)) regulates how to modify the "Cell State" according to previous (t-1) and current (t) timestep.
where: denotes the Hadamard product (elements-wise vector multiplication), σ denotes sigmoidal function, W are the network weights and b are the biases, x(t) is the input at current timestep, while h(t-1) is the "Hidden State" (same as for an RNN cell) of the previous timestep. Next, the "Input Gate" (Equations (6) and (7)) regulates how much of the current information (timestep t) will be added to the "Cell State". The "Cell State" gets modified by the "Forget" and "Input Gates" following the Equation (8). The "Output Gate" (Equation (9)) regulates the outflow of the data from the cell, which is affected by the current input data and the "Cell State". The final output of a LSTM cell is calculated following the Equation (10) and is also the new "Hidden State" which is passed to the next timestep.

Model Evaluation Metrics
To enable comparison with similar research, the results of the models are evaluated by two very popular statistical evaluation metrics in regression tasks, namely, the Mean Squared Error (MSE) and the coefficient of determination (R 2 ). The MSE evaluates the relative distance of the predicted and observed variable, where values closer to zero indicate lower bias between the predicted and observed value, indicating better model performance.
where: y i is the observed variable at time t, y i is the predicted variable at time t, and n is the length of the dataset. The coefficient of determination R 2 describes the square of the sample correlation coefficient, ranging from 0 to 1, and measures the observed variance explained by the model. Values near 0 indicate no correlation, while values near 1 indicate good model performance, meaning that the model can explain the variance in the observed data.
where: y i is the mean of the observed data and y is the mean of the predicted data. Relative squared error (RSE) is a relative measure, denoting it what it would have been if a simple predictor was used. It takes the total squared error and normalizes it by dividing it by the total error of a simple predictor. where:ŷ The lower the RSE, the better the model, so for a perfect fit model, the values of RSE andŷ are both equal to zero.
Mean bias error (MBE) denotes the mean forecast error, due to overcast or undercast, meaning it measures how much does the model overpredict or underpredict the values.
Lower values are desirable, where negative MBE means the model predicts lower than measured, while a positive MBE means the model predicts larger values, respectively. The leverage statistic in regression analysis is often used to quantify the influence of certain data points on the regression line, referred as high leverage points. The studentized residuals are calculated as follows:ε = y − y The average leverage calculated as: where: p is the number of predictor variables. The specific leverages of each predicted data point are calculated from the weight matrix, by following: Or written as linear combination: where: h in are the leverages of each data point. Each data point whose leverage is 3 times higher than is flagged as "high leverage point". The results of leverage analysis are then plotted as leverage versus studentized-residuals plots.

Model Structure and Training Process
Research is based on open-source software and libraries. Python [31] is the programming language of choice. In addition, NumPy [32], Pandas [33], statsmodels [34], Matplotlib [35] and Seaborn [36] libraries are imported for processing, management and visualization data. The model based on an ANN was developed in Pytorch [37], an open source library available online at: https://pytorch.org (accessed on 17 June 2021).
The ANN model used in this study consists of two parts: the first part handles the sequential (time series) meteorological data and is based on two LSTM layers and two fully connected layer. The goal of this part of the model is to provide the "state" of the clay according to the meteorological data. The output of the first part serves as an input to the second part of the model, consisting of two fully connected layers with ReLU activation, handling clay blasting data input data (8 features).
The output of the whole model is a 2D vector (two predicted features), for the purpose of this research the model is trained to predict the Stemming Length [m] and Explosive Charge Mass [kg]. However, the approach is flexible and enables on-the-go changing of input and/or output feature number, meaning, the number of input and/or output features can be changed by desire. For the purpose of this research in order to evaluate the influence of meteorological data on blasting in clay soil, we ran two versions of the model. The "model_1" had only 2 meteorological input features (precipitation and mean daily air temperature from meteorological station Varaždin), while the "model_2" had 8 input features (precipitation, minimum and maximum air temperature recorded at meteorological stations Varaždin and Novi Marof, and additionally relative humidity recorded at meteorological station Varaždin). Both models were run in Jupyter Notebooks on an PC with Intel i5 8400 processor, 16 GB of RAM and RTX 2060 Super discrete GPU with 8 GB of video memory. Both models were trained on the GPU, with same hyperparameters, 4 batches of data per iteration, 5 x 10 −5 learning rate, 120 hidden nodes in the LSTM layers with 0.3 dropout and 5 nodes in the fully connected hidden layers. Due to a regression problem, Mean Squared Error (MSE) loss function was used, together with the AdaMax [38] optimizer handling the gradient descent during backpropagation process when the network was trained.
To acquire best possible model results, several trial and error runs were undertaken, with the goal of finding best possible combination of hyperparameters. Training was considered successful right before the model started to overfit to training data, so model_1 was trained for 50 epochs while model_2 was trained for 60 epochs; due to more input features, the network needed more iterations to adapt to training data. To enable results reproducibility, torch.manual_seed has been used and set to 14.

Results
The The charts of explosive charge mass versus the resulting borehole expansion (Figures 6 and 7) and volume of resulting cavity (Figures 8 and 9) clearly indicate that it is possible to use ANN-s to predict the required mass of explosive charge, considering a desired expansion of the borehole or volume of the resulting cavity. The good visual correlation is further substantiated by the statistical metrics in the form of MSE and R 2 ( Table 2) where only the predicted results for the testing period of model_2 is below 0.8. Similarly, the MSE worst value is 0.054 kg 2 during validation period of model_2. The overall low values of the RSE for the independent validation and testing periods indicate the good generalization ability of the proposed models in terms of predicting charge mass. The MBE statistic in case of model_1 indicates that charge mass in all three periods gets underpredicted, meaning that the model predicted slightly lower values than expected. Model_2 corrects this behaviour, where only for validation period we measure slightly lower underpredictions, and for training and test period slight overpredictions are encountered, although with slightly higher error margin than for model_1 (RSE metric).    On the other hand, it is considerably harder to identify a clear correlation between the predicted and observed values for both the cases of scatter charts of stemming length versus borehole expansion (Figures 10 and 11) and similarly for borehole cavity (Figures 12 and 13). However, model_2 yields significantly better results in terms of the stemming length predictions (visually) and supported by the statistical metrics. Overall, the performance of model_2 in terms of stem length is slightly better than model_1, through all four evaluated metrics. Considering the MBE metric, model_1 underpredicts slightly during training phase, while it overpredicts the stem length for both validation and testing phase. The behaviour changes in model_2, where during the test phase a slight underprediction occurs.     Considering the average leverage of 0.681 for model_1 and 1.09 for model_2, a conclusion can be made that although the plots show a relative high leverage point during the validation phase, with leverage of 0.653 for explosive charge (Figures 14 and 15) and 0.518 for the stem length (Figures 16 and 17), both the predicted variables and all four models, are under the 3 times leverage threshold.

Discussion
The presented results clearly indicate the importance of meteorological parameters during blasting in clay. It is important to emphasize that adding more input features to an ANN does not necessarily yield better prediction results. Adding additional input data also adds to complexity of the model, which consequently leads to more complicated training procedure of the ANN. This is visible, where adding more input features required increased training duration (more epochs). Interestingly, adding more meteorological data improved the stemming length predictions, but somewhat lowered the explosive charge mass prediction quality, yielding more than satisfactory results, nonetheless.
The weak results of stemming length prediction have to be attributed to the input data, since most of the blasting experiments were done with 0.5 m or 1 m stemming length. Such results are somewhat expected and reflected by the statistical metrics, especially in high values of RSE during training phase for stem lengths. The error decreases during validation and testing phases due to larger range of different stem lengths used during those two phases. For a classical data interpretation, this can be favourable; however, in case of predicting this exact parameter using ANN-s this is not the best possible input parameter. The coefficient of determination values are not particular high; therefore, the results of stem length predictions should be taken with care. Additionally, in situ experiments with different stem lengths will be considered in possible future research. Nonetheless, the network managed to produce results of some value. Such results can be a guideline for new blasting, in terms of which stemming lengths could be used as experiments to acquire valuable data for improving such models.
On the other hand, the predictions result of the explosive charge mass are very satisfactory and prove that such a modelling approach using ANN-s should be considered before field experiments, in order to avoid the high cost of experimental borehole drilling and blasting. The coefficient of determination results of >0.8 indicate that the predictions of charge mass are more than satisfactory. The bias values of model_1 indicate slight underprediction of values but are fixed with the enhanced model_2. The RSE values are overall low, for both models, so even a model with just one meteorological station with just precipitation and mean air temperatures can be considered sufficient for charge mass predictions. It is also proved that weather, and meteorological data in particular, has significant effects on blasting. This study should serve as an idea for similar studies in this field due to its cost effectiveness and simple data requirements.

Conclusions
Presented in this research is a novel approach in the field of blasting by implementing Artificial Neural Networks as model for predicting the required explosive charge mass and stem length to achieve the desired volume of borehole cavity or borehole expansion. A crucial part of our model is the influence of meteorological data which has an impact on the state of the clay, and therefore impacts the blasting procedure and results, respectively. By creating two models forced with a different number of input features in the form of meteorological data-precipitation, air temperature, and relative humidity-the goal was to increase the prediction accuracy of blasting related parameters, in this case stem length and charge mass.
The results of charge mass prediction are very satisfactory, both models performed good, with coefficient of determination R 2 values ranging from 0.77 for test phase of model_2 to 0.99 for validation phase of model_1. Interestingly, model_1 underpredicted the charge mass for all three phases (MBE values, Table 2), while the same behaviour was only noticed in validation phase for the second model with more meteorological data. Adding more features of meteorological data lowered the overall model bias slightly, albeit causing some overpredictions in case of model_2. Such results are probably due to lower count of blasting related data points in the validation phase, where one point shows somewhat higher leverage in relation to other points from the set, although lower than the overall threshold value of leverage.
Stem length prediction results show lower overall reliability, reflected in the coefficient of determination values R 2 ranging from 0.106 for validation to 0.753 for training phase of model_1, respectively. The values of odel_2 range from 0.224 for the test phase to 0.544 for the training phase. Furthermore, the stem length results show high results of squared errors, which is another factor indicating lower model performance. Such results are probably due to the lower number of different stem lengths available for training and can easily be fixed with addition of more in situ data. Currently, this is not feasible due to financial reasons; however, it will be certainly investigated in future research regarding this topic.
Overall, the results prove that it is possible to use ANN-based data-driven models in the field of blasting in terms of predicting required explosive charge mass and stem lengths. Care has to be taken to choose data with a proper range of different values for the training phase, which is a common drawback of data-driven models solely dependent on input/output combinations of data. Further, the results show that meteorological data and conditions influence the blasting in clay and should definitely be considered before any blasting experiments in order to maximize the effectiveness and avoid any potential risks. With the study, we have shown that adding more meteorological data in the form of more features does not necessarily enhance model performance. However, in the case of stem lengths, it resulted in more consistent results, with lower range of R 2 values. This research can and should serve as an idea for further investigations and implementations of ANN predictor models in the field of blasting.