Temperature Prediction Model in the Main Ventilation System of an Underground Mine

A model to forecast the underground temperature in a mine ventilation circuit was developed on the basis of a case study and actual data describing temperature, airflow, and drift length collected over several years. A mathematical model featuring seven variables with interactions provided reliable predicted temperatures, achieving a correlation of R2 = 0.933 with an estimation error of ±2 °C. Its soundness was proven using both the node-to-node analysis and the multi-node approach. The multi-node approach was shown to be an interesting option to model underground mining environments. This model can be very useful to predict the temperature evolution along the main ventilation system, determine the best workplace conditions in terms of temperature, and analyze different planning scenarios of the mine. Moreover, some recommendations are presented for obtaining reliable data when using temperature sensors and the model in a U-shaped ventilation system.


Introduction
The environmental conditions in an underground mine are crucial in terms of efficiency and health and safety conditions [1,2]. Furthermore, these factors have to be taken into account for mine planning due to the restrictions imposed by proper ventilation conditions. This becomes increasingly relevant as shallow deposits run out and mineral resources are extracted at greater depth, and as mines increase in size with more infrastructure, with an intensification in the use of mining equipment [3,4] and human resources [5]. This creates a situation where the potential effects must be defined [6], with a need to determine the best possible conditions related to performance in the workplace [7][8][9][10].
One of the most important variables regarding environmental conditions in an underground mine is temperature [11][12][13]. It not only affects the concentration levels of the workers [14], but also creates unsafe conditions [15,16]. Currently, there is an increasing interest in the characteristics of underground spaces [17], with previous studies creating a model to predict the temperature at the bottoms of shafts in metal mines [18], simulating the temperature at the working face of a coal mine [19], and determining the main parameters of heat from the intrinsic properties of the deposit [20]. In this regard, airflow is crucial to maintain adequate temperatures in the workplace, ensuring its supply during the whole lifespan of the mine and any future changes in the ventilation circuit [21].
There were even some studies focusing on potash mining [3] or applying cooling systems [22][23][24]. However, it is necessary to not only identify the current conditions in terms of temperature, but also to predict the future conditions in order to plan the necessary actions. Some recent studies focused

Case Study
The research was based on an underground mine that extracts potash using an irregular room and pillar method, at a depth between 700 and 1000 m. The layout of the mine involves two main levels with a U-shaped ventilation circuit. The principal ventilation circuit is placed on the lower level, while the higher level contains the workings faces, with several splits of airflow from the level below, as well as the auxiliary ventilation system. The main ventilation system consists of a single exhausting fan. Extraction is achieved via three shifts lasting eight hours each, working continuously throughout the day.
The different exploitation panels include an extension of the lower level and, subsequently, an extension of the same panel above. Therefore, it is important to know, in advance, the temperature evolution along the principal circuit, to determine the environmental conditions at the entrance of the panel and to be able to analyze different scenarios.

Data Collection
Determination of the necessary data required to define the ventilation characteristics of the drifts in a mine was previously considered in a case study [27,34]. Seven control points representing the main ventilation system were established, from the entry to the end of the intake, collecting data every 3 min between January 2016 and June 2019. The following variables were collected: -Time of the measurement (date and hour); -Control point (X, Y, and Z coordinates); -Actual distance between control points (m); -Dry temperature ( • C); -Wet temperature ( • C); Appl. Sci. 2020, 10, 7238 3 of 11 -Effective temperature ( • C); -Airflow (m 3 /s).
Wet and dry temperatures are crucial to understanding the environmental conditions in an underground mine [35]; however, an analysis of the humidity evolution along the intake does not show appreciative changes. Therefore, the effective temperature, described by Spanish law (RGNBSM, ITC 04.7.02), was taken into account as a predictive variable, which includes wet and dry temperatures combined into a single value with a certain proportion, as shown in Equation (1).
where T ef is the effective temperature ( • C), T w is the wet temperature ( • C), and T d is the dry temperature ( • C).

Data Processing and Analysis
The equipment to record and collect data was calibrated regularly to avoid inaccurate measurements [36]. The effective temperature was transformed into daily average values, resulting in a total of 1155 daily records. Figure 1 shows the temperatures measured in each node, corresponding to the fixed control points representing the ventilation conditions of the intake. As can be seen in Figure 1, there was a higher dispersion in the initial nodes due to the influence of outside conditions, which were quite variable depending on the season of the year.
Wet and dry temperatures are crucial to understanding the environmental conditions in an underground mine [35]; however, an analysis of the humidity evolution along the intake does not show appreciative changes. Therefore, the effective temperature, described by Spanish law (RGNBSM, ITC 04.7.02), was taken into account as a predictive variable, which includes wet and dry temperatures combined into a single value with a certain proportion, as shown in Equation (1).
where Tef is the effective temperature (°C), Tw is the wet temperature (°C), and Td is the dry temperature (°C).

