Validation of an analytical model of groundwater velocity based on laboratory test and numerical simulation

Surface water–groundwater exchange affects the material and energy transfer of rivers and adjacent riparian zones. As an intuitive carrier of energy, temperature can effectively reflect the spatial and temporal variation of surface water–groundwater exchange process. In this paper, the influence of water head variation and sand sample uniformity on its temperature field and seepage field is studied through a one-dimensional sand column laboratory test. To verify the accuracy of the one-dimensional vertical heat analysis model, the vertical seepage velocity measured in the indoor test is compared with the vertical submerged exchange rate calculated by the four analysis models. The results show that the Hatch analytical solution, Keery analytical solution, McCallum analytical solution and Luce analytical solution calculated by VFLUX2 through MATLAB are reliable for calculating the vertical undercurrent exchange rate of the heterogeneous sand column.


Introduction
As one of the most crucial water resources on earth, rivers have always been the focus of research in the fields of hydraulic engineering, environmental engineering, and earth science, especially the relationship between surface water and groundwater (Dun et al. 2014;Kanduč et al. 2014;Lu et al. 2018) are gaining importance in water resource management policies. Continuous, dynamic interactions between these two water bodies sustain the biogeochemical and ecological processes at the water-groundwater interface in rivers, lakes, and wetlands (Hayashi and Rosenberry 2002;Vogt et al. 2012), and these dynamic interactions occur mainly in the subsurface of the riverbed and riparian zone, where water flows across the riverbed and banks into shallow subsurface flow paths and back into the river (Trauth et al. 2013;Kiel and Cardenas 2014;Lane et al. 2018;Chow et al. 2019;Lu et al. 2020), this process is always accompanied by the transport and exchange of water, heat and solute, which is very important for flood storage, significant impact on water quality and ecosystem health (Qureshi and Harrison 2001;Wilding et al. 2014;Swanson et al. 2017).
Irrigation water in rivers, streams, lakes, wetlands, irrigation water in bare canals, offshore seawater and other porous media passing through the bottom of water bodies will exchange with groundwater in most cases. The exchange process between surface water and groundwater is always accompanied by the transfer and exchange of energy. Among them, temperature, as an intuitive carrier of energy, is a characterization factor that can reflect the temporal and spatial changes in the undercurrent exchange process. As a natural tracer, it is easy to observe and noncontaminating. (Engelhardt et al. 2011(Engelhardt et al. , 2013. In recent years, In recent years, with the development of integrated sensor and computer coding, it is possible to record a large number of temperature time series records at multiple locations and multiple depths (Swanson and Cardenas 2010;Gordon et al. 2012). Applying a 1D heat transport model to a high-resolution temperature time-series data set collected using distributed temperature sensing (DTS) enabled the mapping of the spatial variation of vertical flux with depth (Vogt et al. 2010). Automated computational methods have facilitated the more flexible recording of high-resolution DTS temperature data to characterize fine-scale variations in depth-dependent flux at multiple sites (Fritz and Arntzen 2007;Gordon et al. 2012). Exchange rates are derived from temperature time series generated by groundwater and heat transport (Cardenas 2010;Lautz 2010). Vogt et al. (2012) used distributed optical fiber temperature sensors (DTS) to continuously observe the temperature of shallow groundwater in the riparian zone. By analyzing the temperature time series of the vertical section, they found that the groundwater temperature in the riparian zone is not spatially related. Uniform distribution, but non-uniform distribution with groundwater velocity at different depths. Chen (2020) uses the real-time monitoring data of water temperature and water level to study the distribution characteristics of the temperature field of the subsurface layer in the riparian zone in different seasons and different spatial locations, and uses the water temperature data to calculate the groundwater flow rate by the Hatch amplitude method (Hatch et al. 2006) and the COMOSL Multiphysics numerical solution method comparison of the computed results: although the two are different in magnitude, the overall variation law of the computed results is the same. Zhu et al. (2013) used the river bed temperature as a tracer variable, based on the analytical solution of the one-dimensional steady-state vertical heat transport problem, and revealed by continuous monitoring of the temperature of different vertical profiles of the river bed in the field. The distribution law of the depth of the undercurrent zone is revealed. Ji et al. (2018) took the continental beach in the Manwan Reservoir area of the Lancang River as the research object. They monitored the changes in water level and water temperature in the monitoring well of the continental beach in real-time during the operation of the reservoir. Based on the water level and water temperature data, the flow rate was accounted for by the hydrodynamic and temperature tracing methods, respectively. Scholars such as Irvine et al. (2015) used the sand column test, adopted the one-dimensional heat transport theory, and studied the exchange volume between surface water and groundwater through the temperature amplitude ratio and phase difference; Halloran et al. (2016) analyzed the measured temperature data by arranging temperature sensors and using the one-dimensional heat transport theory in field experiments, and concluded that the method of calculating seepage velocity using the amplitude method is the most reliable method for calculating seepage velocity in the subsurface zone; Conant (2004) identified the existence of the sub current area and delineated the scope of the subsurface zone by the temperature tracing method; Alexander and Caissie (2003) and Briggs et al. (2012) adopted the temperature tracing method to calculate, the water exchange between surface water and groundwater in the subsurface zone was calculated through the temperature data of continuous dynamic monitoring, and the relevant hydrogeological parameters in the area were estimated. Wondzell et al. (1996) quantitatively studied the undercurrent convection flow and groundwater inflow in river gravel ridges by monitoring data from field experiments and establishing groundwater models. The basic theory of the exchange of undercurrents in the riverbed was established.
Surface water-groundwater exchange rates vary significantly in time and space, making them difficult to quantify, especially at the multi-spatial and long-term resolution required for many field applications (Kalbus et al. 2006;Fritz and Arntzen 2007;Fanelli and Lautz 2008;Briggs et al. 2012). With the accumulation of riverbed sediments over time, the surface water-groundwater exchange rate is affected by solar radiation and air temperature, and the temperature field of river surface water and undercurrent zones are characterized by diurnal, monthly and seasonal fluctuations (Kasahara and Hill 2006;Rosenberry and Pitlick 2009). Phenomena such as melting snow, blizzard, tide, and reservoir operation are the main factors that cause the water level fluctuation in the river reach. Fluctuating of the river water level is an important factor in the confluence of river water and groundwater in the riparian zone. Along with the change of the upstream river flow, the river water level fluctuates. The fluctuation of the water level increases the wetting area of the riparian zone and thus expands the range of river water and riparian groundwater exchange (Taniguchi 2003;Harvey et al. 2013); such fluctuations contribute to the changing hydraulic gradient between river water and groundwater, increasing the rate of surface water-groundwater exchange (Edwardson et al. 2003). In the past few years, there have been some notable advances in the application of 1D heat transport models for surface water-groundwater exchange, and there are many practical examples of applying these results (Fanelli and Lautz 2008;Anderson et al. 2010;Lautz et al. 2010). Hatch et al. (2006) and Keery et al. (2007) applied a one-dimensional heat transfer model to infer surface water-groundwater exchange rates from records of subsurface temperature over time.
Through the research of many scholars above, we find that the existing research is to obtain a relatively accurate one-dimensional vertical heat analytical model by comparing the analytical solution of the one-dimensional vertical heat transport problem with the numerical solution obtained by the hydraulic method. Still, there is a lack of corresponding test verification. Therefore, in this paper, a one-dimensional homogeneous sand column indoor test is carried out based on a one-dimensional thermal transport analytical model. The head changes are simulated by the indoor tests. The effects of water level rise and fall processes on the temperature and seepage fields of homogeneous and inhomogeneous sand columns are investigated. Comparing the temperature data measured in the laboratory test to calculate the groundwater flow rate and the measured seepage water. The vertical seepage velocities measured in the laboratory test were compared with the vertical underflow exchange rate calculated by the four analytical models. To verify the accuracy of the one-dimensional vertical heat analysis model.

Methods
Temporal and spatial variations affect the rate of surface water-groundwater exchange. Hatch et al. (2006) and Keery et al. (2007) apply a one-dimensional heat transfer model Surface water-groundwater exchange rates are derived from recordings of subsurface temperature over time. The existing research have to obtained a more accurate analysis model of one-dimensional vertical heat transport by comparing the analytical solution of the one-dimensional vertical heat transport problem with the numerical solution obtained by the hydraulic method. However, there is a lack of corresponding experimental validation. Therefore, in this paper simulates the head change by comparing indoor tests. The groundwater flow rate calculated by the measured temperature data is compared with the groundwater flow rate obtained by the measured water infiltration to verify the accuracy of the one-dimensional vertical thermal analysis model.

One-dimensional thermal transport equation
The exchange modes of heat transport in surface water-groundwater exchange (shown in Fig. 1) mainly consist of two modes: heat conduction (heat transfer through solids in sediments and water flow) and thermal convection (heat transfer through water flow). Stallman (1963) assumed that the geological body is a saturated homogeneous isotropic porous medium; In porous media, solid-liquid two-phase transport is isotropic, where water flow is stable, and heat transport is variable. The solid phase skeleton and pore fluid are in a state of local thermal equilibrium, and their thermal characteristics remain unchanged in time and space. A fundamental equation describing the transport of water and heat in saturated porous media was proposed (Shanafield et al. 2011). Therefore, the mathematical model of heat transport in the subsurface flow zone is as follows: In the formula, T is the temperature (°C); t is time (s); ρ w and C w are the density (kg/m 3 ) and mass heat capacity (J/kg·°C) of water, respectively; ρ and C are the equivalent density (kg/m 3 ) and equivalent mass heat capacity (J/kg °C) of the saturated porous medium, respectively; u x , u y and u z are the seepage velocity (m/s) in the x, y and z directions; λ is the equivalent thermal conductivity of the saturated porous medium [w/(m s)].
In addition, due to the complex structure of the undercurrent zone, the surface water-groundwater exchange has a high degree of variability in time and space (Vogt et al. 2010;Liu Dongsheng et al. 2017). In addition to the overall material properties of saturated porous media, the exchange rate also needs to comprehensively consider the effects of seasons, rainfall, solar radiation, and snowfall (Song et al. 2019).
In 1965, Stallman (1965) proposed a one-dimensional heat transfer equation based on the assumption of onedimensional vertically non-uniform steady migration of flow in porous media (Lautz 2012).
In the formula, K = ρCπ/Pk; V = ρ w C w u z /2λ; a and b are constants.

Existing models
Stallman's pioneering work laid the foundation for subsequent developments in the field of surface and groundwater exchange. Due to the improvement of computer data processing performance, innovation of temperature measurement instruments and acquisition methods, the many scholars later improved the Stallman analysis model based on this. The improved model is more suitable for complex boundary conditions and hydrogeological environments (Lautz 2012).

Hatch solution
In the formula, υ t is the temperature front-end movement speed (m/s): The right side of Eq. (5) represents the changes in the amplitude and phase of temperature fluctuations with the increase of sediment depth, respectively. Hatch et al. (2006) divided Eq. (5) into two components and solved the amplitude and phase of the temperature between the measurement points at two different depths, respectively thus obtaining: In the formula, u zAr and u zΔφ are the vertical groundwater seepage velocity (m/s) calculated from the amplitude and phase, respectively; Δz is the vertical distance of the temperature monitoring point (m); A r is the temperature amplitude A d of the deeper measurement point The ratio of the temperature measurement point amplitude A s at the shallow measurement point, As that is, A r = A d /A s (as shown in Fig. 2b); Δφ is the phase difference between the two monitoring points, which is the time lag (s). In Fig. 2, a represents a schematic diagram of the layout of common temperature monitoring instruments, and b represents the temperature time series temperature measurement curve of the deep and shallow measurement points. Keery et al. (2007) ignored the effect of thermal dispersion on heat transport in porous media, improved the Stallman model, and gave the analytical solution of one-dimensional heat transport equation: Fig. 2 Sketch for the temperature time-series curves of deep and shallow measurement points and the commonly layout approaches of temperature monitoring instrument (adapted from Schmidt et al. 2006) in the formula: H = ρ w C w /λ。

McCallum solution
Although the Hatch analytical method and the Keery analytical method can calculate the vertical seepage velocity relatively accurately, there is a certain difference between the vertical seepage velocity calculated by the amplitude ratio method and the phase lag method. The calculation results are affected by the thermophysical parameters of the porous medium, and the thermophysical parameters have certain uncertainties, which have a significant impact on the accuracy of the analytical model (Shanafield et al. 2011;. In response to the above problems, McCallum et al. (2012) optimized the equation: k e = Δz 2 P 2 ln A r 4 2 Δ 2 − P 2 ln 2 A Δ P 2 ln 2 A r + 4 2 Δ 2 P 2 ln 2 A r − 4 2 Δ 2 where u zArφ is the vertical seepage velocity (m/s) calculated by the amplitude-phase combination method.

Luce solution
Luce et al. (2013) The values of parameters a and b in the undetermined coefficients of the Stallman analytical solution can be approximately expressed as the amplitude ratio and phase lag of the day-night temperature signal divided by the vertical distance between the measuring points, namely: And from this, it follows: in the formula: η = -lnA r /Δφ; ω = 2π/P, in here k e = λ/ρC.

Test apparatus
The test device (shown in Fig. 3) consists of a plexiglass water tube, a high-precision temperature sensor (Omega), a water pump, a data collector, a computer, a measuring cup, an electronic scale, a rubber tube, a heating rod, a water bucket and a sand column. The plexiglass water pipe is divided into upper and lower parts, the inner diameter is 10 cm, the outer diameter is 12 cm, the upper length is 100 cm, and the lower length is 20 cm, which is connected by column-shaped connecting sections. The upper plexiglass water pipe is arranged every 10 cm from 5 cm above the bottom. The temperature sensor is arranged with a filter gauze at the bottom, the fine sand prepared in advance is placed, and the temperature sensors are compacted and placed at an interval of 10 cm. The height of the sand column is 70 cm. There are three water outlets and one water inlet above the sand column. The water in the bucket is heated by the heating rod, and the water flows from the bucket into the water pipe by the water pump, and flows into the bucket through the water outlet at different heights to form stable circulating water. A water outlet is set at the lower part, and the amount of water seeping from the sand column is measured. The purpose of the rubber hose is to connect the water pump, bucket, and plexiglass hose and to receive the water seeping out from the sand column.

Test method
Using sand column experiment, a one-dimensional flow of water and heat is generated in the porous medium. The upper end of the sand column is provided with a water flow with quasi-sinusoidal temperature change by a heating rod. We obtained through pre-experiments that it takes 15 min to heat 1 ℃ when the bucket is full of water. Continuous heating for 2 h can make the water temperature of the bucket rise to 30 ℃. After that, ice is added to lower the bucket water temperature by 1 °C every 15 min. In this way, the upper end of the sand column is provided with a sinusoidal variation of water flow from 22 to 30 ℃. Before the test began, the sand column was saturated with water flow through the sand column for several days until the flow measured by the measuring cup stabilized and the sand column had a uniform temperature distribution. The experiments were carried out for two different water head changes (rising and falling), two different soil types (homogeneous and heterogeneous), three different water heads (10 cm, 20 cm and 35 cm) and sinusoidal temperature changes. The test plan is shown in Table 1. Plan 1: The water pump provided tap water with an initial temperature of 22 °C to the plexiglass water pipe through the water inlet. The rubber pipe with the 35 cm water outlet was placed in the bucket to form a circulating water with a water head of 35 cm. When the water output from the bottom water outlet is stable, the temperature of the bucket is heated by the heating rod to provides a sinusoidal variation of temperature. The temperature data are read by the data collector. Then, place the 20 cm outlet rubber hose in the bucket to form circulating water with a head of 20 cm. After the water temperature is heated and cooled the rubber tube of 10 cm outlet is placed in the bucket to form circulating water with a head of 35 ~ 20 ~ 10 cm.
Plan 2: The water pump provided tap water with an initial temperature of 22 °C to the plexiglass water pipe through the water inlet. The rubber pipe with the 10 cm water outlet was placed in the bucket to form a circulating water with a head of 10 cm. Then, place the 20 cm outlet rubber hose in the bucket to form circulating water with a head of 20 cm. After the water temperature is heated and cooled the rubber tube of 35 cm outlet is placed in the bucket to form circulating water with a head of 10 ~ 20 ~ 35 cm.
Plan 3: The upper 10 cm of fine sand is dug out and filled with coarse sand of equal height. The water pump provides tap water with an initial temperature of 22 °C to the plexiglass water pipe through the water inlet, and the rubber pipe of the 10 cm water outlet is placed in the bucket. Circulating water with a head of 10 cm was formed. When the water output from the bottom outlet was stable, the temperature of the bucket was controlled by a heating rod. Then, place the 20 cm outlet rubber hose in the bucket to form circulating water with a head of 20 cm. After the water temperature is heated and cooled the rubber tube of 10 cm outlet is placed in the bucket to form circulating water with a head of 35 ~ 20 ~ 10 cm.

Results and discussion
For the one-dimensional thermal transport analytical model, the commonly used solvers are VFLUX2  and Ex-Stream (Swanson and Cardenas 2011). In this experiment, VFLUX2 was used to solve the vertical seepage velocity. VFLUX2 was developed through MATLAB computing language programming. The Hatch analytical solution, Keery analytical solution, McCallum analytical solution and Luce analytical solution can be used to calculate and analyze the vertical seepage velocity of the sand column. The rate of submerged flow exchange in the riparian zone derived from the above four analytical models is obtained. Table 2 lists the parameters used by the four analytical models in the VFLUX2 calculation program. The parameters are mainly based on Lautz et al. (2012) and laboratory measurements. Figure 4B shows that the dotted line represents the temperature amplitude of each monitoring point, and the solid line represents the temperature time series data after filtering. In the real environment, various noise signals can have an impact on the measured temperature data. The VFLUX2 program can filter natural and human-generated noise with an finite impulse response (Lautz 2012).

Influence of dropping water head of the homogeneous sand column
In this study, the sand column experiment generates a one-dimensional vertical flow in the sand column. In these one-dimensional flow fields, temperature changes are monitored by temperature sensors at different positions, and the flow velocity at different depths is calculated by VFLUX2. Since the phase method can only calculate the size, the direction cannot be determined, and only the seepage velocity calculated by the amplitude method is considered. The calculation results have the following characteristics: in the flow velocity calculation results of the amplitude method, the flow velocity at 0.1 m is also calculated, and the velocity range value calculated by the Hatch amplitude method is 3.66 × 10 -5 -8.85 × 10 -5 m/s, the velocity values calculated by the Keery amplitude method are 3.55 × 10 -5 -8.18 × 10 -5 m/s. For the amplitude-phase combination method, the flow velocity calculated by the McCallum solution and the Luce solution coincide, and the calculated flow velocity value is 7.30 × 10 -5 -8.90 × 10 -5 m/s, which is larger than that of the amplitude method. When the water head is 35 cm, the measured seepage velocity is about 6.9 × 10 -5 m/s; when the water head is 20 cm, the measured seepage velocity is about 5.9 × 10 -5 m/s; when the water head is 10 cm, the measured seepage velocity is 5.1 × 10 -5 m/s. As the water head decreases, the calculated results of each model gradually decrease. As we can see that the calculated results of the Hatch amplitude method at 0.1 m are in good agreement with the measured values, which is consistent with the previous research results of Lautz (2012). To compare the groundwater flow velocity at different depths, the seepage velocity calculated by the Hatch amplitude method at 0.10 m, 0.20 m and 0.30 m was selected as the research object, Fig. 4C. The value ratio is relatively small. In the depth direction, the seepage velocity decreases, and the flow velocity at 0.2-0.3 m is close to parallel, which is consistent with the field test results of Li et al. (2016) match.

Influence of rising water head on the heterogeneous sand column
According to the analysis, the four analytical models have high accuracy in both the calculated and measured values of the homogeneous sand column vertical flow velocity. The Hatch amplitude method has the highest agreement between the calculated results and the measured values. Figure 6 shows the curve of variation of vertical temperature and seepage velocity under the condition of rising head of nonhomogeneous sand column. Because the 0.05 m and 0.15 m temperature sensors are in different sand grains, the flow velocity at 0.2 m is selected for this working condition. The calculated value obtained by the Hatch amplitude method. The e velocity value is 9.23 × 10 -5 ~ 9.35 × 10 -5 m/s, and the velocity value calculated by the Keery amplitude method is 9.01 × 10 -5 ~ 9.11 × 10 -5 m/s, for the amplitude-phase combination method, the flow velocity calculated by the McCallum solution and the Luce solution coincide, and the calculated flow velocity value is 1.37 × 10 -4 ~ 1.38 × 10 -4 m/s, which is larger than the amplitude method. When the water head is 10 cm, the measured seepage velocity is about 1.18 × 10 -4 m/s; when the water head is 20 cm, the measured seepage velocity is about 1.33 × 10 -4 m/s; when the water head is 35 cm, the measured seepage velocity is 1.49 × 10 -4 m/s. As we can see from Fig. 6B that the flow velocity at 0.2 m calculated by the analytical model is smaller than the measured flow velocity, which is consistent with the results of working condition 1. The flow velocity calculated by McCallum solution and Luce solution is the largest, followed by the Hatch amplitude method, and Keery's approach was the lowest. On the whole, the magnitude of the flow rate calculated by the amplitude ratio method and the amplitude combination method is relatively consistent with the measured flow rate size and the overall change law of the calculated result. Analytical solution, McCallum analytical solution and Luce analytical solution are reliable in analyzing vertical seepage velocity. For the amplitude ratio method, the velocity results calculated by Hatch solution and Keery solution maintain a certain numerical difference. Keery et al. (2007) conducted a comparative analysis of the Hatch solution and the Keery solution, and concluded that the Hatch solution is more accurate than the Keery solution due to the consideration of thermal dispersion effects. From the above analysis, as we can see that it is more reliable to use the four models to calculate the vertical seepage velocity of the heterogeneous sand column. Therefore, in the subsequent calculation and analysis of the vertical flow velocity of the heterogeneous underflow zone, four analytical models can be used for analysis.

Conclusion
In this paper, the temperature data measured by the laboratory test and the measured seepage flow rate are verified by VFLUX2 to calculate and analyze the vertical seepage velocity through four MATLAB analytical models: Hatch analytical solution, Keery analytical solution, McCallum analytical solution and Luce analytical solution. The main conclusions are as follows: 1. Three kinds of indoor sand column tests were conducted.
The effect of decreasing head on homogeneous sand column; the effect of increasing head on homogeneous sand column and the effect of increasing head on nonhomogeneous sand column were studied, respectively. The vertical seepage velocity calculated by the measured seepage flow is compared with the solutions of the four analytical models, and the solution of the vertical flow velocity of the analytical model in the heterogeneous underflow zone is analyzed. 2. The vertical seepage velocity is calculated from the measured seepage volume.

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

Declarations
Conflict of interest The authors declare that are no conflict of interest. Fig. 6 Variation of temperature and seepage velocity of heterogeneous sand column with rising water head (A represents the temperature change curve of each monitoring point under different water heads; B represents the vertical seepage velocity calculated by various analytic models of VFLUX2)