A correction method for calculating groundwater flow by considering the influence of environmental factors on heat sources

Abstract

This invention discloses a correction method for estimating groundwater flow considering the influence of environmental factors on heat sources, relating to the fields of hydrogeological exploration and groundwater resource assessment. The method includes the following steps: Step 1, constructing a one-dimensional vertically homogeneous saturated porous medium assumption, obtaining basic parameters through field and laboratory measurements; Step 2, constructing a basic heat transport equation, using Stallman's one-dimensional vertical transient heat conduction-convection equation to describe the groundwater temperature field without rainfall disturbance; Step 3, rainfall disturbance correction, employing a dual correction method that introduces a rainfall heat source term and improves surface boundary conditions. This invention uses the above-mentioned correction method for estimating groundwater flow considering the influence of environmental factors on heat sources. Since rainfall heat is mainly concentrated at the surface for exchange, the boundary condition correction is simpler; considering the continuous vertical release of heat, the internal heat source correction is more optimal, effectively characterizing the disturbance of shallow ground temperature by rainfall events and improving the accuracy of groundwater vertical flux inversion.

Description

Technical Field

[0001] This invention relates to the field of hydrogeological exploration and groundwater resource assessment technology, and in particular to a correction method for estimating groundwater flow considering the influence of environmental factors on heat sources. Background Technology

[0002] Groundwater is a precious resource of the earth, a high-quality geological resource, a rich information carrier, and an active ecological environment factor. At the same time, it is also a disaster-causing factor that cannot be ignored. Therefore, knowing the actual flow of groundwater is of great value in water resource assessment, ecological environment protection, and engineering construction.

[0003] Currently, methods for calculating groundwater flow tend to simplify boundary conditions or fail to consider transient environmental changes such as common natural weather phenomena like rainfall. For example, the steady-state analytical method assumes a steady-state heat flow and directly solves for vertical flux using the curvature of the temperature-depth curve. It completely ignores surface temperature changes caused by factors such as climate warming / land use change, which can produce instantaneous fluctuations in the heat conduction profile, making the interpretation of vertical groundwater flux estimates more complex; and the transient process of instantaneous thermal shocks from precipitation and subsequent rebalancing. The transient analytical method introduces a linearly time-varying surface temperature boundary condition into the equations to simulate long-term climate warming trends. However, its boundary conditions are overly simplified: they can only handle linear or simple periodic changes and cannot characterize the "nonlinear, pulsed" temperature abrupt changes caused by precipitation events. Numerical inversion methods and other numerical models, based on numerical solutions of one-dimensional convection-conduction equations, invert flow by fitting observed temperature time series. While this method can handle complex boundary conditions, it is essentially a reverse deduction based on a black box of measured data, and does not explicitly reflect transient changes such as rainfall in the actual correction terms of the formula.

[0004] Current technology still has shortcomings in the following three aspects: Steady-state and simplified assumptions are divorced from dynamic reality: the actual heat flow field changes in real time under the combined influence of global warming, land use, precipitation, snow accumulation and other meteorological factors. The root of the error lies in the error of the theoretical premise.

[0005] Distortion of idealized physical property parameters: Even if transient changes are considered, they are only simple linear changes. Replacing the variable land-atmosphere interface temperature with only coarse assumptions will ignore the nonlinear changes and step recovery brought about by rainfall, and will lose key thermal disturbance information, resulting in bias in the model driving source itself.

[0006] (3) Signal inversion has equivalence: the inversion method will have the dilemma of multiple solutions. Completely different physical processes may produce extremely similar temperature patterns, making the inversion results lack uniqueness and the physical interpretation unclear.

[0007] The above methods are all based on the transient conduction-convection heat flow equation in the Cartesian z-direction (depth) (Stallman, 1965). However, most people treat this equation as an identity, ignoring the fact that sudden changes such as rainfall can cause new heat source terms, which leads to the aforementioned problems. Summary of the Invention

[0008] The purpose of this invention is to provide a correction method for calculating groundwater flow rate by considering the influence of environmental factors on heat sources. This correction method considers the influence of environmental factors such as rainfall on the temperature field and disturbances, and improves the accuracy of groundwater flow velocity inversion by introducing a heat source correction term.

