Thermodynamic-Based Rock Temperature Response and Geological Structure Evolution Simulation

Figure 1 shows the principle of rock temperature response under the influence of the surface heat flow field. Surface heat flow is not only affected by the geothermal gradient but is also closely related to underground heat source activities and crustal structural changes. In traditional static grid methods, since the computational grid remains unchanged between different time steps, the complex variations in the surface heat flow field during the geological structure evolution process are often overlooked. In practical applications, especially when influenced by solar radiation or crustal movement, the computational grid tends to change. For example, the movement of the solar radiation center can lead to grid updates, causing the grids at different times to be inconsistent. In the Mckenzie model, if the computational grids at two time steps are different, the heat flow information generated at the previous time step may not be accurately recognized by the new grid at the subsequent time step, which in turn affects the accuracy of temperature response prediction. Therefore, using the moving grid method can dynamically update the computational grid and capture changes in the surface heat flow in real-time, providing more accurate and timely heat flow data for subsequent temperature response predictions.

By continuously updating the grid positions and distribution according to the changes in geological structure, it ensures that each computational grid is adapted to the current geological background. This dynamic grid can flexibly respond to surface heat flow fluctuations caused by geological activities, capturing the processes of heat flow propagation, energy generation, and dissipation in a timely manner. It is especially capable of reflecting changes in heat sources and how heat propagates in different directions, particularly in complex tectonic belts and volcanic activity regions. When the grid is updated, the heat flow information from the previous time step may not be directly compatible with the new grid, so it is necessary to transfer and adjust the heat flow data based on the relationship between the two grids. This process involves topological analysis of grids from different time steps, and through interpolation, the heat flow data from the previous time step is adapted to the current grid to ensure the continuity and accuracy of the temperature response.

1.png

Figure 1. Principle of rock temperature response under the influence of surface heat flow field

To achieve real-time monitoring of surface heat flow propagation and energy generation and dissipation processes, this paper combines the Mckenzie model with a spectral interpolation numerical method to ensure the accurate transfer and update of surface heat flow information between different time steps, thus providing effective support for real-time monitoring of rock temperature response. In the Mckenzie model, the surface heat flow field is described by the two-dimensional distribution of the spectral density V(δ, ϕ). This spectral density not only describes the energy distribution of surface heat flow at different frequencies and directions but also reflects the impact of dynamic geological structure changes on heat flow propagation. Therefore, accurately updating the spectral energy density information generated at the previous time step to the new grid nodes at the next time step is crucial to ensuring the accuracy of real-time monitoring.

In the specific numerical method, the spectral density data of the first set of grid nodes at time S needs to be determined. For each grid node, the distribution of the spectral density at different frequencies and directions has already been effectively described in the calculation results at time S. Then, at time S+1, by comparing the new grid node positions with the geographical locations of the grid nodes at time S, we can select the two nearest grid nodes. The spectral density data of these two grid nodes will be used to estimate the spectral density of the new grid node. Specifically, the average frequency and direction angle of the two selected grid nodes can help calculate the spectral density of the new node. Assuming the grid node is numbered by u and the spectral frequency time step is represented by j, the average frequency of these points can be obtained by the following formula:

$l_{u, j}=\int V_u(\delta, \phi) \delta^j f \delta f \phi$           (1)

To calculate the spectral density of the new node, the average frequency of these two grid nodes must first be obtained. This is done by weighted averaging the frequencies of the two nodes, with the weights typically determined by the distance relationship between the nodes. With the calculated average frequency and direction angle, the spectral density of the new grid node can be determined. Specifically, the average frequency of node i can be obtained by the following formula:

$\bar{\delta}_u=\frac{l_{u, s+1}}{l_{u, s}}$         (2)

Assuming that the geographical distance from node 1 to the grid point is represented by q₁, and the geographical distance from node 2 to the grid point is represented by q₂, the average frequency of the new grid node at time S+1 is:

$\bar{\delta}=\left(q_2 l_{1, s+1}+q_1 l_{2, s+1}\right) /\left(q_2 l_{1, s}+q_1 l_{2, s}\right)$          (3)

The average direction angle of these points can be obtained by the following formula:

$\begin{aligned} & l_{u, a}=\int V_u(\delta, \phi) COS(\phi) f \mathscr{f} \phi \\ & l_{u, b}=\int V_u(\delta, \phi) SIN(\phi) f \mathscr{\mathscr { f }} \phi\end{aligned}$           (4)

