#### Abstract

In research using heat tracing technology to investigate the lateral hyporheic exchange in the shallow geological body of the riparian zone, the accurate estimation of temperature changes can provide a scientific basis for quantifying the process of lateral hyporheic exchange. To improve the accuracy of estimating temperature changes in the riparian zone, a hydrothermal coupling model considering parameter heterogeneity was established based on existing models of the relationship between thermal conductivity and saturation. The model was verified by temperature data from laboratory experiments, and the effect of the thermal conductivity prediction models was compared with that of the partial differential equation (PDE) modeling approach. The results show that the established hydrothermal coupling model can effectively characterize the temperature changes observed in a generalized laboratory model of the riparian zone, and the model simulation effects vary with the equivalent thermal conductivity models. In addition, several thermal conductivity empirical models are suggested for further application. The model parameter sensitivity analysis indicated that the hydraulic conductivity , VG model parameters (*α* and *β*) and heat capacity of soil have a relatively large effect on the temperature output of the model. The results of this study will provide reference for the selection of equivalent thermal conductivity model for simulating temperature variations in the riparian zone.

#### 1. Introduction

The riparian zone is an important transitional area for the exchange of matter, energy, and information between the river ecological system and the land ecological system [1, 2]. With yearly increases in water shortages, water pollution, and flood disaster events, the riparian zone plays a significant role in the protection of the structure and function of the aquatic ecosystem and land-and-water ecotone system [3–5]. Therefore, an increasing number of scholars have focused on work in this field, which has become a hotbed of interdisciplinary research.

Heat transport constantly occurs during the lateral hyporheic exchange process. Therefore, observations of the temporal and spatial temperature changes in the riparian zone can be used quantify the hyporheic exchange flux and even help us understand the processes of lateral hyporheic exchange occurring in the riparian zone [6, 7]. Sawyer et al. [8] evaluated the water flux and maximum distance of lateral hyporheic exchange by monitoring the water temperature and water level in the riparian zone downstream of the Longhorn dam. Gerecht et al. [9] selected Sawyer’s et al. [8] study area and established a series of water pressure and temperature sensors in the riparian zone to study the summer water level attenuation and temperature distribution trends in the riparian zone. Zhu et al. [10] selected a section of a river between the downstream tail water of Xin’an River Dam (located in Zhejiang Province, China) and Fuchunjiang Reservoir (located in Zhejiang Province, China) as the field study area; they monitored the water level and temperature in the channel and riparian zone along the river transect to quantitatively describe the spatial heterogeneity and dynamic characteristics of the riparian lateral hyporheic exchange and explained the response relationship between the water stage and temperature in the riparian zone and the channel. Ren et al. [11] quantitatively analyzed the effects of water temperature, water level, and radiation temperature (i.e., the ground surface temperature, which is influenced by solar radiation) variations on the water and heat exchange processes in the riparian zone by constructing a laboratory sand tank experiment. Watson et al. [12] continuously monitored the temperature and water level in the riparian zone along the lower Colorado River (located in Texas, USA) and studied the lateral hydrothermal exchange and flow processes under stable and unstable river water level conditions.

Many studies have used temperature data to analyze hydrothermal exchange in the riparian zone [13–16]. However, most research methods involve laboratory or field experiments, which have relatively high costs. Although a large number of models have been proposed [11, 17–19], there is still room for further improvement. Therefore, it is necessary to find a suitable coupling method. The soil equivalent thermal conductivity is an important driving factor of hydrothermal coupling models and can influence the simulation accuracy [20]. However, in some software that simulates the heat transfer of a geological body, the equivalent thermal conductivity is usually considered a fixed value or is predicted according to the embedded empirical formula [17, 21, 22]. In theory, this approach is not accurate. On one hand, the equivalent thermal conductivity is affected by many physical soil factors, and porosity and saturation are considered the most important factors [20, 23, 24]. The riparian zone consists of a saturated zone and an unsaturated zone, and the equivalent thermal conductivity varies with the degree of saturation, which is different from that of the riverbed hyporheic zone. On the other hand, although some empirical formulas of the equivalent thermal conductivity have been embedded in the relevant simulation software, with the continuous development of thermal conductivity prediction models, more reasonable and accurate models may be ignored. Therefore, we need to find a suitable modeling approach, so that different equivalent thermal conductivity models can be easily considered in hydrothermal coupling methods.