Data Processing and Analysis
The equipment to record and collect data was calibrated regularly to avoid inaccurate measurements [36]. The effective temperature was transformed into daily average values, resulting in a total of 1155 daily records. Figure 1 shows the temperatures measured in each node, corresponding to the fixed control points representing the ventilation conditions of the intake. As can be seen in Figure 1, there was a higher dispersion in the initial nodes due to the influence of outside conditions, which were quite variable depending on the season of the year. According to Navarro et al. [37], variations in temperature in an underground mine are caused by several intrinsic and extrinsic characteristics, as shown in Equation (2).
where ΔTac represents autocompression (°C), ΔTrp represents thermal rock properties (°C), ΔTee represents heat emission from electrical equipment (°C), ΔTde represents heat emission from diesel equipment (°C), ΔThm represents human metabolism (°C), and ΔTw represents thermal water (°C). The underground atmospheric temperature, at any location inside the mine, can be expressed as a function of the underground initial temperature, eliminating the autocompression factor from Equation (2), since the temperature at the bottom of the shaft is taken into account, thereby obtaining Equation (3).
where t2 is the underground atmospheric temperature (°C), and t1 is the underground initial temperature (°C). According to Navarro et al. [37], variations in temperature in an underground mine are caused by several intrinsic and extrinsic characteristics, as shown in Equation (2).
where ∆T ac represents autocompression ( • C), ∆T rp represents thermal rock properties ( • C), ∆T ee represents heat emission from electrical equipment ( • C), ∆T de represents heat emission from diesel equipment ( • C), ∆T hm represents human metabolism ( • C), and ∆T w represents thermal water ( • C). The underground atmospheric temperature, at any location inside the mine, can be expressed as a function of the underground initial temperature, eliminating the autocompression factor from Equation (2), since the temperature at the bottom of the shaft is taken into account, thereby obtaining Equation (3).
where t 2 is the underground atmospheric temperature ( • C), and t 1 is the underground initial temperature ( • C). The intake drifts from the case study had very steady conditions in terms of equipment usage, working force, and thermal rock properties, due to them being opened for a long time. It was also considered that the thermal water had no influence, due to the intrinsic characteristics of the type of Appl. Sci. 2020, 10, 7238 4 of 11 mineral extracted in the case study. The physicochemical properties of the air conditions were not considered either. Therefore, the more developed form of Equation (3), with simplifications, can be written as Equation (4).
where ∆T intake depends on the dry and wet temperatures expressed as equivalent temperature ( • C). However, the airflow is not steady throughout the main ventilation circuit, with some leakages existing along the intake and/or return, especially in U-shaped circuits as considered in this study, as well as airflow splits to different extracting panels.
In this case, the cross-section did not have a substantial variation along the intake drift. Therefore, the heat contribution was similar along the intake drifts and the airflow was mainly reduced due to its distribution to the workplace. Distance was also very important to consider because of the interaction between the air and heat sources along the drifts. Equation (5) describes the factors influencing temperature variation along a drift.
where C is the airflow (m 3 /s), and D is the distance (m). A mathematical model is a simplified representation of reality with predictive variables (x 1 , x 2 , . . . , x n ) and target variables that need to be estimated [38]. The model presented herein was based on the variables detailed in the previous paragraph, with the aim of forecasting the temperature using in situ measurements of temperature and airflow along the ventilation circuit, taking into account seven preset control points (see Figure 2).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 11 The intake drifts from the case study had very steady conditions in terms of equipment usage, working force, and thermal rock properties, due to them being opened for a long time. It was also considered that the thermal water had no influence, due to the intrinsic characteristics of the type of mineral extracted in the case study. The physicochemical properties of the air conditions were not considered either. Therefore, the more developed form of Equation (3), with simplifications, can be written as Equation (4).
where ΔTintake depends on the dry and wet temperatures expressed as equivalent temperature (°C). However, the airflow is not steady throughout the main ventilation circuit, with some leakages existing along the intake and/or return, especially in U-shaped circuits as considered in this study, as well as airflow splits to different extracting panels.
In this case, the cross-section did not have a substantial variation along the intake drift. Therefore, the heat contribution was similar along the intake drifts and the airflow was mainly reduced due to its distribution to the workplace. Distance was also very important to consider because of the interaction between the air and heat sources along the drifts. Equation (5) describes the factors influencing temperature variation along a drift.
where C is the airflow (m 3 /s), and D is the distance (m). A mathematical model is a simplified representation of reality with predictive variables (x1, x2, …, xn) and target variables that need to be estimated [38]. The model presented herein was based on the variables detailed in the previous paragraph, with the aim of forecasting the temperature using in situ measurements of temperature and airflow along the ventilation circuit, taking into account seven preset control points (see Figure 2). An analysis of the real data suggested a strong linear relationship among the effective temperatures in the nodes, especially in the case of adjacent nodes, as can be seen in Figure 3. Here, t n is the temperature, whether measured or estimated ( • C), c n is the airflow measured (m 3 /s), x n is the relative position of the control points or nodes (m), and n i is the node number.
An analysis of the real data suggested a strong linear relationship among the effective temperatures in the nodes, especially in the case of adjacent nodes, as can be seen in Figure 3.

