{'fraction_non_alphanumeric': 0.030355908207947847, 'fraction_numerical': 0.014947486254238963, 'mean_word_length': 4.782219683990101, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical Model to Assess the State and Increase of Temperatures in Underground Mine Galleries: A Tool to Support Heat Recovery Projects

Underground mining is facing growing challenges related to the need to mine deeper and at higher temperatures, to operational expenditures associated with energy consumption, lower grade ores, environmental constraints, and social pressures. In this scenario, a new numerical model is proposed to estimate temperature increase inside mining galleries to provide specific criteria for heat recovery projects, which may consider heat extraction from abandoned mines using closed-loop geothermal systems or from operating mines using the exhaust ventilation air. This model couples different approaches from previous models and include key parameters unemployed until this moment, such as wall roughness and velocity profile modeling, what would allow for a more realistic estimation of convective heat transfer phenomena, which is critical to predicting heat exchange in ventilation air due to the turbulent nature of the airflow. The model also includes other heat sources that could be present inside galleries and should be accounted for, such as machinery, once the heat dissipated to the environment might be substantial depending on the equipment and gallery geometry. The general intention of this project is to account for every heat source that may contribute to increasing the temperature inside the gallery, so it becomes tangible to harness as much heat energy as possible, preventing energetic losses and stimulating an increase of thermodynamic efficiency in underground mining operations. The model is not validated yet with real temperature data, but preliminary results agree with the ones from previous models.


Introduction
Mining has been facing growing pressures of environmental, economic, social, and technical nature in the last years. Full legal compliance with state regulations has become an increasingly insufficient means of satisfying society's expectations concerning mining issues (Prno, 2013). In this scenario, emerges the need to obtain a social license to operate, strongly linked to society's perception of legitimacy.
One of the possible perspectives regarding legitimacy is the moral field, which depends on whether an action is viewed as acceptable by the stakeholders. Procedural legitimacy, for instance, emphasizes the means of achieving a particular goal, and not the outcome itself: it is expected that the company adopts the latest technologies or processes in compliance with codes of best practice, to demonstrate it adheres to adequate methods of production (Saenz, 2019).
One of the many actions that may be considered to suit the public expectation is to increase efficiency and reduce waste. The roadmap for transforming the EU into a competitive, lowcarbon economy by 2050, highlights a potential for CO2 emission reduction through an increase in energetic efficiency (European Commission, 2018), that would result in less energy consumption and lower environmental impacts related to power generation.
From an energetic point of view, the mining industry is exceptionally electro-intensive. In the most prominent mining countries, power consumption by the mining industry alone usually represents considerable fractions of the total power generated, as in the United States (3%) and South Africa (6%) (Carvalho and Millar, 2012). Worldwide, mining consumes approximately 4% of all the power generated (International Energy Agency 2017).
Mining's strong power dependency is also demonstrated by its impacts on the operational costs of a mine, which ranges International Journal of Terrestrial Heat Flow and Applied Geothermics. VOL. 4, NO. 1 (2021); P. 110-118. from 15 to 40% of a mine's OPEX (ARENA, 2018). Even though it is necessary to reduce energy consumption, the tendency is that global demand increases in the next years. The increase in the demand for energy in mining can be caused by numerous factors, including the reduction of ore grades and the growing difficulties associated with accessing new orebodies, impacting the costs (Moran et al., 2014).
Nevertheless, the technical viability of new mining projects is a core issue, especially when it comes to underground mining, where the production tends to occur at greater depths (Thrybom et al., 2015). Deep mining happens in a challenging context, where innovative solutions and the development of new practices are required, and also the adoption of additional safety measures to overcome difficulties and ensure adequate working conditions and the expected profits (Ranjith et al., 2017).
To address these challenges, a numerical heat transfer model was developed to support an underground mining operation's energy efficiency through a better understanding of the thermal energy flows inside the mining galleries and how these flows relate to ventilation operations. The model performs calculations to estimate heat exchange inside an underground mining gallery, which allows for an approximation of the temperature increase of ventilation air along its course through the galleries network. This temperature increase is the critical parameter for optimizing the airflow and evaluating the possibility of recovering part of the thermal energy from the exhaust air. The model could also be applied to estimate heat recovery in closed-loop geothermal systems and to provide criteria for viability studies.

Methodology
The use of geothermal technologies to recover thermal energy from underground environments is well established outside the mining industry. However, until the present day, thermal energy has been overlooked inside the sector, with less than 20 documented cases of geothermal systems in mines around the world until the year 2013 (Preene and Younger, 2014). Among the documented cases, some projects have been developed to recover geothermal energy from rock masses and can be divided into two categories: I) projects to recover heat from abandoned galleries after the mine's closure; II) projects to recover heat simultaneously to the mine's operation.
In both cases, many different viability studies and mathematical models have been developed to plan and support the implementation of such projects, allowing for the estimations of heat exchange, of the temperature increase of the fluid (water or air) flowing through the galleries, or even long term exhaustion of the geothermal reservoir. Thus, a literature review is essential to formulate a numerical model, so it´s possible to assess what has been accomplished and what is yet to be done in heat modeling.