In the history of soil thermal conductivity model development, many scholars have made many contributions. Some of researches obtained formulas by analyzing and fitting test data (i.e., empirical models), and some of them established models through theoretical analysis to predict the thermal conductivity of the soil (i.e., theoretical models) [25, 26]. Empirical models are generally obtained by establishing the relationship between the equivalent thermal conductivity and saturation based on normalizing the difference between the thermal conductivity of the soil in the saturated state and dry state [27]. However, theoretical models are typically adapted from predictive models of other physical properties (e.g., hydraulic conductivity, dielectric permittivity, and electronic conductivity) obtained from mathematical algorithms when the thermal conductivity of each component and the corresponding volume fractions are given [28]. Compared with empirical models, theoretical models have some shortcomings, such as low prediction accuracy, complex formulas, and difficult-to-determine parameters [20]. They are difficult to be directly applied to the actual soils. Therefore, this section summarizes several representative equivalent thermal conductivity empirical models of soils (Table 1).

On the basis of our previous research works [18, 33], the major objectives of this study are to (1) provide a new modeling approach for hydrothermal coupling models to improve the estimation accuracy of lateral temperature changes in the riparian zone; (2) compare and analyze the simulation effects of different equivalent thermal conductivity models based on numerical simulation and laboratory temperature records; and (3) investigate the sensitivity of the model parameters to reduce the workload of model correction.

#### 2. Mathematical Model of a Hydrothermal Coupling Model

##### 2.1. Seepage Field

The Richards equation based on conservation of mass and Darcy’s law is used for the numerical simulation of water movement in saturated-unsaturated zones [34]:where *p* is the pressure, Pa; is the acceleration of gravity, m/s^{2}; is the density of water, kg/m^{3}, which is assumed to be constant; *μ* (*T*) is the dynamic viscosity of water, (Pa·s), which is a function of temperature [35]; *z* is the elevation, *m*; *C*_{m} is the negative reciprocal of the slope of the soil water characteristic curve (SWCC), m^{−1}; *S*_{e} is the relative saturation of the unsaturated zone; *k*_{s} and *k*_{r} (*θ*) are the hydraulic conductivities of the saturated and unsaturated states, respectively, m/s; *S*_{s} is the elastic water storage rate, Pa; ∇ is the Laplace operator; and *Q*_{m} is the source term of the seepage field, kg/(m^{3}·s).

The SWCC of the unsaturated zone is described by the VG model [36]:where *α*, *β*, and *m* are the VG model parameters, of which *α* is the reciprocal of the intake value of the SWCC, m^{−1}; *β* is the indicator parameter of the slope of the SWCC; it reflects the distribution of soil pores and is obtained by fitting the SWCC; *m* = 1 − 1/*β*; *θ*_{r} is the residual water content, m^{3}/m^{3}; *θ*_{s} is the saturated water content, m^{3}/m^{3}; and *h*_{p} is the pressure head (*h*_{p} = *p*/, where *p* is the pressure, Pa, and is the volume weight of water, N/m^{3}), *m*.

Equations (2)–(5) reflect the relationship between relative saturation (*S*_{e}), volumetric water content (*θ*), the pressure head (*h*_{p}), the water capacity (*C*_{m}), and the relative hydraulic conductivity of the unsaturated zone (*k*_{r} (*θ*)); these equations can be directly applied in mathematical models of the unsaturated zone.

##### 2.2. Temperature Field Model