[0009] To achieve the above objectives, the present invention provides a correction method for estimating groundwater flow based on the influence of environmental factors on heat sources, comprising the following steps: Step 1: Construct a one-dimensional vertically homogeneous saturated porous medium hypothesis, and obtain temperature variation curves at various depths over time through field monitoring and laboratory measurements. T ( z,t Hourly rainfall intensity I ( t ), thermal diffusivity saturated hydraulic conductivity Ks Changes in moisture content ; Step 2: Based on the assumption of a one-dimensional vertically homogeneous saturated porous medium, a basic heat transport equation is constructed. The one-dimensional vertical transient heat conduction-convection equation proposed by Stallman (1965) is used to describe the underground temperature field without rainfall disturbance. ; Step 3: Correct for rainfall disturbances by introducing a rainfall heat source term and improving surface boundary conditions.

[0010] Preferably, in step one, one-dimensional vertical thermal transport must be satisfied, the soil is regarded as a homogeneous isotropic medium, the thermal properties are uniformly distributed in space and do not change with depth; the soil is always in a saturated state within the study depth range, and the influence of the unsaturated zone is ignored.

[0011] Preferably, the step three, which involves introducing a rainfall heat source term correction, includes the following steps: S311. Introduce an explicit heat source term S on the right side of the equation in step two. i (z,t), assuming there are no other significant heat sources besides rainfall, we can deduce from the law of energy conservation: in, This indicates new sources and sinks of heat generated by rainfall. The density of water, The specific heat capacity of water, For the first i The rainfall intensity per unit time of a rainfall event. For the first i The temperature of the rainwater at a certain moment in this rainfall event. h ( z Let be the vertical distribution function of the heat source, using an exponential decay form: in, The characteristic attenuation depth represents the effective depth of the heat effect; S312. Integrating the equations, we obtain the modified heat transport equations. Substituting the heat source term from S311 into the basic equations, we obtain the one-dimensional heat transport governing equations considering rainfall corrections: Numerical solutions are obtained using the finite difference method or the finite element method, employing an explicit Euler scheme. The time step must satisfy the CFL stability condition to ensure computational convergence. ; Where Δz is the spatial grid spacing, which is generally taken as 0.01m; S313, Using measured temperature time series Inversion of unknown parameters q The global optimization algorithm is used to solve the problem, and the inversion objective is to minimize the root mean square error (RMSE) between the model-calculated temperature and the observed temperature. S314. Method validation and accuracy evaluation: Validation is performed using synthetic data and measured data.

[0012] Preferably, in S311 Controlled by peak intensity and If is a constant, then: in, Where I is the thermal diffusivity, and Ipeak is the maximum instantaneous intensity during a rainfall event. Take a value of 0.1 to 0.4.

[0013] Preferably, S313 includes the following steps: S3131. Construct a forward model. Using the modified heat transport equation as the forward model, numerical discretization is performed using the finite difference method. The computational domain of the model is depth z∈[0, Lz The spatial step size Δz is 0.01m; the time step size Δt must satisfy the stability condition of the explicit scheme, while the implicit scheme is unconditionally stable, thus ensuring computational accuracy and convergence. S3132. Define the objective function, using the measured temperature time series. To achieve the calibration objective, the root mean square error (RMSE) is chosen as the objective function to measure the simulated values. Degree of deviation from observed values: in, For sensor depth, N The total number of observations at all depths and all time points; S3133, Optimization algorithm selection and parameter range setting. v For the scalar to be inverted, the differential evolution algorithm is used, with the following parameters set: the population size is 15–20 times the number of parameters to be inverted, the maximum number of iterations is set to 100–200, and the mutation factor and crossover probability are set to default values. v Usually taken - m / s; S3134. Inversion Execution and Result Output: Forward Model Calculation of Current... v The value corresponds to the simulated temperature, and the objective function value is calculated. The output value minimizes the RMSE. v The value is used as the inversion result.

[0014] Preferably, in S314, the synthetic data uses the modified equation of the present invention to generate synthetic temperature data containing rainfall events, and adds measurement noise of ±0.02℃. The flow rate is inverted using the heat-free term and the method of the present invention, respectively, and the relative error is calculated. The measured data verifies the combination of real measurement conditions and real weather conditions by comparing the degree of fit between the original formula and the improved formula with the real values.

[0015] Preferably, step three, the improvement of surface boundary condition correction, includes the following steps: S321. Improved surface boundary conditions, considering the energy budget at the surface (z=0), including air convection heat exchange and precipitation heat flux, derived according to energy conservation: in, H The convective heat transfer coefficient, = T (0, t ) represents the surface temperature. The heat flow from the surface into the ground via conduction is equal to the sum of the convective heat exchange between the surface and the air and the net heat flux from rainfall. S322. Complete Mathematical Model and Numerical Solution: Combine the formula in S321 with the equation in step two, and add the boundary conditions and initial conditions to form a complete mathematical model. ; S323. Inverting vertical groundwater flux using measured temperature time series. Inversion of unknown parameters The inversion objective is to minimize the root mean square error (RMSE) between the model-calculated temperature and the observed temperature. ; S324. Based on actual observations of rainfall intensity, daily soil temperature variation, seasonal cooling trend, and station characteristics, simulated data are generated to compare the original formula and the improved formula with the observed data.

[0016] Preferably, S321 includes the following steps: The external heat flux at the Earth's surface is This is typically provided by air convection heat exchange: in, Tsurface = T (0, t ( ) represents the surface temperature. H The convective heat transfer coefficient is the lower boundary of the thin layer. z = ε The internal heat flow at that location is: Thin layer due to δ The total heat generated by the source term is: Ignoring the change in internal energy of the thin layer as ε approaches zero, the law of conservation of energy states: Let ε→0, then Substituting, we get: The boundary conditions for surface heat flux are obtained as follows: .

[0017] Preferably, in S322, the finite difference method or finite element method is used for numerical solution. At each time step, the internal nodes are first calculated based on the current temperature field, and then the surface temperature is updated through discretized boundary conditions. .

[0018] Therefore, this invention employs the aforementioned correction method for calculating groundwater flow considering the influence of environmental factors on heat sources. Based on the transient conduction-convective heat flow equation in the Cartesian z-direction (depth), and incorporating the Green-Ampt influence on soil infiltration capacity limitations, it is the first to add precipitation factors to the transient conduction-convective heat flow equation. Combined with transient analytical methods (Carslaw, Jaeger, Taniguchi), groundwater flow can be calculated more accurately. This method not only significantly improves the accuracy of groundwater flow calculation but also reduces the possibility of multiple solutions for other methods (numerical inversion methods), indirectly improving the accuracy of this method. It abandons the original steady-state assumption calculation method and, based on the transient analytical method, adds a heat source correction term for the influence of rainfall. Furthermore, it can calculate different groundwater flow rates based on rainfall intensities under multiple rainfall events over a period of time. The value is used to achieve accurate and reasonable superposition of heat source terms.

[0019] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0020] Figure 1 This is a comparison chart of the accuracy improvement rate of the present invention under different characteristic attenuation constants L; Figure 2 The curves showing the variation of the characteristic attenuation constant L for different rainfall intensities according to this invention; Figure 3 The curves showing the change in the accuracy improvement rate of this invention under different rainfall intensities are shown. Figure 4 The following are comparison charts of error analysis under different rainfall intensities: a) is a comparison chart of error analysis under a specific rainfall intensity of 30 mm / h; b) is a comparison chart of error analysis under different rainfall intensities. Figure 5 This is a comparison chart of the residuals obtained by fitting the temperature inversion data using the original formula and the present invention. Figure 6 The present invention employs a method for optimizing boundary conditions, and presents an error analysis comparison chart for a rainfall intensity of 25 mm / h. Detailed Implementation