Projects of heat recovery after the mine is closed
In this category are the projects that are based in the flooding of abandoned galleries, that act as a heat exchanger heating water that flows at a particular rate trough the networks, in a closed-loop geothermal system.
A well-documented case of such kind is the one reported by Jessop et al. (1995), developed to recover heat from an abandoned coal mine in Nova Scotia, Canada, where the heated water is used by a town located in the vicinities of the mine to heat buildings and houses. Banks et al. (2004) also report other operations where closed-loop geothermal systems are applied in abandoned mines in the United States, Norway, and Scotland.
The main issue regarding the implementation of these projects is to demonstrate the heat-recovery potential, which is based on the study of analog situations, empirical models, and mathematical models. Rodríguez and Díaz (2009) proposed a semi-empirical model to estimate the temperature increase in the water flowing through galleries, and Ghoreishi-Madiseh et al. (2012) developed a mathematical model based on transient heat exchange differential equations to estimate both the heat increase in the water and the depletion of thermal energy in the rock mass. Those models might be an essential means to make feasible more projects involving heat recovery in underground mines.

Projects to recover heat simultaneously to the mine operation
During the lifetime of a mining operation, with the presence of machinery, personal, electrical equipment, and so forth, one alternative to recover heat from the rock mass is to use the ventilation air as a sort of geothermal fluid.
Some projects of this type are based on the recirculation of ventilation air that has already been heated up during its course through the galleries and can be used as a way to increase the intake air temperature and reduce pre-conditioning costs associated to the burning of fuels. That is the strategy used by a potash mine in Canada (Hall et al., 1990) and Williams mine in the Northwest of Ontario (Kalantari et al., 2019). In a different approach, Ghoreishi-Madiseh et al. (2015) proposed a model of injecting air trough broken rock mass in the bottom of an open pit, based on the seasonal thermal storage capacity of the rocks, also to save costs with the pre-conditioning of air.
Other authors refer to the possibility of using the heat in the exhaust air and not directly the heated or cooled air, as previously mentioned. Morkun et al. (2015) discuss the possibility to recover heat from the ventilation air in heat pumps, and Kalantari and Kalantari et al. (2019) created a mathematical model in MATLAB to recover heat from the exhaust air and then use this heat to increase the temperature of the ventilation's intake air in heat exchangers. Sbarba et al. (2012) developed a software to estimate economic gains from projects based on heat recovery from exhaust air. Other methods (Zhu et al., 2015) can be used to estimate the heat increase in the ventilation air along its course through the mine, supporting this kind of heat recovering project.

Numerical Model
The proposed model can be used to support projects in both the categories mentioned, to any fluid flowing through the galleries, at any flow regime. It incorporates parameters previously neglected as average gallery wall roughness, admits velocity profile variation according to the position and flow rate, includes heat exchange by convection, and deals with other punctual heat sources that may be present inside the gallery, such as machinery and personnel. All of those variables are used in an algorithm written in Octave and Python to solve a transient heat exchange differential equation in two dimensions via the Finite Volume Method.