Referring to the report of Healy and Ronan [17], the saturated-unsaturated heat transfer model is expressed by the following equation:where *θ* is the volumetric water content, m^{3}/m^{3}, which is equal to the porosity in the saturated zone; is the heat capacity of water, J/(m^{3}°C); *C*_{s} is the heat capacity of soil, J/(m^{3}°C); *T* is the groundwater temperature, °C; *K*_{eq} (*θ*) is the equivalent thermal conductivity, W/(m^{3}°C); *D*_{H} is the hydrodynamic dispersion coefficient, m^{2}/s; *u* is the average velocity of groundwater (, where is the Darcy seepage velocity), m/s; and *Q*_{s} is the source term of the temperature field, W/m^{3}.where *α*_{T} and *α*_{L} are transverse dispersion and longitudinal dispersion, respectively, *m*; *δ*_{ij} is a kriging constant equal to 1 when *i* = *j* and 0 otherwise; is the vector of velocity in direction *i*; is the vector of velocity in direction *j*, m/s; and is the magnitude of velocity vector, m/s.

#### 3. Implementation of the Hydrothermal Coupling Model in COMSOL

In the hydrothermal coupling problem involving the riparian zone, the interaction between the land and water fields includes two components: (1) water movement will cause a change in temperature in the riparian zone and (2) a lateral temperature change in the riparian zone will cause a change in the physical state of the water flow. These two processes are highly coupled. In this paper, COMSOL software (COMSOL Multiphysics 5.2a, COMSOL Inc., Stockholm, Sweden) was used to solve this strong coupling problem. Notably, COMSOL can be used to solve linear and nonlinear problems, as well as time-dependent steady-state and transient problems. The material properties and boundary conditions in the model can be defined as constants or functions. In addition, the software contains predefined modules (such as fluid flow, heat conduction, and structural mechanics) and a self-defined mathematical module, which makes it particularly suitable for solving complex physics-related coupling problems [37].

In this paper, the temperature field model is a highly nonlinear partial differential equation (PDE), and there is no fixed calculation module to determine the coupling effect between the temperature (*T*) and the equivalent thermal conductivity (*K*_{eq} (*θ*)). Therefore, the PDE equation with temperature (*T*) as the dependent variable needs to be obtained by entering the derived partial differential equation into the self-defined mathematical module. The standard form of the coefficient PDE in the self-defined mathematical module is given as follows:where *u* is the dependent variable; *e*_{a} is the mass coefficient, kg/(m·s°C); *d*_{a} is the damping coefficient, J/(m^{3}°C); *c* is the diffusion coefficient, W/(m°C); *α* is the conservation flux convection coefficient, W/(m^{2}°C); *γ* is the conservative flux source term, W/m^{2}; *β* is the convection coefficient, W/(m^{2}°C); a is the absorption coefficient, W/(m^{3}°C); and *f* is the source term, W/m^{3}. These coefficients are defined by the users according to the governing equations.