The average direction angle of node u at time S can be obtained by the following formula:

$\bar{\phi}_u=x T A N\left(\frac{l_{n, b}}{l_{u, a}}\right)$         (5)

The average direction angle of the sought point can be obtained by the following formula:

$\bar{\phi}=x T A N\left(\frac{q_2 l_{1, b}+q_1 l_{2, b}}{q_2 l_{1, a}+q_1 l_{2, a}}\right)$           (6)

Based on the calculations of the node's average frequency and average direction angle, the spectral density of the sought point can be further determined:

$\begin{aligned} & V(\delta, \varphi)=q_2 V_1\left[\bar{\delta}_1 \delta / \bar{\delta}, \varphi-\left(\bar{\varphi}-\bar{\varphi}_1\right)\right]+q_1 V_2\left[\bar{\delta}_2 \delta / \bar{\delta}, \varphi-\left(\bar{\varphi}-\bar{\varphi}_2\right)\right]\end{aligned}$             (7)

Next, spectral density interpolation is performed using these data. The goal of spectral interpolation is to accurately map the spectral density matrix at time step S to the new grid nodes at time step S+1 through topological relationships. Since the grids at time steps S and S+1 may be different, traditional interpolation methods may not be directly applicable. Therefore, this paper adopts a distance-weighted interpolation method, combining the geographical position relationships and spectral density characteristics of the grid nodes before and after to perform the interpolation process. The interpolated spectral density matrix not only effectively updates the node information but also maintains the consistency and continuity of spectral energy in both frequency and direction.

The result of spectral density interpolation is a new two-dimensional V×L matrix, where V represents the frequency segments and L represents the direction segments. This matrix provides accurate surface heat flow information for the new grid nodes at time step S+1. At this point, the surface heat flow field of the grid nodes at S+1 fully reflects the energy distribution from the previous time step while adapting to the new grid's geographical location and structural characteristics. Therefore, spectral density interpolation not only solves the information transfer issue caused by grid updates but also provides reliable heat flow data for subsequent temperature response predictions.

Real-time monitoring of rock temperature response requires efficient and accurate data processing methods. Considering that the rock mountain area is a large model area, the generated grid files contain a large number of node data, and traditional numerical methods may face problems such as large computational load and low computational efficiency when handling large-scale data. In order to achieve real-time monitoring of rock temperature responses, this paper chose to interpolate the spectral energy density matrices between the two grids, which can accurately transmit the heat flow field information from the previous time step to the next time step while ensuring efficiency. The interpolation method utilizes the relationship between the grids at the previous and next time steps to effectively estimate the spectral density distribution of the new grid nodes, avoiding data inconsistency caused by grid changes.

In order to achieve real-time monitoring of rock temperature response, three main interpolation methods were considered in this paper. The nearest point assignment method directly assigns the heat flow information of the new grid node to the heat flow data of the nearest old grid node. This method is simple, requires little computation, and can quickly obtain surface heat flow information for new nodes, making it suitable for scenarios where high computational speed is required. However, the nearest point assignment method is too simple and does not fully consider the spatial relationship between the new node and multiple nodes, so it may lead to large interpolation errors, especially in areas with complex geological structures or drastic heat flow changes. The linear interpolation method uses inverse distance weighting, where nearby nodes have a greater influence on the interpolation result. Specifically, using inverse distance weighting, nearby nodes have a larger effect on the result, thus more accurately reflecting local heat flow changes. To improve computational efficiency, this paper adopts a linear interpolation strategy based on the "window method." In this method, when the program is initially compiled, it filters all nodes and eliminates those far from the interpolation node, retaining only the nearby nodes for interpolation calculations.

The triangular area method is a geometry-based interpolation method, particularly suitable for handling grid nodes that are triangular elements. This method not only considers spatial topological relationships but also accurately reflects local variations of triangular grids. This method is suitable when there are significant changes in grid topology, especially in cases where high precision is needed to describe the heat flow distribution of new nodes, thus effectively improving the accuracy of interpolation results. Specifically, assume the energy density of the new node is represented by V_j, and the energy densities of the existing unit nodes X, Y, Z are represented by V_x, V_y, V_z, and the areas of the triangular elements are T_x, T_y, T_z. The surface heat flow information for the newly inserted node can be obtained as follows:

$V_j=V_x T_x+V_y T_y+V_z T_z$          (8)