Domain and Geometry
In order to simplify the problem's geometry, the gallery was conceived as a cylinder, wrapped by an external cylinder, representing the rock mass, so it is possible to adopt cylindrical coordinates and to use the same approach to solve the differential equations as proposed by Ghoreishi-Madiseh et al. (2012). Figure 1 shows the considered domain that admits temperature variation in the radial and longitudinal directions and assumes a uniform behavior in the angular direction. As the z-axis is a symmetry axis for the figure, it is possible to consider only the top part in the numerical model.

Physical description of the phenomena
Through the proposed simplified geometry, it is possible to subdivide the studied region into four distinct areas, as shown in Figure 2. Region A corresponds to the rock mass that revolves the gallery, except for the gallery walls; region B is the rock portion in direct contact with the fluid that flows through the galley; region C is the area inside the gallery where the fluid touches the rock wall; region D corresponds to the fluid that fills the gallery but does not touch the gallery wall. At regions C and D, the fluid flow is considered, beginning at section = ( 1 − 2 )/2 with a positive direction in z. The flow velocity is a function of the distance from the wall and depends on the flow regime. The flow is assumed incompressible, and the fluid's behavior is assumed to be the same as a Newtonian Fluid.
When the fluid moves with a laminar flow, the velocity profile inside the gallery could be assumed as in Eq. (1) (Kays, 1980): Where is the local velocity, 0 is the gallery radius, is the radial coordinate of the point, and is the mean velocity. When the flow regime is turbulent, it is considered that the internal gallery walls are rough, and the variable , representing the average roughness height, is introduced to the problem. The procedure to determine the velocity profile is based on equations (2) to (7) as described by White (2006).
First, the Prandtl Equation (2) is used to calculate Darcy's Friction Factor. Then, Eq. (3) can be written to the flow velocity: Where is Darcy's Friction Factor, is the Reynolds number for a circular duct with diameter , + is the flow velocity according to the distance to the duct wall, and = 0.41 is the Kármán constant. With , the Fanning's Friction Factor is calculated by = 4 , which relates to the shear stress at the duct wall by the equation (4): Where is the fluid's density. With the shear stress, Eq.5 allows for the calculation of the friction velocity at the duct wall * : Eq. (6) relates ̅, mean velocity in time, with + and * , and its written by: The velocity variation in time is neglected in this analysis, so ̅ = is the function for the velocity that defines the velocity profile in a turbulent flow. Replacing (3) and (5) in (6), Eq. (7) is obtained, which corresponds to the velocity profile function.
In all the domain, the phenomena of heat exchange by heat conduction occurs. Conduction is the mechanism that describes energy transportation from the region with higher temperature to the region with lower temperature by electron drift (Bayazitoglu and Ozisik, 1988), happening inside a solid or a fluid, between different solids, or between a solid and fluid as long as the fluid is not moving.
This mechanism is mathematically translated by the Fourier's Law of transient heat conduction, expressed by Eq. (8) with cylindrical coordinates, neglecting energy sources inside the solid and the energy variation in the direction ϕ.
Where ⁄ is the temperature variation with time, ⁄ e ⁄ are the temperature gradients, and = ⁄ is the thermal diffusivity of the media, being its thermal conductivity and the specific heat. The higher the value of the thermal diffusivity, the more quickly heat is diffused through the media (Holman, 2002).
The heat convection mechanism occurs at the regions B and C, once the fluid moves in contact with the rock at a different temperature. The flow is directly associated with heat transfer between solid and fluid and happens simultaneously to heat conduction (Bayazitoglu and Ozisik, 1988). This mechanism can be expressed through Eq. (9), called the Newton Law for cooling.
Where is the temperature at the solid's border, is the fluid's temperature, is the contact area between solid and fluid, and ℎ is the convection coefficient.
Towards estimating the value for ℎ, it is necessary to make a few assumptions: the wall temperature is assumed constant (meaning that the time intervals considered in the solution need to be short and the ℎ coefficient needs to be recalculated throughout the whole time interval); fluid flow is assumed fully developed; the gallery length is at least 60 times the gallery diameter. With that, it is possible to estimate ℎ by (10), according to Janna (Janna, 2000): Where is the Nusselt number. For laminar flows, = 3.66 (Bayazitoglu and Ozisik, 1988). For turbulent flows, can be estimated by the Pethukov Equation, which is an empirical formula valid for rough pipes, given by Eq. (11) for liquids and Eq. (12) for gases: Where µ is the fluid's dynamic viscosity at the fluid's temperature, µ is the fluid's dynamic viscosity at the rock's temperature, is the fluid's mean temperature, is the rock wall's mean temperature, is the Reynold's number, is the Prandtl's number, and is given by: As the fluid flows inside the gallery, a temperature gradient develops along the gallery's cross-section. Liquids are affected by temperature variation, which interferes with the fluid's physical properties, and gases are more influenced by absolute temperature (Bayazitoglu and Ozisik, 1988). The exponent of equations (11) and (12) are correction factors for the Nusselt number for initially assuming isothermal conditions for the fluid. The same correction is applied to Darcy's Friction Factor, for liquids and gases, respectively: Where is the corrected value for the friction factor and is the friction factor for isothermal conditions. The coefficients e assume values according to Bayazitoglu e Ozisik (Bayazitoglu and Ozisik, 1988 In case the mine is still operating, heat-dissipating elements may be present inside the galleries, such as machinery, auxiliary fans, and personnel. For personnel, it is assumed that the heat dissipated is 117.24W per individual (Bauer and Kohler, 2009).
For machinery, equation (18), adapted from Zhu et. al (2015), can be applied to estimate the heat dissipated to the air: Where is the heat generated from machinery in J/s; 1 is the installation factor (ratio of motor maximum consumption power to rated power), usually taken around 0.7; 2 is the simultaneous factor (ratio of summation of used motor rated power to the total rated power); 3 is the load factor, usually taken around 0.4-0.5, and is the machinery rated power, in W.
The heat dissipated by these sources cause local temperature to increase inside the gallery at the areas they are placed, and calculations are made within the model to evaluate these temperature changes in time.
In case the gallery is being ventilated with air injected from the surface, it is necessary to consider auto-compression effects associated with the air's energy gain due to the pressure increase caused by a vertical descending movement in mining shafts. Part of this energy is converted into heat, which Where is the heat transferred in the mining shafts in J/s, is the air's mass flow in kg/s, is the surface air intake temperature in Celsius, ℎ is the shaft's depth, and H is the shaft entrance elevation. It is also possible to adapt eq. (19) to account for gallery inclination, which would also cause a temperature gradient due to the auto-compression of air.