Model Used
On the basis of data analysis, a generalized linear model with four parameters (see Equation (6)) and a linear model with interactions and seven parameters (see Equation (7)) were chosen for this study. A sensitivity analysis was carried out to assess the significance of the parameters from the model [39].
= 0 + 1 + 2 + 3 + 4 + 5 + 6 , where ti is the temperature at the beginning of the stretch (°C), tf is the temperature at the end of the stretch (°C), Dx is the length of the stretch (m), and Dc is the airflow change along the stretch (m 3 /s). The model estimated the daily mean temperature at the end of the stretch as a function of the daily mean temperature at the beginning of the stretch, the length of the stretch, and airflow variation. The prediction was achieved using information from node to node, as well as with a multi-node approach. Despite the quantity of data being very large, it was necessary to obtain all parameters in all nodes to use the models proposed in Equations (6) and (7). The final amount of data could be substantially reduced due to various causes, such as sensors being broken or lost due to aggressive operation conditions, abnormal values due to heavy equipment working close to the sensor, or sensor malfunction. The use of a multi-node approach with a lag step allows obtaining more data. In this study, the following lag steps were applied: -Lags step of one node, using six stretches created by nodes 1, 2, 3, 4, 5, 6, and 7; -Lag step of two nodes, using five stretches created by nodes 1-3, 2-4, 3-5, 4-6, and 5-7; -Lag step of three nodes, using four stretches created by nodes 1-4, 2-5, 3-6, and 4-7.
In addition, the multi-node system reduces the tendency toward the prediction of mean values. Figure 4 displays the variables used in the model, while Table 1 details the length, Dx, and the airflow variation along the intake, Dc.

Model Used
On the basis of data analysis, a generalized linear model with four parameters (see Equation (6)) and a linear model with interactions and seven parameters (see Equation (7)) were chosen for this study. A sensitivity analysis was carried out to assess the significance of the parameters from the model [39].
where t i is the temperature at the beginning of the stretch ( • C), t f is the temperature at the end of the stretch ( • C), D x is the length of the stretch (m), and D c is the airflow change along the stretch (m 3 /s). The model estimated the daily mean temperature at the end of the stretch as a function of the daily mean temperature at the beginning of the stretch, the length of the stretch, and airflow variation. The prediction was achieved using information from node to node, as well as with a multi-node approach. Despite the quantity of data being very large, it was necessary to obtain all parameters in all nodes to use the models proposed in Equations (6) and (7). The final amount of data could be substantially reduced due to various causes, such as sensors being broken or lost due to aggressive operation conditions, abnormal values due to heavy equipment working close to the sensor, or sensor malfunction. The use of a multi-node approach with a lag step allows obtaining more data. In this study, the following lag steps were applied: -Lags step of one node, using six stretches created by nodes 1, 2, 3, 4, 5, 6, and 7; -Lag step of two nodes, using five stretches created by nodes 1-3, 2-4, 3-5, 4-6, and 5-7; -Lag step of three nodes, using four stretches created by nodes 1-4, 2-5, 3-6, and 4-7.
In addition, the multi-node system reduces the tendency toward the prediction of mean values. Figure 4 displays the variables used in the model, while Table 1 details the length, Dx, and the airflow variation along the intake, Dc. Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 11

Stretch Nodes Dx (m) Dc (m 3 /s)
n6-n7 530 −22 Airflow variation was caused by its allocation depending on mine planning and air requirements in different workplaces, which could vary over time. Moreover, stretches S3 and S4 had a small reduction in airflow due to leakages in the circuit.