To achieve accurate rock temperature response prediction, this paper first uses the inverse distance weighted average method for surface heat flow information topology. This method can effectively capture local variations in the heat flow field while maintaining high computational efficiency. By introducing the "window method" in the interpolation process, nodes close to the new node are selected for weighted averaging, reducing computational load and thus improving computational speed. To further improve the model's operational efficiency, this paper introduces parallel computing technology, using multi-threading to process grid node information. In this way, although single-threaded computation may take longer, parallel computing can control topology time within 2 minutes, greatly improving overall computational efficiency.

Based on the above efficient computational methods, once accurate surface heat flow information is obtained in real time, it can be input into the temperature response model, and temperature changes at each node can be calculated through numerical simulation and physical equations. Specifically, using the heat conduction equation and considering boundary and initial conditions, the evolution process of the internal temperature of the rock with time and space can be simulated. In this process, the generation, propagation, and dissipation of heat flow are accurately considered and reflected, ensuring the accuracy of the simulation results. By continuously updating surface heat flow information and performing real-time temperature response calculations, this paper can dynamically monitor and predict temperature changes in rocks, promptly capturing the impact of geological events or structural activities on the temperature field, thereby achieving precise prediction of rock temperature response.

From Figure 5, it can be observed that the average temperature in the rock model changes with axial strain in a certain regularity under different thermal expansion coefficients. As the axial strain increases, the temperature of the rock gradually rises, especially in the high-strain zone (such as when the strain exceeds 0.02), where the temperature increase becomes more significant. For example, when the axial strain is 0, the average temperature of the rock is 0, and as the strain gradually increases, the temperature also shows a gradual rise. Under different thermal expansion coefficients, the temperature response of the rock exhibits different rates and amplitudes of change, especially at the thermal expansion coefficients of 0.0455×10^-6 and 0.0551×10^-5, where the amplitude of temperature change with strain is relatively large. This change can be attributed to the influence of the thermal expansion coefficient on the accumulation and conduction of heat inside the rock. The larger the thermal expansion coefficient, the more significant the increase in internal temperature under strain. Furthermore, as axial strain increases, the rise in temperature also reflects the stress and thermodynamic changes the rock undergoes during the geological structure evolution process, especially when the strain reaches a critical value, after which the temperature rise becomes more gradual, indicating that the rock structure begins to stabilize.

From Figure 6, it can be observed that the range of the rock model changes with axial strain under different thermal expansion coefficients shows an obvious trend. Under low thermal expansion coefficients, as the axial strain increases, the range gradually increases from 0, and the increase is relatively smooth, reaching a maximum of 0.45. This indicates that under low thermal expansion conditions, the temperature response of the rock is more moderate, and the temperature change is less sensitive to strain. In contrast, under higher thermal expansion coefficients, as axial strain increases, the range rises significantly, especially when the strain exceeds 0.6, the range quickly jumps to 1.6, and continues to increase at a high rate as strain continues to increase. This indicates that a larger thermal expansion coefficient causes a greater temperature change in the rock under strain, and the thermal response of the rock becomes more sensitive, showing more significant changes in temperature range. The rock models with different thermal expansion coefficients exhibit a strong nonlinear relationship between temperature range and axial strain, especially in the larger strain range, where the temperature response is more pronounced.

5.png

Figure 5. Change of average temperature with axial strain in the rock model under different thermal expansion coefficients

6.png

Figure 6. Change of range with axial strain in the rock model under different thermal expansion coefficients

7.png

Figure 7. Temperature and depth relationship curve of the rock mountain under steady-state heat conduction

Based on the data provided in Figure 7 for the temperature and depth relationship curve under steady-state heat conduction conditions, it can be clearly observed that with the increase in depth, the temperature of the rock gradually increases. This trend has been verified under different temperature settings. Specifically, under the condition of 120℃, the temperature of the rock increases from 20℃ to 270℃ with increasing depth, and the temperature gradient remains relatively steady, with the largest change occurring in the depth range of 1.6 meters to 1.8 meters, showing a strong correlation between depth and temperature. Under conditions of 200℃ and 260℃, the temperature of the rock exhibits a more obvious warming trend, especially at deeper depths (such as 1.4 meters and above), where the temperature change is more pronounced, indicating that the intensity and depth of the heat source have an important impact on the temperature distribution. From the trendline data, it can be seen that the rock temperature increases more significantly under higher temperature settings, and at deeper depths, the temperature rise is more substantial, particularly when the heat source temperature is high, leading to a significant increase in the temperature gradient.