Energy balance
The heat transfer mechanisms previously discussed happen simultaneously, so it is necessary to associate its equations in a single energy balance equation, given by (20).
Where is the temperature variation in time, is the temperature variation due to the fluid flow, is the heat transferred by conduction, and is the convective heat transfer. As different domain regions may not be subjected to all the heat transfer mechanisms and effects, they have to be calculated separately. The strategy was to create a point grid, so each point is subjected to the conditions related to its position; thus, it is possible to apply the Finite Volume Method to solve the differential equation.

Finite Volume Method
Each grid point is the center of a volume cell, and its boundaries are the midpoints between the grid points. Figure  3 shows a grid point P, which is the center of the volume cell limited by the red boundaries, where the heat flow occurs. The temperature variation is assumed to be linear from one grid point to another in the directions and , as suggested by Patankar (Patankar, 1980). The velocity variation within the volume cell is considered negligible; once the volume cell's dimensions in the direction are minimal compared to the scale of the problem, the velocity function is dealt with as a constant for a particular volume cell. Then, the equation can be solved by first integrating the differential equation (20) in and : which is equivalent to The coefficients 1 to 6 from are expressed in Table 1, considering that the grid is regular in and , Equation (21) is the explicit discretization of the transient energy balance given by Eq. (20). The explicit scheme was adopted, so the calculations are simplified. If implicit methods were applied, it would be necessary to solve a set of simultaneous equations for each grid point, which would represent extra computational costs. As per Patankar (1980), if the scheme adopted is explicit and time steps are not small enough, physically unrealistic results may emerge, so it is also necessary to evaluate the maximum time interval to guarantee convergence.  ( Function is given by Eq. (1) or Eq. (7), depending on the flow regime. The harmonic mean method (Patankar, 1980) is utilized to calculate the local values for , using the thermal conductivity of the point and its neighbors', where: The same method is applied to calculate , e .