The temperature field model (equation (6)) should be rearranged according to the form of equation (8); then, the above coefficients are defined in the corresponding edit box. Therefore, in equation (8), da ∂*u*/∂*t* can be set as the time-varying term in equation (6) (i.e., the left side of equation (6), with terms and , can be considered a combination of the conduction term (i.e., the first term of the right side of equation (6), ∇·*K*_{eq} (*θ*) ∇*T*), dispersion term (i.e., the second term of the right side of equation (6), ∇·*θ**D*_{H}∇*T*), and convection term (i.e., the third term of the right side of equation (6), ∇·*θ**uT*) in equation (6). Additionally, *f* can be set as the temperature field source term (i.e., *Q*_{s}) in equation (6) (where *e*_{a} = 0, *d*_{a} = *θ* + (1 − *n*) *C*_{s}, *c* = *K*_{eq} (*θ*) + *θ**D*_{H}, *α* = 0, *γ* = 0, *β* = *θ**u*, *a* = 0, and *f* = 0). Thus, the temperature field model can be equivalent to the coefficient PDE, and the equivalent thermal conductivity *K*_{eq} (*θ*) can be easily adjusted in the edit box. The equivalent temperature field model is coupled with the equivalent thermal conductivity model by *K*_{eq} (*θ*), and the water content *θ* in *K*_{eq} (*θ*) is used to couple the temperature field model and the seepage field model (the seepage field model is implemented by the groundwater flow module built-in COMSOL).

#### 4. Laboratory Investigations and Model Settings

The experiment was performed in a 3-D sand tank with a length of 60 cm, a width of 20 cm, and a height of 80 cm. The upstream and downstream sections of the sand tank were separated into an upstream water inlet tank and a downstream water outlet tank with a width of 10 cm by a plexiglass plate. The upstream water inlet tank had overflow ports of 30 cm and 50 cm to maintain a stable infiltration head, and an outlet was arranged at a height of 5 cm of the downstream water outlet tank. A schematic diagram of the entire test device structure is shown in Figure 1(a). The soil material used in the test is medium and fine sand with a bulk density of 1560 kg/m^{3}, and the particle size distribution of the studied sands is shown in Figure 2. In addition, 30 PT 100 temperature sensors (accuracy: ±0.15°C) were arranged on the front of the sand tank (Figure 1(b)) and connected to a computer through an A-D (analog-digital converter) box to obtain the real-time temperature data series at each monitoring point along the sand tank. Additional details on the laboratory test device and methods were provided by Ren et al. [38]. The infiltration water temperature was 9.5°C, the infiltration head was 25 cm, and there was no radiation condition (i.e., the influence of solar radiation on the soil surface temperature was not considered) to study the dynamic hydrothermal characteristics in the riparian zone under the condition of low-temperature water infiltration.

**(a)**

**(b)**

Based on the established hydrothermal coupling model and observed data from laboratory experiments, the finite element solution of the hydrothermal coupling model was obtained using a combination of the groundwater flow module and PDE module in COMSOL software. A mapping grid was used for the computational domain, and the grid size was 1.5 cm. The whole computational domain included 2160 elements and 2255 grid nodes. The computational domain and mesh layout are shown in Figure 3. In terms of boundary conditions, for the seepage field, the infiltration boundary on the left side (OD) was assumed to be a constant head boundary (i.e., infiltration head = 25 cm); the boundary on the right side (AB) was a permeable boundary, and a detailed explanation of this type of boundary can be found in the literature [39]. The other boundaries (OA, BC, and CD) for the seepage field were set as nonflux boundary. Because the temperature field model was equivalent by the coefficient PDE, two types of boundary conditions are applicable, namely, Dirichlet and Newman boundary conditions, which can be expressed as follows:where *h* is a custom coefficient; *r* and are the source terms at the boundary; and *n* is the outer normal direction of the boundary. The infiltration boundary on the left side (OD) was used a Dirichlet boundary (i.e., water temperature = 9.5°C, achieved from equation (9)). The other boundaries (OA, AB, BC, and CD) for the temperature field were Newman boundaries (i.e., adiabatic boundary conditions, achieved by equation (10)).

**(a)**

**(b)**

For the initial condition of the seepage field, it is assumed that the pressure head is 0 m and the initial condition of the temperature field is set to 20°C according to the initial temperature of sand in the experiment.

According to the laboratory experiment and related literature, the calculation parameters of the hydrothermal coupling model are given in Table 2. Among them, the values of parameters *k*_{s}, *θ*_{s}, *θ*_{r}, and *ρ*_{b} were obtained from the laboratory experiment, and the values of parameters *C*_{s} and were taken from the research of Naranjo and Simth involving sandy soil [19]. The values of the parameters *α*, *β*, *C*_{clay}, *C*_{sand}, *C*_{silt}, and *C*_{om} were obtained from Carsel and Parrish [40], and *n* and *q* were taken from of He et al. [41]. The soil density, saturated water content, residual water content, and soil properties of the abovementioned sandy soils were obtained from the previous studies [19, 40, 41] of similar sandy soil. Therefore, other parameters (*C*_{s}, , *α*, *β*, *C*_{clay}, *C*_{sand}, *C*_{silt}, and *C*_{om}) are also considered similar and can provide a reference for the findings of this study. The thermal conductivity of water is set to 0.58 W/(m°C), and the density of water is 1000 kg/m^{3}.

#### 5. Results and Discussion

##### 5.1. Model Validation

To quantitatively evaluate the simulation effect of the hydrothermal coupling model, the root-mean-squared error (RMSE), coefficient of determination (*R*^{2}), and relative error (Re) are used as the model evaluation indices [42]. Figure 4 shows a comparison of the observed and simulated temperature values at different monitoring points based on the equivalent thermal conductivity model of Ren et al. [33], and the results of the evaluation indices at different monitoring points are given in Table 3. During the test, the data at monitoring points T5, T6, T12, T14, T17, T22, T29, and T30 failed to be obtained because of sensor failure. Figure 4 and Table 3 show that the RMSE values ranged from 0.10–2.02°C with an average value of 0.68°C, and there is a small deviation between the simulated and observed values. In addition, we found that the large RMSE values for monitoring points T2, T4, T7, T9, and T10, which indicates that there is a large deviation between the simulated and observed values at these monitoring points. The value of *R*^{2} varied from −1.17 to 1.00, and small *R*^{2} values were observed at monitoring points T20, T21, T23, T24, T25, T26, T27, and T28. These monitoring points were all located near the upper surface of the sand tank, which implies that the consistency between the simulated temperature and the actual temperature is relatively poor when the established model is used to estimate the temperature of the shallow riparian zone. However, the RMSE values at these monitoring points were relatively small, which indicates that the deviations between the observed and simulated temperatures are small. The Re value ranged from 0.46–10.89% with an average value of 3.58%, and large values were only occurred at monitoring points T4 and T10. Therefore, the established hydrothermal coupling model can reasonably reflect the lateral temperature changes in the riparian zone during the water infiltration.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

##### 5.2. Evaluation of the Simulation Results of Different Equivalent Thermal Conductivity Models

The simulation effects of different models vary with the monitoring points. Therefore, it is necessary to comprehensively evaluate the effect of each equivalent thermal conductivity model. In this section, the simulated and observed temperature data collected at all the monitoring points in the sand tank are considered for global comparison, and the results are shown in Figure 5. Notably, the observed and simulated temperature values of the Ren model, Lu model, and Côté and Konrad model are mainly concentrated near the 1 : 1 line, which indicates that the lateral temperature changes simulated by the hydrothermal coupling model based on these equivalent thermal conductivity models are close to the observed temperature variation overall, and the deviation is small. However, the observed and simulated temperature values of the Campbell model were scattered near the 1 : 1 line, and these scattered data are mainly in the high-temperature period, which indicates that the temperature curve simulated by the hydrothermal coupling model under the Campbell model deviates from the actual temperature curve in the initial period and yields values higher than the actual temperature values.

To quantitatively evaluate the simulation effects of the hydrothermal coupling model with different equivalent thermal conductivity models from a global perspective, the statistical results of the simulated and observed temperature values at 22 monitoring points in the sand tank for different thermal conductivity models are given in Table 4. According to the results in Table 4, the average RMSE, *R*^{2}, and Re values for the Ren model, Lu model, and Côté and Konrad model were 0.86°C, 0.95, and 3.57%, respectively, and each evaluation index reflected a satisfactory result. The Lu model is a thermal conductivity model that has been widely used in recent years [43, 44] and provides a good simulation effect [18]. Research by Ren et al. [18] revealed that coupled hydrothermal modeling based on the Lu model yielded better results than the traditional modeling technique, suggesting that changes in the equivalent thermal conductivity model can promote the accuracy of coupled riparian zone modeling. The results of this study not only show that the Lu model has some advantages in the estimation of lateral temperature changes in the riparian zone but also further indicate that the thermal conductivity model proposed in our previous work [33] yields excellent performance.

Figure 6 compares the simulated and observed mean temperature variations at 22 monitoring points in the sand tank for different equivalent thermal conductivity models. Notably, during the infiltration of low-temperature water, the temperature of the soil inside the sand tank is decreasing; however, the temperature curve simulated based on the Campbell model plots among the uppermost temperature curves simulated by the other equivalent thermal conductivity models, which indicates that the values of the equivalent thermal conductivity are underestimated by the Campbell model. Overall, from the mean temperature curves simulated based on different equivalent thermal conductivity models in Figure 6, the values of equivalent thermal conductivity predicted by the models display the following order from large to small: Ren model > Lu model > Côté and Konrad model > Johansen model > Campbell model.

From the above analysis results, the simulation effects of the hydrothermal coupling model vary with the model of equivalent thermal conductivity used. Therefore, the results suggest that it is necessary to select a reasonable equivalent thermal conductivity model when using a hydrothermal coupling model to analyze the hydrothermal dynamics in the riparian zone. The traditional methods of coupled hydrothermal modeling are generally based on the default module of heat transfer in porous media and focus on the internal heat transfer of objects, which limits the selection of an equivalent thermal conductivity prediction model. However, the PDE module is used in this paper to model the heat transfer equation, which makes the selection of the equivalent thermal conductivity model more flexible. These results can provide a reference for the promotion and reasonable selection of an equivalent thermal conductivity model to simulate temperature variations in saturated-unsaturated sandy soils.

##### 5.3. Sensitivity Analysis

The established hydrothermal coupling model contains many parameters and requires extensive model calibration. Therefore, the Morris method was selected to analyze the global sensitivity of the parameters in the hydrothermal coupling model with the Ren model to the temperature output. The specific calculation principle of the Morris method can be found in [45].

The model parameters selected in this paper include: hydraulic conductivity (*k*_{s}), saturated water content (*θ*_{s}), residual water content (*θ*_{r}), porosity (*n*), VG model parameters (*α*, *β*), and heat capacity of soil (*C*_{s}). In the numerical calculation process, each parameter was evaluated based on a 10% increase, and the calculation time was 500 minutes. The difference between the temperature field and the rated temperature field due to changes in the parameters was determined for each finite element node, and the average change was used to represent the change in temperature. According to the calculation principle of the Morris method, the combination of parameters was input into the hydrothermal coupling model, and then, the sensitivity of the parameters and the corresponding interactions was obtained by solving the model in turn. Table 5 shows the calculation results for 28 groups of parameter combinations and the global sensitivity obtained by the Morris method.

To facilitate the analysis, the temperature fluctuation results for each parameter combination in Table 5 are presented in a histogram (Figure 7). Table 5 and Figure 7 show that when only a single parameter changes (Group 1–7), *k*_{s}, *β*, *C*_{s}, and *α* have obvious effects on the temperature output of the model, followed by *θ*_{s}, *n*, and *θ*_{r}. When multiple parameters change (Group 8–28), groups 8, 11, 14, 16, 18, 19, 22, and 25 are the parameter combinations that cause the largest temperature fluctuations, and all of these parameter combinations include *k*_{s}, *β*, *C*_{s}, or *α*, which suggests that *k*_{s}, *β*, *C*_{s}, and *α* are still the main factors that influence the temperature output when multiple parameters are considered. As shown for group 28, although all parameters increased by 10%, the temperature fluctuation was 0, which reflects a mutual influence among the parameters in the model, which was not simple numerical superposition.

In addition, Table 5 and Figure 7 show the temperature fluctuations caused by different parameter combinations with both positive and negative effects. The positive effects indicate that the increases and decreases in the parameters are consistent with the changes in temperature in the process of parameter correction, and the negative effects reflect the opposite trends. Therefore, these positive and negative effects can provide guidance for model parameter correction. In this study, if the simulated temperature values were less than the experimentally observed data, we increased the positive-effect parameters or decreased the negative-effect parameters so that the simulation process approached the actual situation.

Although the sensitivity of hydrothermal coupling model parameters to the calculation output results has been assessed by some scholars in previous studies [11, 18, 19], different modeling approaches and calculation methods may have an impact on the results. Naranjo and Smith [19] used VS2DH software to assess the sensitivity of parameters in a hydrothermal coupling model they established, and the results showed that the hydraulic conductivity (*k*_{s}) and VG model parameter (*α*) were the sensitive parameters that affect the model output; notably, *k*_{s} was the most important and sensitive parameter. However, their results were based on a single-factor sensitivity analysis and did not consider the influence of multiple factors on the results. Ren et al. [11] used the HYDRUS-2D model to analyze the dynamics of water flow and heat transport in riparian zones, and a sensitivity analysis of single factors and multiple factors was performed. The results of their studies not only revealed the internal factors (hydraulic conductivity (*k*_{s}) and VG model parameter (*α*)) that are sensitive to the temperature output of the model but also showed that the water head and temperature (external factors) are sensitive to the results. On the basis of a previous study [11], Ren et al. [18] further used the orthogonal test method to perform a single-factor sensitivity analysis of model parameters, but the results were different from the previous findings. The sensitivity analysis results of Ren et al. [18] showed that in addition to the hydraulic conductivity (*k*_{s}) and VG model parameter (*α*), the heat capacity of soil (*C*_{s}) and VG model parameter *β* also considerably influence the model results. This difference may be caused by the difference between the mathematical model and the equivalent thermal conductivity model. Therefore, different modeling approaches or calculation methods will indeed have an impact on the sensitivity analysis results. However, combined with the results of this paper, in general, the hydraulic conductivity (*k*_{s}) and VG model parameter (*α*) are still the most important and sensitive parameters, and these two parameters mainly affect water movement in the soil. Other parameters, such as porosity (*n*) and heat capacity of soil (*C*_{s}) influence the heat transport process.

#### 6. Conclusions

Based on the seepage theory of unsaturated soil and the principle of energy conservation, a coupled solution to the PDE of the temperature and seepage fields is achieved through COMSOL software and the appropriate equivalent thermal conductivity model. Combined with temperature series data from laboratory experiments, the validity of the model was verified, and the results showed that the established model can effectively reflect the dynamic temperature changes in the soil based on the results of a sand box test.

A reasonable selection of the equivalent thermal conductivity model will improve the effect of hydrothermal coupling simulations. The temperature data observed in the laboratory test were compared with the simulated temperatures based on each equivalent thermal conductivity model, and the Ren model, Lu model, and Côté and Konrad model are recommended for future hydrothermal coupling analysis in saturated-unsaturated soils. However, the experimental material used in this paper was sandy soil; therefore, the results are only applicable to the sandy soil. The suitability for other soil types must be further confirmed by applying the proposed modeling approach in subsequent studies.

The results of sensitivity analysis revealed that the hydraulic conductivity (*k*_{s}), VG model parameters (*α* and *β*), and heat capacity of soil (*C*_{s}) greatly affect the temperature output of the hydrothermal coupling model. Therefore, in the application of the model, the observation accuracy of these parameters should be improved; thus, the numerical model that yields simulation results similar to the measured values can be obtained, even if the secondary factors are not fully considered or ignored.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Authors’ Contributions

W. Z. and J. R. designed the framework and wrote the manuscript. W. Z. and G. C. collected the data. W. Z. and Z. S. were mainly responsible for the simulation analysis. W. Z. and L. X. verified the results of the model. J. R. and Z. S. provided funding support.

#### Acknowledgments

This research was funded by the National Natural Science Foundation of China (Grant nos. U1765205 and 51679194), the Natural Science Foundation of Shaanxi Province (Grant no. 2020JM-448), the Planning Project of Science and Technology of Water Resources of Shaanxi (Grant no. 2019slkj-12), the Belt and Road Special Foundation of the State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering (Grant no. 2019490711), and the Priority Academic Program Development of Jiangsu Higher Education Institutions (YS11001).