Based on the above experimental data and the relationship between temperature and depth, it can be inferred that there may be regions in the geological structure where the heat source is more concentrated, especially in the depth range of 1.2 meters to 1.6 meters, where temperature changes are more dramatic. This may indicate uneven distribution of local heat sources, leading to a significant nonlinear change in the temperature field. This phenomenon is consistent with the geological occurrence of uneven heat source distribution and the accumulation of heat at depth. Furthermore, the change in temperature gradient suggests that the structural evolution of the region may include compression, folding, or faulting activities in the rock layers, which can influence the distribution of heat and its heating effect on the rock, further promoting changes in the temperature gradient.

According to the inversion experiment results of the upper boundary layer of Rock Mountain provided in Table 1, the temperature variation trends at the boundary layer under active and passive temperature fields can be observed. Specifically, when the heat source intensities are 115, 216, and 257, there are certain changes in the surface temperature and upper boundary layer temperature under the active temperature field. When the heat source intensity is 115, the surface temperature is 25.6℃, the passive boundary layer temperature is 16.5℃, and the active upper boundary layer temperature is 18.9℃. As the heat source intensity increases to 216, the surface temperature rises to 28℃, and the boundary layer temperature and upper boundary layer temperature are 16.5℃ and 22.6℃, respectively. When the heat source intensity reaches 257, the surface temperature slightly increases to 28.4℃, while the boundary layer temperature remains at 16.5℃, and the active upper boundary layer temperature is 21.5℃. These data show that with the increase of the heat source intensity, the surface temperature and upper boundary layer temperature in the active temperature field gradually increase, while the passive boundary layer temperature remains relatively constant, showing a clear influence of the heat source on the upper boundary layer. From the experimental results, it can be concluded that the heat source intensity significantly impacts the temperature field of the Rock Mountain boundary layer. As the heat source intensity increases, the surface and upper boundary layer temperatures gradually increase, especially in the active temperature field, where the temperature variation is more obvious. This indicates that the heat source intensity is directly positively correlated with the rock temperature variation. The stability of the passive boundary layer temperature might suggest a relatively stable geological structure background, where the heat source does not directly affect the deep rock structure or the heat conduction process is somewhat suppressed.

Table 1. Rock mountain upper boundary layer inversion experiment results

According to the Rock Mountain heat source temperature inversion experiment results provided in Table 2, the matching of the original heat source intensity and the inferred heat source temperature can be observed. For a heat source intensity of 115, the original value is 115, the inferred value is 124.23, with an error of 0.82, and the relative error is only 0.67%. When the heat source intensity is 216, the original value is 216, the inferred value is 213.26, the error is 1.56, and the relative error is 0.81%. For a heat source intensity of 257, the original value is 257, the inferred value is 256.89, the error is 6.48, and the relative error is 2.56%. From these results, it can be seen that the heat source temperature inversion model can accurately invert the heat source temperature in most cases, especially under lower heat source intensities, where the relative error is small. As the heat source intensity increases, the relative error of the inversion result slightly increases, particularly for higher heat source intensities (257), where the error increases, but it remains within a reasonable range.

Table 2. Rock mountain heat source temperature inversion experiment results

From the heat source temperature inversion experiment results, it can be seen that the model demonstrates high inversion accuracy under most heat source intensities. When the relative error is small, the inversion method is accurate and can restore the actual heat source temperature well, indicating that the heat source inversion method has high reliability in predicting temperature fields. However, at higher heat source intensities, the relative error increases, possibly due to the non-uniformity of the heat source at higher temperature gradients and the limitations of the model's calculations in high-temperature regions. This result suggests that, for larger heat source intensities, there may be more complex geological structures or heat source distributions, which require further optimization of the model's adaptability.

Based on the heat source temperature inversion error data provided in Figure 8, the variation trend between the inferred values and the true values can be observed. In the heat source temperature range from 120℃ to 263.5℃, the inferred values are very close to the true values, and the difference between them is very small. Specifically, at different heat source temperatures, the difference between the inferred values and true values is relatively small, and the error value gradually increases, but remains within 1℃. For example, when the true value is 120℃, the inferred value is 120℃, with an error of 0℃; when the temperature increases to 140.5℃, the inferred value is 141℃, with an error of 0.5℃; at the maximum value of 263.5℃, the inferred value is 267℃, with an error of 3.5℃. Overall, the inversion error increases slightly with the temperature, but these errors are relatively small, indicating that the heat source temperature inversion method can predict the heat source temperature well in most cases.

8.png

Figure 8. Heat source temperature inversion error chart