Boundary conditions
The boundary conditions defined to solve Eq. (21) were the same as proposed by Ghoreishi-Madiseh et al. (2012): an isothermal boundary is assumed at = 0, = 2 and = 2 /2, so the points located at these locations remain with constant temperature; an adiabatic frontier is assumed at = 0 due to the symmetry, so the heat flows at < 0 and > 0 are equivalent and cancel each other at the border which divides both regions.
At the instant = 0, air and rock are in thermal equilibrium, so all the points are at a temperature except for the points located at the gallery's inlet, which have temperature . At region A, 5 = 6 = 0 (there is no convection nor fluid flow); at B, 6 = 0 (there is no fluid flow); at C, all the coefficients are different from zero; at D, 5 = 0 (there is no convection).

Numerical Algorithm
An algorithm was created in GNU Octave and Python to interactively solve Eq. (21) for all the grid points after a time interval . The algorithm includes a function that automatically adjusts the time steps for the higher value possible for which the solution is convergent.
First, the input parameters shown in Table 2 are read, and all the preliminary calculations necessary to solve Eq. (21) are executed.

Parameter
Rock density ( ) Rock's specific heat ( ) Rock's thermal conductivity ( ) Rock's initial temperature ( ) Air density ( ) Air's specific heat ( ) Air's thermal conductivity ( ) Air's dynamic viscosity (µ) Prandtl's number ( ) Inlet air temperature ( ) Mean air velocity ( ) Gallery diameter ( 1 ) Gallery length ( 1 ) Outer diameter ( 2 ) Outer length ( 2 ) Mean wall roughness ( ) Number of grid points in z (ngpz) Number of grid points in r (ngpr) Each variable is assigned to all the grid points according to their position, coordinates, and media. According to its coordinates, one specific function of the algorithm determines the flow regime and the point's velocity. All boundary conditions are defined, and all coefficients from 1 to 6 are calculated beforehand, so the iterative function that solves Eq. (21) perform all calculations from pre-defined parameters, except for the convection coefficient h, that is adjusted at every iteration. Convergence is assured by a specific function verifying if any unrealistic condition emerges, forcing the calculations to stop and correcting the time step. Binary variables are attributed to differentiate rock from fluid and to distinguish the regions showed in Figure 2.
The temperature evolution contour plots are stored, and the calculations are performed until the changes in the temperature state are negligible or during a specified time interval. A table containing mean fluid temperature and its corresponding position in is created and mean exhaust temperature stored.
The temperature evolution contour plots are stored, and the calculations are performed until the changes in the temperature state are negligible or during a specified time interval. A table containing mean fluid temperature and its corresponding position in is created and mean exhaust temperature stored.

Results and discussion
The model could be used for multiple purposes. In some cases, the heat exchange process itself and the temperature variation in time might be the main object of study. In other cases, the stabilized temperature state could be the main focus, so the mean exhaust temperature and the final temperature state inside the gallery would be known, allowing for heat recovery estimating and underground environment control. Some tests regarding these different approaches are shown in the next topics.

1. Numerical process and temperature evolution
The first test assumes an empty gallery filled by air. At the instant = 0, rock and air inside the gallery are in thermal equilibrium at 31 ⁰C, and at > 0, air starts flowing at = International Journal of Terrestrial Heat Flow and Applied Geothermics. VOL. 4, NO. 1 (2021); P. 110-118. 0.1 / and 15 ⁰C. The physical properties assumed for the rock were = 1500 / 3 , = 1260 / ⁰ and = 3.1 / ⁰ , with an average roughness with an average roughness = 0.3 , equivalent to the roughness of an intake adit with medium roughness, uneven floor, and 1% area deviation (Montecinos and Wallace, 2010). Figure 4 shows the temperature evolution when cold air begins to flow until a stabilized temperature state is observed. The mean exhaust temperature is 15.75 ⁰C, showing a temperature increase of 0.75 ⁰C in the air after travelling 2000 inside the gallery. If only the average roughness was altered to = 0.1 , the mean exhaust temperature would have been 15.6 ⁰C.