[0021] The technical solution of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0022] Unless otherwise defined, the technical or scientific terms used in this invention shall have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0023] Example Please see Figures 1-6 This invention provides a correction method for estimating groundwater flow based on the influence of environmental factors on heat sources, comprising the following steps: Step 1: Construct a one-dimensional vertically homogeneous saturated porous medium hypothesis, and obtain temperature variation curves at various depths over time through field monitoring and laboratory measurements. T ( z,t Hourly rainfall intensity I ( t ), thermal diffusivity saturated hydraulic conductivity Ks Changes in moisture content .

[0024] This method is applicable to sites that meet the following conditions: ① one-dimensional vertical thermal transport, the soil can be regarded as a homogeneous isotropic medium, and the thermal property parameters ( ① ρc) is spatially uniformly distributed and does not change with depth; ② the study depth range remains saturated, ignoring the influence of the unsaturated zone (applicable to scenarios where the vadose zone above the groundwater level is thin or where surface water accumulates during rainfall). Based on this, the following basic parameters were obtained through field monitoring and laboratory measurements: 1. Temperature time series: A high-precision temperature sensor array is deployed vertically along the target site at typical depths of 0.05m, 0.1m, 0.2m, 0.3m, and 0.5m. Temperature data is continuously recorded at 1-10 minute intervals to obtain temperature variation curves over time at each depth. T ( z , t The sensor accuracy should be no less than ±0.02°C.

[0025] 2. Rainfall data: Synchronous rainfall data, including hourly rainfall intensity, is obtained from on-site automatic weather stations or nearby weather stations. I (t (Unit: mm / h) and rainwater temperature T rain ( t (Unit: °C). Rainfall intensity needs to be converted to SI units (m / s). ; 3. Soil thermal properties, thermal conductivity (W·m) -1 ·K -1 ): can be determined on-site by thermal probe or obtained by referring to tables based on soil texture (e.g., sandy soil 1.5-2.5, clay soil 1.0-1.8).

[0026] Volumetric heat capacity ρc (J·m -3 ·K -1 ): This can be determined by the thermal pulse method or calculated based on component weighting; the typical value range is 2.0 × -3.0× The thermal diffusivity is calculated from the above parameters: ; in, 4. Soil hydraulic parameters, saturated hydraulic conductivity Ks (m / s): can be obtained through field permeability tests or laboratory soil column tests; typical value for sandy soil. - loam - .

[0027] Moisture content change Δ (Dimensionless): Represents the increase in soil volumetric water content before and after rainfall. It can be measured by a soil moisture sensor or estimated based on rainfall intensity and soil water holding capacity. Based on experience, the value range is 0.1 to 0.4.

[0028] Step 2: Based on the assumption of a one-dimensional vertically homogeneous saturated porous medium, a basic heat transport equation is constructed. The one-dimensional vertical transient heat conduction-convection equation proposed by Stallman (1965) is used to describe the underground temperature field without rainfall disturbance. ; Step 3: Correct for rainfall disturbances by introducing a rainfall heat source term and improving surface boundary conditions.

[0029] The correction for the precipitation heat source term includes the following steps: S311. Introduce an explicit heat source term S on the right side of the equation in step two. i (z,t), assuming there are no other significant heat sources besides rainfall, we can deduce from the law of energy conservation: in, This indicates new sources and sinks of heat generated by rainfall. The density of water, The specific heat capacity of water, For the first i The rainfall intensity per unit time of a rainfall event. For the first i The temperature of the rainwater at a certain moment in this rainfall event. h ( z Let be the vertical distribution function of the heat source, using an exponential decay form: in, The characteristic attenuation depth represents the effective depth of the heat effect; The value can be taken in one of the following two ways: A reasonable value is directly taken based on experience (usually 0.05–0.3 m); the value proposed in this patent is based on thermal diffusivity. Changes in moisture content With peak rainfall intensity The expression is jointly determined. Assuming a short-duration heavy rainfall event... Controlled by peak intensity and If is a constant, then: Where Ipeak represents the maximum instantaneous intensity (m / s) of the rainfall event. The physical meaning of this formula is: the greater the thermal diffusivity and the greater the change in water content, the greater the thermal impact depth; the greater the rainfall intensity, the faster the heat is carried into deeper layers, resulting in a smaller attenuation depth. The parameter range selected in this experiment is: Take 0.5× ~3.0× , Using a value of 0.1 to 0.4, the typical value of L calculated is 0.05 to 0.3 m, which is consistent with field experimental observations.

[0030] S312. Integrating the equations, we obtain the modified heat transport equations. Substituting the heat source term from S311 into the basic equations, we obtain the one-dimensional heat transport governing equations considering rainfall corrections: This equation can be solved numerically using the finite difference method or the finite element method. If an explicit Euler scheme is used, the time step must satisfy the CFL stability condition (condition ⑧) to ensure computational convergence. ; , .

[0031] S313, Using measured temperature time series Inversion of unknown parametersq Assuming that during the analysis period q It can be considered a constant (or a piecewise constant) and solved using a global optimization algorithm (such as differential evolution). The inversion objective is to minimize the root mean square error (RMSE) between the model-calculated temperature and the observed temperature.

[0032] S3131. Construct a forward model. Using the modified heat transport equation as the forward model, numerically discretize it using the finite difference method (or finite element method). The computational domain of the model is depth z∈[0, Lz The spatial step size Δz is 0.01m; the time step Δt must satisfy the stability condition of the explicit scheme, while the implicit scheme is unconditionally stable, ensuring computational accuracy and convergence; the input data includes: the initial values ​​of the measured temperature at each depth. T ( z ,0) (usually the temperature profile at the first observation time is taken); surface boundary conditions of the time series; parameters related to the precipitation heat source term: characteristic attenuation depth L (Calculated according to formula (20) or taken from empirical values), thermal properties of water Soil thermal diffusivity wait.

[0033] S3132. Define the objective function, using the measured temperature time series. To achieve the calibration objective, the root mean square error (RMSE) is chosen as the objective function to measure the simulated values. Degree of deviation from observed values: in, For sensor depth, such as 0.05m, 0.10m, 0.20m, 0.40m, N The total number of observations at all depths and all time points; the smaller the objective function value, the better the simulation results match the measured data.

[0034] S3133, Optimization algorithm selection and parameter range setting. v For the scalar to be inverted, a global optimization algorithm can be used to search for the optimal solution within a preset physical reasonable range. The differential evolution algorithm is used, with the following parameter settings: the population size is 15–20 times the number of parameters to be inverted, the maximum number of iterations is set to 100–200, and the mutation factor and crossover probability use their default values. v The search scope is set based on the hydrogeological background of the study area, and usually takes... - m / s.

[0035] S3134. Inversion Execution and Result Output: In each iteration, the algorithm calls the forward model to calculate the current... vThe value corresponds to the simulated temperature, and the objective function value is calculated. The output value minimizes the RMSE. v The values ​​are used as inversion results. To assess the uncertainty of the inversion, multiple independent runs or posterior analysis can be performed, and statistics can be collected. v The mean and standard deviation.

[0036] S314. Method Validation and Accuracy Evaluation: To verify the superiority of this method, the following two types of validation are required: validation with synthetic data and validation with measured data.

[0037] The synthetic data uses the modified equation of this invention to generate synthetic temperature data including rainfall events, and adds measurement noise of ±0.02℃. The flow rate is inverted using the heat-free term and the method of this invention, respectively, and the relative error is calculated. The measured data verifies the combination of real measurement conditions and real weather conditions by comparing the degree of fit between the original formula and the improved formula with the real values.

[0038] Step 3, the improvement of surface boundary condition correction, includes the following steps: S321. Improved Surface Boundary Conditions: To characterize the transient disturbances of rainfall events on the surface thermal state, this invention proposes a boundary condition based on surface heat flux balance, replacing the traditional fixed temperature boundary. It considers the energy budget at the surface (z=0), including air convection heat exchange and rainfall heat flux. Based on energy conservation, the following is derived: in, H The convective heat transfer coefficient, = T (0, t The surface temperature is represented by the heat flux from the surface into the ground via conduction (left side), which is equal to the sum of the convective heat exchange between the surface and the air (first item on the right) and the net heat flux from rainfall (second item on the right). When cold rain occurs, the second item becomes negative, causing the surface temperature to drop. S322. Complete Mathematical Model and Numerical Solution: Combine the formula in S321 with the equation in step two, and add the boundary conditions and initial conditions to form a complete mathematical model. This model can be solved numerically using either the finite difference method or the finite element method. When using the finite difference method, the depth z is uniformly discretized into N nodes with a node spacing Δz and a time step Δt. These nodes must satisfy the stability conditions of the explicit scheme (unconditional stability is not possible with the implicit scheme). At each time step, the internal nodes are first calculated based on the current temperature field, and then the surface temperature is updated using the discretized boundary conditions. This formula is derived from the first-order difference approximation of formula (23) and can be directly used for explicit time advancement.

[0039] S323. Inverting vertical groundwater flux using measured temperature time series. Inversion of unknown parameters And possibly unknown H Assuming that during the analysis period It can be considered a constant (or a piecewise constant), and solved using a global optimization algorithm (such as differential evolution). The inversion objective is to minimize the root mean square error (RMSE) between the model-calculated temperature and the observed temperature. The parameter search range is set based on the hydrogeological background of the study area, and the boundary conditions during the optimization process are... H It can be used for inversion or fixed based on experience.

[0040] S324. Based on actual observations of rainfall intensity, daily soil temperature variation, seasonal cooling trend, and station characteristics, simulated data are generated to compare the original formula and the improved formula with the observed data.

[0041] Detailed derivation process: The transient conduction-convective heat flow equation in the Cartesian z-direction (depth) (Stallman, 1965): Based on this, add a heat source item: in, i Indicates the number of rainfall events. This indicates new sources and sinks of heat generated by rainfall.

[0042] Now let's derive... The expression for is defined as the rate of increase in net heat due to precipitation infiltration per unit time and unit volume. We assume a tiny soil cubic unit with a base area of... Δy, height Δz. Rainwater infiltrates vertically through it. The volume of rainwater flowing into the unit per unit time: in, For the first i Volumetric flow rate of the rainfall event For the first i The rainfall intensity per unit time of a rainfall event. A Given the base area, the heat energy carried by the rainwater in this portion is the heat energy flowing in: in, For the incoming heat energy, The density of water, The specific heat capacity of water, For the first i The temperature of the rainwater at a certain moment during a rainfall event.

[0043] Assuming that after rainwater undergoes a complete heat exchange with the soil instantaneously, the water temperature when it flows out of the unit is equal to the soil temperature. T ( z , t Therefore, the heat energy carried away during outflow is: in, For the outflowing heat energy, Since the soil temperature is constant, the net increase in heat is: because The definition is the increase in heat per unit time and per unit volume. Therefore, we divide the net increase in heat by the volume of the soil unit to obtain: Therefore, With 1 / Δ z They are directly proportional; if we subdivide the soil unit infinitely, that is, when Δ z Approaching positive infinity, It will approach infinity, which physically corresponds to the exchange of heat within an infinitely thin layer, which is precisely what Dirac... function This describes the situation. Therefore, we define a new function. h ( z Its physical meaning is: at depth z The proportion of heat released per unit volume of infiltrated water. h ( z Satisfying the normalization condition Therefore, we can conclude that: Let us now explain in detail the formula (8) In a simple and rough assumption, when rain falls to the ground, all the heat exchange between the rainwater and the soil occurs at the surface ( z =0) This is completed instantly on an infinitely thin layer, and the rainwater, carrying a new temperature, immediately and completely enters the soil. Afterwards, as it flows in the soil, it no longer exchanges new heat with the surroundings. =Dirac function Right now: Of course, this assumption is too idealistic. In actual rainfall, heat is not released at a single point on the surface, but continuously along the path of the infiltration front, and the intensity of the release decreases with depth. Therefore, we introduce an exponential decay function: in, The characteristic attenuation depth is represented by , and its larger value indicates a deeper extent of heat influence. A is the normalization coefficient. To satisfy the normalization condition, Will h ( z Substituting into equation (11) yields: visible: Therefore, equation (10) can be expressed as: Now let's derive... The conventional expression assumes that all rainfall infiltrates (no runoff, no surface water) and that the soil moisture profile exhibits an exponential decay distribution after the rainfall ends: Define a feature time scale τ This indicates the transfer of heat from the Earth's surface to depth. L The time required. τ It should be determined by both thermal diffusion and infiltration. Based on dimensional analysis, a reasonable form for the unit of length can be derived: in, denoted as thermal diffusivity.

[0044] Now, let's determine the time scale. τ ,because τ This reflects how the water flow carries heat to the surface. L How long does it take to reach the depth? The water flow velocity can be measured using the Darcy flow rate. q It indicates, but q This is precisely the target we need to invert; it cannot appear in [the following context]. L In the definition, therefore we use (Saturated hydraulic conductivity) is used to approximate q The magnitude is because the surface is nearly saturated during rainfall, and the infiltration rate is close to [a certain value]. Therefore, we have: Substituting into equation (16) yields: After simplification by squaring, we get: Add Factors affecting (moisture content change) What should be affected is the actual infiltration rate. Instead Ks According to the Green-Ampt infiltration model: use Substitute (18) You will then receive: (19) After testing and correction, it passed To replace saturated hydraulic conductivity This better reflects the physical reality that infiltration is controlled by rainfall intensity under heavy rainfall conditions, namely: According to equations (1), (2), (8), and (14), we can obtain: However, in general calculations, Difficult to obtain precisely, usually If we can directly determine the formula, it becomes: We found that when the characteristic decay depth L→0, the heat source term in equation (22) When the expression approaches the Dirac delta function δ(z), equation (22) becomes: This equation describes the concentration at the Earth's surface ( z The influence of an instantaneous heat source with thickness ε on the temperature field (0 ≤ ε) is considered. To transform this into a boundary condition, a thin layer of thickness ε near the Earth's surface (0 ≤ ε) is also considered. z ≤ ε Energy balance analysis was performed on the thin layer.

[0045] Let the external heat flux at the Earth's surface (positive downwards) be... This is typically provided by air convection heat exchange: in Tsurface = T (0, t ( ) represents the surface temperature. H The convective heat transfer coefficient. The lower boundary of the thin layer ( z = ε The internal heat flow at point (positive downwards) is: Thin layer due to δ The total heat generated by the source term is: Neglecting the change in internal energy of the thin layer (which tends to zero as ε→0), energy conservation gives: make ε →0, then Substituting, we get: After processing, the boundary conditions for surface heat flux are obtained: How to add a heat source item: Step 1: This method is applicable to sites that meet the following conditions: ① One-dimensional vertical thermal transport, the soil can be considered a homogeneous and isotropic medium, and the thermal properties (κ, ρc) are spatially uniformly distributed and do not change with depth; ② The study depth is always in a saturated state, and the influence of the unsaturated zone is ignored (applicable to scenarios where the vadose zone above the groundwater level is thin or where surface water accumulates during rainfall). Based on these conditions, the following basic parameters are obtained through field monitoring and laboratory measurements: 1. Temperature Time Series: A high-precision temperature sensor array is vertically deployed at the target site, with typical depths of 0.05m, 0.1m, 0.2m, 0.3m, and 0.5m. Temperature data is continuously recorded at intervals of 1-10 minutes to obtain the temperature change curve T(z,t) at each depth over time. The sensor accuracy should be no less than ±0.02℃.

[0046] 2. Rainfall data: Synchronous rainfall data is obtained from on-site automatic weather stations or nearby weather stations, including hourly rainfall intensity I(t) (unit: mm / h) and rainfall temperature Train(t) (unit: °C). Rainfall intensity needs to be converted to SI units (m / s). 3. Soil thermal physical properties, thermal conductivity κ (W·m) -1 ·K -1 ): can be determined on-site by thermal probe or obtained by referring to tables based on soil texture (e.g., sandy soil 1.5-2.5, clay soil 1.0-1.8).

[0047] Volumetric heat capacity ρc (J·m -3 ·K -1 ): This can be determined by the thermal pulse method or calculated based on component weighting; the typical value range is 2.0 × ~3.0× The thermal diffusivity is calculated from the above parameters: in, 4. Soil hydraulic parameters, saturated hydraulic conductivity Ks (m / s): can be obtained through field permeability tests or laboratory soil column tests; typical value for sandy soil. - loam - .

[0048] Change in soil moisture content Δθ (dimensionless): represents the increase in soil volumetric moisture content before and after rainfall. It can be measured by a soil moisture sensor or estimated based on rainfall intensity and soil water holding capacity. Based on experience, the value ranges from 0.1 to 0.4.

[0049] Step 2: Based on the assumption of a one-dimensional vertically homogeneous saturated porous medium (conditions ① and ②), construct the basic heat transport equation, and use the one-dimensional vertical transient heat conduction-convection equation proposed by Stallman (1965) to describe the underground temperature field without rainfall disturbance: Where q is the vertical Darcy velocity (m / s, positive downwards). , The volumetric heat capacity of water (J·) · ), take 4.18× This equation neglects thermal dispersion and source / sink terms, and is applicable to homogeneous saturated porous media.

[0050] Step 3: To quantitatively describe the disturbance of the temperature field by rainfall, an explicit heat source term S_i(z,t) is introduced on the right-hand side of the basic equation. Assuming there are no other significant heat sources besides rainfall (condition ③), the following can be derived based on energy conservation: Where h(z) is the vertical distribution function of the heat source, which adopts an exponential decay form: in, The characteristic attenuation depth (m) represents the effective depth of the heat effect. The value can be taken in one of the following two ways: A reasonable value is directly taken based on experience (usually 0.05–0.3 m); the method proposed in this patent, which combines thermal diffusivity α, water content change Δθ, and peak rainfall intensity, is used instead. The expression is jointly determined. Assuming a short-duration heavy rainfall event... If the peak intensity is controlled and Δθ is constant (conditions ⑤ and ⑥), then: Where Ipeak represents the maximum instantaneous intensity (m / s) of the rainfall event. The physical meaning of this formula is: the greater the thermal diffusivity and the greater the change in water content, the greater the thermal impact depth; the greater the rainfall intensity, the faster the heat is carried into deeper layers, leading to a decrease in the attenuation depth. The parameter range selected in this experiment: α is 0.5 × ~3.0× With Δθ ranging from 0.1 to 0.4, the typical value of L calculated is 0.05 to 0.3 m, which is consistent with field experimental observations.

[0051] Step 4: Integrate to obtain the modified heat transport equation. Substitute the heat source term from Step 3 into the basic equation to obtain the one-dimensional heat transport governing equation considering rainfall correction: This equation can be solved numerically using the finite difference method or the finite element method. If an explicit Euler scheme is used, the time step must satisfy the CFL stability condition (condition ⑧) to ensure computational convergence. ; Where Δz is the spatial grid spacing, which is generally taken as 0.01m.

[0052] Step 5: Utilize the measured temperature time series Invert the unknown parameter q. Assuming q can be considered a constant (or piecewise constant) during the analysis period, a global optimization algorithm (such as differential evolution) is used to solve it (condition ⑦). The inversion objective is to minimize the root mean square error (RMSE) between the model-calculated temperature and the observed temperature.

[0053] 1. Construct a forward model, using the modified heat transport equation as the forward model, and perform numerical discretization using the finite difference method (or finite element method). The computational domain of the model is depth z∈[0,Lz], with a spatial step size Δz0.01m; the time step size Δt must satisfy the stability condition of the explicit scheme (if the implicit scheme is used, it is unconditionally stable) to ensure computational accuracy and convergence. Input data includes: the initial value of the measured temperature T(z,0) at each depth (usually the temperature profile at the first observation time); the surface boundary conditions of the time series; parameters related to the precipitation heat source term: characteristic attenuation depth L (calculated according to formula (20) or taken as an empirical value), and the thermal properties of water. Soil thermal diffusivity wait.

[0054] 2. Define the objective function, using the measured temperature time series. To achieve the calibration objective, the root mean square error (RMSE) is chosen as the objective function to measure the simulated values. Degree of deviation from observed values: in, Let N be the sensor depth (e.g., 0.05m, 0.10m, 0.20m, 0.40m), and N be the total number of observations at all depths and time points. The smaller the objective function value, the better the simulation results match the measured data.

[0055] 3. Optimization Algorithm Selection and Parameter Range Setting: Since v is a scalar to be inverted, a global optimization algorithm can be used to search for the optimal solution within a pre-defined physical reasonable range. Differential Evolution is used because it is robust to nonlinear problems and does not require gradient calculation. The algorithm parameters are set as follows: the population size is 15–20 times the number of parameters to be inverted, the maximum number of iterations is set to 100–200, and the mutation factor and crossover probability use default values. The search range of v is set according to the hydrogeological background of the study area, typically taking [value missing]. ~ m / s (corresponding to 8.6–86 mm / day). If there is prior knowledge of the recharge in the study area, the range can be appropriately narrowed.

[0056] 4. Inversion Execution and Result Output: In each iteration, the algorithm calls the forward model to calculate the simulated temperature corresponding to the current v value and calculates the objective function value. The final output is the v value that minimizes the RMSE as the inversion result. To assess the inversion uncertainty, multiple independent runs or posterior analysis can be performed to statistically analyze the mean and standard deviation of v.

[0057] Step Six: Method Validation and Accuracy Evaluation To verify the superiority of this method, the following two types of validation are required: 1. Synthetic data verification, based on known real traffic. Under these conditions, synthetic temperature data incorporating rainfall events is generated using the modified equation of this invention, with measurement noise of ±0.02℃ added. Flow rates are inverted using both the conventional method (without a heat source term) and the method of this invention, and the relative errors are calculated. The comparison shows that when rainfall intensity is ≥5 mm / h, the error of the method of this invention is reduced by more than 50% (see...). Figure 3 , Four ).

[0058] 2. Verification with measured data: Based on the actual measurement of the geothermal gradient in Zhonggu Village, Kangding, and the real weather conditions at the same time, the results of comparing the fit between the original formula and the improved formula with the actual values ​​show that the improved formula can improve the performance.

[0059] Specific steps for modifying boundary conditions: Steps one and two are the same as above.

[0060] Step 3: Improved Surface Boundary Conditions. To characterize the transient disturbances of rainfall events to the surface thermal state, this invention proposes a boundary condition based on surface heat flux balance, replacing the traditional fixed temperature boundary. Considering the energy budget at the surface (z=0), including air convection heat exchange and rainfall heat flux, the following is derived based on energy conservation: in, =T(0,t) represents the Earth's surface temperature. The physical meaning of this formula is: the heat flow from the Earth's surface into the ground via conduction (left side) equals the sum of the convective heat exchange between the surface and the air (first term on the right) and the net heat flux from rainfall (second term on the right). When cold rain occurs, the second term becomes negative, causing the surface temperature to drop.

[0061] Step 4: Complete Mathematical Model and Numerical Solution Combine equation (29) with the internal equation (1), and add the lower boundary conditions and initial conditions to form a complete mathematical model: This model can be solved numerically using either the finite difference method or the finite element method. When using the finite difference method, the depth z is uniformly discretized into N nodes with a node spacing Δz and a time step Δt. These nodes must satisfy the stability conditions of the explicit scheme (unconditional stability is not possible with the implicit scheme). At each time step, the internal nodes are first calculated based on the current temperature field, and then the surface temperature is updated using the discretized boundary conditions. This formula is derived from the first-order difference approximation of formula (23) and can be directly used for explicit time advancement.

[0062] Step 5: Invert vertical groundwater flux using measured temperature time series. Inversion of unknown parameters (And possibly unknown H). Assuming that during the analysis period... This can be treated as a constant (or piecewise constant) and solved using a global optimization algorithm (such as differential evolution). The inversion objective is to minimize the root mean square error (RMSE) between the model-calculated temperature and the observed temperature. The parameter search range is set according to the hydrogeological background of the study area, for example, q∈[ , ]m / s, H∈[5,25]W· · During the optimization process, H in the boundary conditions can participate in the inversion or be fixed based on experience.

[0063] Step Six: Based on the actual observation range of the Heihe Huazhaizi Station—including rainfall intensity (≤40mm / h, temperature difference ≤8℃), daily soil temperature variation, seasonal cooling trend, and other real station characteristics—simulated data is generated. The original and improved formulas are then compared with the observed data. It can be seen that when the rainfall intensity is around 25mm / h, the accuracy improvement can reach 28.6%. Figure 6 .

[0064] Therefore, this invention employs the aforementioned correction method for calculating groundwater flow considering the influence of environmental factors on heat sources. Based on the transient conduction-convective heat flow equation in the Cartesian z-direction (depth), and incorporating the Green-Ampt influence on soil infiltration capacity limitations, it is the first to add precipitation factors to the transient conduction-convective heat flow equation. Combined with transient analytical methods (Carslaw, Jaeger, Taniguchi), groundwater flow can be calculated more accurately. This method not only significantly improves the accuracy of groundwater flow calculation but also reduces the possibility of multiple solutions for other methods (numerical inversion methods), indirectly improving the accuracy of this method. It abandons the original steady-state assumption calculation method and, based on the transient analytical method, adds a heat source correction term for the influence of rainfall. Furthermore, it can calculate different groundwater flow rates based on rainfall intensities under multiple rainfall events over a period of time. The value is used to achieve accurate and reasonable superposition of heat source terms.

[0065] Hydrological-thermal dynamic coupling characteristics: Traditional thermal tracing methods (such as inversion models based on the Stallman equation) neglect transient thermal disturbances during rainfall events and assume soil thermophysical parameters to be constants, failing to dynamically correlate temperature changes with the real-time evolution of groundwater vertical flux. This invention proposes an explicit rainfall heat source term and its quantitative relationship with soil physical parameters, combined with a finite difference numerical platform, to achieve, for the first time, real-time dynamic coupling of heat transfer, convection, and rainfall disturbances. Based on the one-dimensional transient heat conduction-convection equation, a method using thermal diffusivity is introduced. Changes in moisture content With peak rainfall intensity The characteristic attenuation depth is determined by the joint factors The model dynamically correlates the real-time responses of rainfall intensity, soil hydraulic parameters, and temperature field. Under typical sandy loam conditions, with rainfall intensity of 30–50 mm / h and soil-air temperature difference >8℃, compared with traditional models that ignore rainfall disturbance, this method reduces the root mean square error of temperature fitting at a depth of 5–20 cm during rainfall by 45%–70%, significantly improving the reliability of groundwater-surface water exchange flux estimation under heavy rainfall conditions.

[0066] Boundary condition optimization: Traditional thermal tracing methods (such as inversion models based on the Stallman equation) typically employ fixed surface temperature boundary conditions, failing to reflect the transient disturbances of surface thermal state caused by rainfall events, leading to significantly increased inversion errors under heavy rainfall conditions. This invention proposes an improved surface heat flux boundary condition, dynamically coupling air convection heat exchange with rainfall heat flux, achieving for the first time a real-time mapping of surface thermal disturbances caused by rainfall events. Based on the one-dimensional transient heat conduction-convection equation, the upper boundary condition is modified. This boundary condition directly couples the atmosphere-surface-rainfall heat exchange process, eliminating the need to introduce empirical internal heat source parameters, resulting in clear physical meaning and simplified parameters.

[0067] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A correction method for estimating groundwater flow rate considering the influence of environmental factors on heat sources, characterized in that, Includes the following steps: Step 1: Construct a one-dimensional vertically homogeneous saturated porous medium hypothesis, and obtain temperature variation curves at various depths over time through field monitoring and laboratory measurements. T ( z,t Hourly rainfall intensity I ( t ), thermal diffusivity saturated hydraulic conductivity Ks Changes in moisture content ; Step 2: Based on the assumption of a one-dimensional vertically homogeneous saturated porous medium, a basic heat transport equation is constructed. The one-dimensional vertical transient heat conduction-convection equation proposed by Stallman (1965) is used to describe the underground temperature field without rainfall disturbance. ； Step 3: Correct for rainfall disturbances by introducing a rainfall heat source term and improving surface boundary conditions.

2. The correction method for calculating groundwater flow based on heat sources considering environmental factors, as described in claim 1, is characterized in that: In step one, one-dimensional vertical thermal transport must be satisfied. The soil is regarded as a homogeneous isotropic medium with thermal properties uniformly distributed in space and not changing with depth. The soil is always in a saturated state within the study depth range, and the influence of the unsaturated zone is ignored.

3. The correction method for calculating groundwater flow rate considering the influence of environmental factors on heat sources according to claim 2, characterized in that, Step three, which involves correcting for the precipitation heat source term, includes the following steps: S311. Introduce an explicit heat source term S on the right side of the equation in step two. i (z,t), assuming there are no other significant heat sources besides rainfall, we can deduce from the law of energy conservation: in, This indicates new sources and sinks of heat generated by rainfall. The density of water, The specific heat capacity of water, For the first i The rainfall intensity per unit time of a rainfall event. For the first i The temperature of the rainwater at a certain moment in this rainfall event. h ( z Let be the vertical distribution function of the heat source, using an exponential decay form: in, The characteristic attenuation depth represents the effective depth of the heat effect; S312. Integrating the equations, we obtain the modified heat transport equations. Substituting the heat source term from S311 into the basic equations, we obtain the one-dimensional heat transport governing equations considering rainfall corrections: Numerical solutions are obtained using the finite difference method or the finite element method, employing an explicit Euler scheme. The time step must satisfy the CFL stability condition to ensure computational convergence. ； Where Δz is the spatial grid spacing, which is generally taken as 0.01m; S313, Using measured temperature time series Inversion of unknown parameters q The global optimization algorithm is used to solve the problem, and the inversion objective is to minimize the root mean square error (RMSE) between the model-calculated temperature and the observed temperature. S314. Method validation and accuracy evaluation: Validation is performed using synthetic data and measured data.

4. The correction method for calculating groundwater flow rate considering the influence of environmental factors on heat sources according to claim 3, characterized in that: S311 Controlled by peak intensity and If is a constant, then we have: in, Where I is the thermal diffusivity, and Ipeak is the maximum instantaneous intensity during a rainfall event. Take a value of 0.1 to 0.

4.

5. The correction method for calculating groundwater flow rate considering the influence of environmental factors on heat sources according to claim 4, characterized in that, S313 includes the following steps: S3131. Construct a forward model. Using the modified heat transport equation as the forward model, numerical discretization is performed using the finite difference method. The computational domain of the model is depth z∈[0, Lz The spatial step size Δz is 0.01m; the time step size Δt must satisfy the stability condition of the explicit scheme, while the implicit scheme is unconditionally stable, thus ensuring computational accuracy and convergence. S3132. Define the objective function, using the measured temperature time series. To achieve the calibration objective, the root mean square error (RMSE) is chosen as the objective function to measure the simulated values. Degree of deviation from observed values: in, For sensor depth, N The total number of observations at all depths and all time points; S3133, Optimization algorithm selection and parameter range setting. v For the scalar to be inverted, the differential evolution algorithm is used, with the following parameters set: the population size is 15–20 times the number of parameters to be inverted, the maximum number of iterations is set to 100–200, and the mutation factor and crossover probability are set to default values. v Usually taken - m / s; S3134. Inversion Execution and Result Output: Forward Model Calculation of Current... v The value corresponds to the simulated temperature, and the objective function value is calculated. The output value minimizes the RMSE. v The value is used as the inversion result.

6. The correction method for calculating groundwater flow rate considering the influence of environmental factors on heat sources according to claim 5, characterized in that, In S314, the synthetic data uses the modified equation of the present invention to generate synthetic temperature data containing rainfall events, and adds measurement noise of ±0.02℃. The flow rate is inverted using the heat-free term and the method of the present invention, respectively, and the relative error is calculated. The measured data verifies the combination of real measurement conditions and real weather conditions by comparing the degree of fit between the original formula and the improved formula with the real values.

7. The correction method for calculating groundwater flow rate considering the influence of environmental factors on heat sources according to claim 2, characterized in that, Step 3, the improvement of surface boundary condition correction, includes the following steps: S321. Improved surface boundary conditions, considering the energy budget at the surface (z=0), including air convection heat exchange and precipitation heat flux, derived according to energy conservation: in, H The convective heat transfer coefficient, = T (0, t ) represents the surface temperature. The heat flow from the surface into the ground via conduction is equal to the sum of the convective heat exchange between the surface and the air and the net heat flux from rainfall. S322. Complete Mathematical Model and Numerical Solution: Combine the formula in S321 with the equation in step two, and add the boundary conditions and initial conditions to form a complete mathematical model. ； S323. Inverting vertical groundwater flux using measured temperature time series. Inversion of unknown parameters The inversion objective is to minimize the root mean square error (RMSE) between the model-calculated temperature and the observed temperature. ； S324. Based on actual observations of rainfall intensity, daily soil temperature variation, seasonal cooling trend, and station characteristics, simulated data are generated to compare the original formula and the improved formula with the observed data.

8. The correction method for calculating groundwater flow rate considering the influence of environmental factors on heat sources according to claim 7, characterized in that, S321 includes the following steps: The external heat flux at the Earth's surface is This is typically provided by air convection heat exchange: in, Tsurface = T (0, t ) represents the surface temperature. H The convective heat transfer coefficient is the lower boundary of the thin layer. z = ε The internal heat flow at that location is: Thin layer due to δ The total heat generated by the source term is: Ignoring the change in internal energy of the thin layer as ε approaches zero, the law of conservation of energy states: Let ε→0, then Substituting, we get: The boundary conditions for surface heat flux are obtained as follows: 。 9. The correction method for calculating groundwater flow based on heat sources considering environmental factors, as described in claim 8, is characterized in that: In S322, the finite difference method or finite element method is used for numerical solution. At each time step, the internal nodes are first calculated based on the current temperature field, and then the surface temperature is updated through discretized boundary conditions. 。

Images