Results and Discussion
The generalized linear model (Equation (6)) and the model with interactions (Equation (7)) exhibited correlations of R 2 = 0.858 and 0.887, respectively. However 5.96% of the data recorded by the sensors displayed irregular behavior compared to the dataset, i.e., extremely high temperatures (greater than ±1.5 IQR (interquartile range)), probably due to the exhaust pipe of the diesel equipment passing close to the sensor during some specific maneuvers in the drift. The data analysis also suggested that some outliers could have been caused by erroneous sensor recording [38]. These outliers were avoided for the creation of the model.
The model fully developed with only regular data is displayed in Equation (8). The R 2 of the model with interactions was higher than that of the generalized linear model; thus, it was used moving forward, attaining R 2 = 0.933 and an estimation error of ±2 °C, with 100% of the residual values lying in the range of ±2.5 °C (Figure 5), which can be considered acceptable for mine ventilation analysis according to the literature [27]. Furthermore, it can be modified according to the needs of the user and characteristics of the mine.   Airflow variation was caused by its allocation depending on mine planning and air requirements in different workplaces, which could vary over time. Moreover, stretches S3 and S4 had a small reduction in airflow due to leakages in the circuit.

Results and Discussion
The generalized linear model (Equation (6)) and the model with interactions (Equation (7)) exhibited correlations of R 2 = 0.858 and 0.887, respectively. However 5.96% of the data recorded by the sensors displayed irregular behavior compared to the dataset, i.e., extremely high temperatures (greater than ±1.5 IQR (interquartile range)), probably due to the exhaust pipe of the diesel equipment passing close to the sensor during some specific maneuvers in the drift. The data analysis also suggested that some outliers could have been caused by erroneous sensor recording [38]. These outliers were avoided for the creation of the model.
The model fully developed with only regular data is displayed in Equation (8). The R 2 of the model with interactions was higher than that of the generalized linear model; thus, it was used moving forward, attaining R 2 = 0.933 and an estimation error of ±2 • C, with 100% of the residual values lying in the range of ±2.5 • C ( Figure 5), which can be considered acceptable for mine ventilation analysis according to the literature [27]. Furthermore, it can be modified according to the needs of the user and characteristics of the mine.  Figure 6 displays the stretches used in the analysis of the estimation error using the node-tonode and multi-node approaches. The x-axis includes the combination of nodes used in the multinode approach, as detailed in Section 3.3.
The highest error was concentrated in nodes 1 and 2, when the lag step was 1, 2 or 3. This was due to the ventilation circuit of the case study (a U-shaped system) having a high variability in temperature at the entrance of the intake, as well as potential problems with connections between the return and intake drifts through doors, and other disturbing factors in areas with a high pressure difference. This issue coincided mainly with nodes 1 and 2, which were within an area with a high airflow and pressure difference. If only nodes 3-6 were considered, the relative error was reduced to <10%. Figure 6. Error distribution as a function of stretch, node-to-node (red), and multi-node lag steps (blue represents two steps and gray represents three steps).
When the estimated temperatures were compared to the observed values in Figure 7 (i.e., those at the end of the stretch (tf)), it can be seen that the estimation was quite good, with a deviation lower than 2 °C in most of the cases. Furthermore, the correlation was also good, whether using node-tonode or multi-node analysis with a lag step of 2 and 3 nodes (see Figure 8); there was important consistency between the temperatures measured and those estimated at the end of the stretches. The  Figure 6 displays the stretches used in the analysis of the estimation error using the node-to-node and multi-node approaches. The x-axis includes the combination of nodes used in the multi-node approach, as detailed in Section 3.3.  Figure 6 displays the stretches used in the analysis of the estimation error using the node-tonode and multi-node approaches. The x-axis includes the combination of nodes used in the multinode approach, as detailed in Section 3.3.
The highest error was concentrated in nodes 1 and 2, when the lag step was 1, 2 or 3. This was due to the ventilation circuit of the case study (a U-shaped system) having a high variability in temperature at the entrance of the intake, as well as potential problems with connections between the return and intake drifts through doors, and other disturbing factors in areas with a high pressure difference. This issue coincided mainly with nodes 1 and 2, which were within an area with a high airflow and pressure difference. If only nodes 3-6 were considered, the relative error was reduced to <10%. Figure 6. Error distribution as a function of stretch, node-to-node (red), and multi-node lag steps (blue represents two steps and gray represents three steps).
When the estimated temperatures were compared to the observed values in Figure 7 (i.e., those at the end of the stretch (tf)), it can be seen that the estimation was quite good, with a deviation lower than 2 °C in most of the cases. Furthermore, the correlation was also good, whether using node-tonode or multi-node analysis with a lag step of 2 and 3 nodes (see Figure 8); there was important consistency between the temperatures measured and those estimated at the end of the stretches. The Figure 6. Error distribution as a function of stretch, node-to-node (red), and multi-node lag steps (blue represents two steps and gray represents three steps).
The highest error was concentrated in nodes 1 and 2, when the lag step was 1, 2 or 3. This was due to the ventilation circuit of the case study (a U-shaped system) having a high variability in temperature at the entrance of the intake, as well as potential problems with connections between the return and intake drifts through doors, and other disturbing factors in areas with a high pressure difference. This issue coincided mainly with nodes 1 and 2, which were within an area with a high airflow and pressure difference. If only nodes 3-6 were considered, the relative error was reduced to <10%.
When the estimated temperatures were compared to the observed values in Figure 7 (i.e., those at the end of the stretch (t f )), it can be seen that the estimation was quite good, with a deviation lower than 2 • C in most of the cases. Furthermore, the correlation was also good, whether using node-to-node or multi-node analysis with a lag step of 2 and 3 nodes (see Figure 8); there was important consistency between the temperatures measured and those estimated at the end of the stretches. The temperature oscillations corresponded to the different seasons of the year, influencing the effective temperature at the bottom of the shaft and along the intake.   Final temperature comparison for different node lag steps as a function of node-to-node (red) and multi-node lag steps (blue represents two steps and gray represents three steps).
If the analysis was carried out with the intermediate or final nodes, the correlation between observed and estimated values increased considerably (see Figure 9). Additionally, there was a substantial reduction in terms of range of deviation (~50%) compared to the analysis of the whole system detailed in Figure 7. This means that the model can be used to predict the temperature in a possible extension of the ventilation circuit, even when only using data from a part of the intake drift, by determining the effective temperature at the entrance of each extracting panel of the mine.   Final temperature comparison for different node lag steps as a function of node-to-node (red) and multi-node lag steps (blue represents two steps and gray represents three steps).
If the analysis was carried out with the intermediate or final nodes, the correlation between observed and estimated values increased considerably (see Figure 9). Additionally, there was a substantial reduction in terms of range of deviation (~50%) compared to the analysis of the whole system detailed in Figure 7. This means that the model can be used to predict the temperature in a possible extension of the ventilation circuit, even when only using data from a part of the intake drift, by determining the effective temperature at the entrance of each extracting panel of the mine. Final temperature comparison for different node lag steps as a function of node-to-node (red) and multi-node lag steps (blue represents two steps and gray represents three steps).
If the analysis was carried out with the intermediate or final nodes, the correlation between observed and estimated values increased considerably (see Figure 9). Additionally, there was a substantial reduction in terms of range of deviation (~50%) compared to the analysis of the whole system detailed in Figure 7. This means that the model can be used to predict the temperature in a possible extension of the ventilation circuit, even when only using data from a part of the intake drift, by determining the effective temperature at the entrance of each extracting panel of the mine. Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 11 Figure 9. Analysis of temperature correlation and deviation between nodes 6 and 7.
Further research should be done to analyze different ventilation circuits and additional variables related to the underground mining environment, such as pollutants or humidity. The prediction of temperatures in other sectors was studied using machine learning techniques [40]. Other models and approaches should also be explored in the future.