2. Stabilized temperature state
The second test assumes the same conditions as the first test, but for a gallery occupied with people and machinery. Two people are located at = 200 , one at = 500 and three at = 700 ; at = 100 two YT28 rock drilling machines with a power of 0.8 kW are placed and a rake loader with a model of P-30 and power of 18.8 kW is positioned at = 500 . For all the machinery 1 = 0.7, 2 = 1 and 3 = 0.4. Figure 5 shows the stabilized temperature state for the described situation. As the equipment is placed far from the air exhaust, heat is dissipated, and heated air is mixed with colder air along the gallery length, which explains the temperature peak surrounding the machinery and the temperature drop near the gallery's end. In this specific test-case, for 500 < < 750 , the temperature is higher than 31 ⁰C at some points, suggesting that the airflow should be adjusted to dissipate the heat from the machinery faster, resulting in lower temperatures inside the gallery and consequently better environmental conditions. The mean exhaust temperature is 17 ⁰C, which represents an increase of 2 ⁰C in the air's temperature along its course through the gallery, in opposition to an increase up to 15.75 ⁰C in the first test case with the same conditions but for an empty gallery. This modeling tool could support future decisionmaking regarding optimal airflow and inlet air temperature, so more suitable conditions are provided for underground workers, and maximum heat recovery is achieved for the exhaust air.
The last test-case illustrates a typical closed-loop geothermal system, where water at 13 ⁰C is injected at 4 l/s in a gallery with a cross-section of approximately 10 2 . The physical properties assumed for the rock were = 2500 / 3 , = 1100 / ⁰ and = 2.27 / ⁰ , and the rock mass temperature is 21 ⁰C; the gallery length is 5000 m. The average wall roughness was neglected, and The Nusselt number was assumed as 1, compatible with a slug flow, once the water flow is laminar with a Reynolds number in the order of 10 2 , meaning that the convection effects occur at a lower magnitude.
These conditions are similar to those of the geothermal system installed in Springhill Nova Scotia, and the results can be compared to those obtained by Rodríguez and Díaz (2009) and Ghoreishi-Madiseh et al. (2012), for similar input conditions. Figure 6 shows the stabilized temperature state for the described situation, and Figure 7 shows a plot of the water temperature versus traveled distance inside the galleries, comparing the results from this model to the mentioned authors'. Due to the low water velocity, the flow regime is laminar, and the temperature increase in water is considerable, once it travels slower and can exchange heat with the rock walls for a more extended time. The mean temperature after running through 5000m of the gallery is 18.5 ⁰C, agreeing with the results obtained by Rodríguez and Díaz (2009) and Ghoreishi-Madiseh et al. (2012). The water's average temperature versus the gallery length is shown in Figure 7. The main divergence from the results obtained by Ghoreishi-Madiseh et al. (2012) is probably due to the convective heat transfer effects, which could cause the rock to cool faster and supply less heat to the water flowing (specially for regions close to L=0, once the gallery inlet is constantly supplied with water at 13 ⁰C) what causes the slope to be less inclined than the compared models'. Once this effect is overcome, as the heat-transfer processes happen at a faster rate due to convection, the temperature increase in water becomes superior to Ghoreishi´s model.
The differences from the results of Rodríguez and Díaz (2009) may be originated in the assumption that the rock mass temperature is constant, which would cause a more significant temperature increase in the water. Once the rock mass cools off in time, it supplies less heat to the water flowing inside the gallery, so the water will reach lower temperatures, what explains why Rodríguez and Díaz (2009) model's slope is above the two others.

Conclusions
The proposed model is a tool to support energy efficiency improvement in underground mining galleries. It brings together different approaches, such as heat recovery in closedloop geothermal systems for abandoned galleries and heat recovery from the ventilation's exhaust air in operating mines. The model also combines different models previously developed, and attempts to increase the level of detail employed in this kind of model until now, including parameters that are relevant in mining but were not considered before, such as wall roughness, and making possible to deal with heat sources inside the galleries, such as machinery and personnel.
It is, ultimately, an effort to bring thermodynamics and fluid dynamics closer to the mining reality, as it intends to adjust well-established concepts from these fields to a different environment and scale.
The results obtained are compatible with other models', but the model is yet to be tested and validated in a real case so that it can be adjusted and improved with the help of empirical parameters.