Conclusions
A model with interactions and seven parameters was found to be the most appropriate to describe the main ventilation system in an underground mine, showing good correlation between estimated and observed temperatures (R 2 = 0.933), with an estimation error of ±2 °C. Furthermore, its soundness was proven in both node-to-node and multi-node analyses. The model obtained can be useful to predict temperature evolution along the principal circuit, find the best workplace conditions in terms of temperature, and analyze different mine planning scenarios. It is simple to use, adaptable to further modifications depending on the data and requirements of the mine, and accurate in temperature predictions.
The higher error seen in nodes 1 and 2 suggests that the analysis and model of a U-shaped ventilation system should begin at least several hundred meters away from the bottom of the shafts or ramp (>2 km in this case study) in order to obtain the highest accuracy (as seen in nodes 3-7); otherwise, the range of the estimation error may be considerable due to the influence of potential disturbances in areas with a high pressure differential and large amounts of airflow. The multi-node approach is an interesting option in underground mining, due to there being many factors that can affect the availability of data from the sensors. The placement of sensors was also found to be crucial for obtaining representative data and, subsequently, an accurate model. It is important to find zones of drift with the lowest possible disturbance, especially regarding the operation of diesel equipment, due to the high temperature of exhaust fumes. In addition, it is necessary to have an accurate calibration protocol to avoid outliers.