Moisture dynamic prediction method based on layered soil improved Green-Ampt redistribution model

Abstract

The invention discloses a dynamic moisture prediction method based on a layered soil improved Green-Ampt redistribution model. The method comprises the following steps: acquiring layered soil parameters and rainfall data of a research object; rainfall segmentation processing: for each calculation time period, based on a composite multi-layer soil parameterization method and rainfall data, determining the time for starting water accumulation at the top of the soil in a calculation time step length; comparing the time for starting water accumulation at the top of the soil with the end time for calculating the time step length, and calculating the increment of the volumetric water content of the soil by adopting different modes according to the relationship between the two times; updating a wet front position based on an improved Green-Ampt redistribution model of layered soil, and realizing dynamic simulation of infiltration and redistribution processes in the multilayer soil in a calculation time period; and calculating runoff production and accumulated infiltration amount data of the rainfall data under the condition of determining the soil parameters, and ending calculation. According to the method, the simulation precision of the layered soil hydrological process is remarkably improved while the calculation efficiency is kept.

Description

Technical Field

[0001] The present invention belongs to the technical field of hydrological models and soil moisture dynamics simulation, and particularly relates to a method for predicting water dynamics of an improved Green-Ampt redistribution model based on layered soil. Background Art

[0002] The accurate simulation of soil moisture dynamics is the core basis for hydrological prediction, agricultural irrigation, and flood disaster prevention and control. Traditional hydrological models generally assume that the soil is a homogeneous medium, but in the actual natural environment, layered soil structures widely exist (such as alternating layers of sand and clay), and there are significant differences in the permeability, water holding capacity, and capillary pressure of each layer. When existing models with the assumption of homogeneous soil simulate the water migration in layered soil, due to the lack of consideration of the sudden change in hydraulic characteristics between layers, there are systematic biases in the prediction of infiltration rate, wetting front advancement depth, and runoff generation.

[0003] The traditional Green-Ampt (GA) model simplifies the wetting front as a discrete object, avoiding complex numerical calculations, but it is only applicable to single-layer homogeneous soil and cannot characterize the water redistribution process in multiple rainfall events. Although the subsequent improved GAR model introduces a multi-wetting front mechanism, there are still significant limitations in the scenario of layered soil: when the wetting front advances across layers, the model does not dynamically adjust the continuity of interlayer parameters, resulting in distorted calculation of water flux at the interface; in addition, the existing methods have insufficient response ability to complex rainfall pulses and are difficult to simulate the secondary infiltration and runoff superposition effects caused by intermittent precipitation. Numerical models based on the solver of the Richard equation can improve the accuracy, but they rely on dense grid discretization and iterative calculations, with low computational efficiency and are difficult to apply to large-scale watershed simulations or real-time disaster warnings.

[0004] The existing research on layered soil models mostly focuses on theoretical expansion, and there are still multiple challenges in practical applications. For example, some models achieve multi-layer characterization by fixing parameter amplification, but the uncertainty of parameter estimation limits their universality; some other models assume a constant surface ponding depth and cannot flexibly adapt to the dynamic changes in rainfall intensity. In addition, the lack of a multi-wetting front merging algorithm leads to mass conservation deviation at the cross-layer interface, further affecting the stability of long-term simulations, highlighting the deficiencies of the existing technologies. Summary of the Invention

[0005] Object of the Invention: The object of the present invention is to provide a method for predicting water dynamics of an improved Green-Ampt redistribution model based on layered soil, which is used to accurately simulate the water infiltration, redistribution, and runoff generation processes in multi-layer soil.

[0006] Technical Solution: The method for predicting water dynamics of the improved Green-Ampt redistribution model based on layered soil according to the present invention includes the following steps:

[0007] S1. Obtain the layered soil parameters and rainfall data of the research object;

[0008] S2. Process the rainfall in segments. For each calculation period, based on the composite multi-layer soil parameterization method and precipitation data, determine the time when water starts to accumulate at the top of the soil within the calculation time step;

[0009] S3. Compare the time t when water starts to accumulate at the top of the soil p with the end time t + △t of the calculation time step. If the end time of the calculation time step is not greater than the time when water starts to accumulate at the top of the soil, go to S4; otherwise, go to S5;

[0010] S4. When there is no water accumulation at the top of the soil, calculate the increment of soil volume water content according to the rainfall magnitude and the calculation time step △t, and go to S6;

[0011] S5. When there is water accumulation at the top of the soil during the calculation period, within the range [t, t p from the start time of the calculation to the time when water accumulation occurs, calculate the increment of soil volume water content according to the rainfall magnitude and the calculation time step t p - t. After water accumulation occurs, within the range [t

[0012] , t + △t] from the time when water accumulation occurs to the end of the calculation period, calculate the surface runoff and the increment of soil volume water content according to the rainfall magnitude, the calculation time step t + △t - t p , and the soil volume water content, and go to S6; p S6. Update the wetting front position based on the improved Green-Ampt redistribution model of layered soil to realize the dynamic simulation of the infiltration and redistribution processes in the multi-layer soil within the calculation period [t, t + △t];

[0013] S7. Judge whether the end time t + △t is greater than the rainfall end time. If it is greater, go to S8; otherwise, return to S2;

[0014] S8. Statistically calculate the runoff and cumulative infiltration data of the rainfall data under the determined soil parameters, and the calculation ends.

[0015] Furthermore, in step S1, the soil is vertically divided into N layers, and the hydraulic parameters are independently defined for each layer, including hydraulic conductivity, saturated hydraulic conductivity parameter of the maximum infiltration capacity, saturated soil volume water content when the soil is completely saturated, initial soil volume water content at the initial simulation moment, residual soil volume water content of the soil, soil thickness, and rainfall rate.

[0016] Furthermore, the composite multi-layer soil parameterization method in step S2 is as follows:

[0017] Furthermore, the composite multi-layer soil parameterization method in step S2 is as follows:

[0018] Dynamically calculate the composite saturated hydraulic conductivity K according to the current soil layer number n where the wetting front is located, the thickness of each layer, and the thickness of each layer s,c Reflecting the comprehensive water conductivity of multi-layer soil, K s,c The calculation formula is:

[0019]

[0020] Among them, Z is the depth reached by the wetting front, K s,i is the saturated hydraulic conductivity parameter of the maximum infiltration capacity of the i-th layer, K s,n is the saturated hydraulic conductivity parameter of the maximum infiltration capacity of the n-th layer, T i is the thickness of the i-th layer of soil.

[0021] Furthermore, determining the starting water accumulation time t at the top of the soil within the calculated time step in step S2 p , includes:

[0022] First, calculate the comprehensive capillary driving force G(θ b , θ) of the wetting front. The formula is:

[0023]

[0024] Among them, θ represents the soil volumetric water content at the wetting front, θ b represents the soil volumetric water content below the wetting front, ψ and ψ b values respectively correspond to the capillary pressures at the soil volumetric water contents θ and θ b , K(ψ) is the hydraulic conductivity of the capillary pressure ψ, K s is the saturated hydraulic conductivity;

[0025] For a certain position Z where the wetting front is located in the n-th layer, the soil infiltration capacity f r is expressed as:

[0026]

[0027] Among them, K s,n is the saturated hydraulic conductivity parameter of the maximum infiltration capacity of the n-th layer, θ s,n is the saturated soil volumetric water content of the n-th layer, θ i,n is the initial soil volumetric water content of the n-th layer, G(θ i,n , θ s,n ) is the comprehensive capillary driving force corresponding to the saturated volumetric water content at the wetting front, I is the cumulative infiltration amount up to the calculation moment, K s,c is the composite saturated hydraulic conductivity;

[0028] If the rainfall rate R is less than the composite saturated hydraulic conductivity K s,c, there will be no ponding within the calculated time step, and the ponding time t p is greater than the end time t + Δt of the calculated time step;

[0029] If the rainfall rate R is greater than the composite saturated hydraulic conductivity K s,c but less than the soil infiltration capacity f r , there may be ponding within the calculated time step [t, t + Δt]. At this time, it is first necessary to calculate the corresponding ponding occurrence time t according to the following formula and the composite saturated hydraulic conductivity K at this time step s,c : p,n :

[0030]

[0031] where I t is the cumulative infiltration amount up to the start time t of the calculation;

[0032] Finally, if the rainfall rate R is greater than the soil infiltration capacity f r , ponding will directly appear at the beginning stage within the calculated time step [t, t + Δt].

[0033] Furthermore, in step S4, calculating the increment of soil volumetric water content according to the rainfall magnitude and the calculated time step Δt includes:

[0034] If the ponding start time t at the soil top p is greater than the end time t + Δt of the calculated time step, there will be no ponding within the calculated time step. At this time, the increment of the cumulative infiltration amount ΔI is calculated by the following formula:

[0035] ΔI = R·Δt

[0036] where R is the rainfall rate.

[0037] Furthermore, in step S5, calculating the increment of soil volumetric water content according to the rainfall magnitude, the calculated time step t + Δt - t p and the soil volumetric water content includes:

[0038] There will be no ponding within the range [t, t p from the start time of the calculation to the ponding occurrence time. At this time, the increment of the cumulative infiltration amount ΔI is calculated by the following formula:

[0039] ΔI = R·(t p - t)

[0040] where R is the rainfall rate;

[0041] For the calculation period [t p , t + Δt], at this time the rainfall rate R is greater than the soil infiltration capacity f r, then:

[0042]

[0043] where I is the cumulative infiltration amount up to the calculation time, and K s,n is the saturated hydraulic conductivity parameter of the maximum infiltration capacity of the nth layer, and θ s,n is the saturated soil volume water content of the nth layer, and θ i,n is the initial soil volume water content of the nth layer, and G(θ i,n , θ s,n ) is the comprehensive capillary driving force at the soil wetting front, and K s,c is the composite saturated hydraulic conductivity;

[0044] Integrating gives:

[0045]

[0046] where t pp is the equivalent time required for infiltration to I t under the parameter conditions in the time step [t, t+Δt], and I t is the cumulative infiltration amount up to the calculation start time t.

[0047] Furthermore, in step S6, the dynamic simulation of the infiltration and redistribution processes in the multi-layer soil during the calculation period [t, t+Δt] includes:

[0048] The composite unsaturated hydraulic conductivity K c (ψ) of the wetting front in the nth layer of soil is calculated by the following formula:

[0049]

[0050] where K j (ψ) is the hydraulic conductivity of the jth layer when it has the same ψ value as the wetting front, and K n (ψ) is the hydraulic conductivity of the nth layer where the wetting front with the ψ value is located, and T j is the soil thickness of the jth layer, and Z is the depth where the wetting front is located;

[0051] Then the redistribution control equation of the soil volume water content θ at the wetting front is expressed as an ordinary differential equation through the cumulative infiltration depth:

[0052]

[0053] where θ o,n is the soil surface volume water content of the nth layer, and θ i,n is the initial soil volume water content of the nth layer, R is the rainfall amount during the calculation period, and K c(ψ) is the composite multi-layer hydraulic conductivity, p is an empirical parameter reflecting the shape characteristics of the wetting front, and K s,n is the saturated hydraulic conductivity parameter of the maximum infiltration capacity of the nth layer, and G(θ i,n , θ o,n ) is the comprehensive capillary driving force corresponding to the volumetric water content on the surface of the wetting front, and I t+△t is the cumulative infiltration amount up to the calculation start time t + △t;

[0054] At the same time, further correct the volumetric water content θ value of the wetting front to obtain the updated volumetric water content θ new of the wetting front:

[0055]

[0056] where θ s,n is the soil saturated volumetric water content of the nth layer where the wetting front is located, and Γ is the redistribution coefficient;

[0057] A new wetting front depth Z needs to be calculated according to the following formula through mass conservation new :

[0058]

[0059] where T j is the thickness of the jth layer of soil, I is the cumulative infiltration amount up to the calculation time, θ o,j is the volumetric water content on the surface of the jth layer of soil, and θ i,j is the initial volumetric water content of the jth layer of soil;

[0060] If the wetting front depth Z calculated above exceeds the lower boundary of the nth layer of soil, it is necessary to repeat the formula for calculating the new wetting front depth Z through mass conservation under the condition that the capillary head ψ of the soil volumetric water content above and below any layer interface penetrated by the wetting front is equal new until or the wetting front reaches the lowest layer boundary of the calculation area, and k is the serial number of the soil layer.

[0061] The present invention also provides a moisture dynamic prediction system based on a layered soil improved Green-Ampt redistribution model, including:

[0062] A data acquisition unit for acquiring the layered soil parameters and rainfall data of the research object;

[0063] A start ponding time calculation unit for rainfall segmentation processing, and determining the start ponding time at the top of the soil within the calculation time step based on the composite multi-layer soil parameterization method and precipitation data for each calculation period;

[0064] A dynamic simulation unit for the time t when water starts to accumulate at the top of the soil p is compared with the end time t+Δt of the calculation time step. If the end time of the calculation time step is not greater than the time when water starts to accumulate at the top of the soil, there is no water accumulation at the top of the soil, and the increment of soil volume water content is calculated according to the rainfall magnitude and the calculation time step Δt; otherwise, there is water accumulation at the top of the soil during the calculation period, and within the range [t, t p from the start time of the calculation to the time when water accumulation occurs, the increment of soil volume water content is calculated according to the rainfall magnitude and the calculation time step t p -t, and within the range [t p , t+Δt] after water accumulation occurs until the end of the calculation period, the surface runoff and the increment of soil volume water content are calculated according to the rainfall magnitude, the calculation time step t+Δt-t p and the soil volume water content; the wetting front position is updated based on the improved Green-Ampt redistribution model for layered soil, and the dynamic simulation of the infiltration and redistribution processes in multi-layer soil within the calculation period [t, t+Δt] is realized;

[0065] A statistical unit for judging whether the end time t+Δt is greater than the rainfall end time. If it is greater, the runoff and cumulative infiltration data of the rainfall data under the determined soil parameters are statistically calculated, and the calculation ends; otherwise, it returns to the start water accumulation time calculation unit.

[0066] The present invention provides a device, including a memory and a processor, wherein:

[0067] The memory is used to store a computer program that can run on the processor;

[0068] The processor is used to execute the steps of the moisture dynamic prediction method based on the improved Green-Ampt redistribution model for layered soil when running the computer program.

[0069] The present invention provides a storage medium, on which a computer program is stored, and when the computer program is executed by at least one processor, the steps of the moisture dynamic prediction method based on the improved Green-Ampt redistribution model for layered soil are realized.

[0070] Beneficial effects: Compared with the prior art, the remarkable technical effects of the present invention are as follows: An innovative parametric method and a dynamic adjustment mechanism are proposed to achieve continuous simulation of any number of soil layers and multiple rainfall pulse events, enabling accurate characterization of the dynamic evolution of the wetting front and the water redistribution process in layered soils; it uses a two-dimensional coordinate system of water content-depth to replace the traditional discrete node division, significantly reducing the number of parameters and improving the calculation efficiency. At the same time, through implicit equation iteration and mass conservation constraints, problems such as poor convergence, large computational overhead, and mass balance errors that are prone to occur in the numerical solution of the Richards equation are effectively solved; the refined characterization of soil spatial heterogeneity and the analysis of multi-layer water transmission mechanisms in the model not only improve the ability to depict surface-groundwater interactions in runoff generation and concentration simulations but also provide a more reliable numerical tool for flood risk assessment and water resource management, especially suitable for the complex scenario requirements of large-scale distributed hydrological models. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 is a flowchart of the method of the present invention;

[0072] Figure 2 is a schematic diagram of the water distribution at the soil layer interface;

[0073] Figure 3 is a comparison schematic diagram of the cumulative infiltration and runoff generated by the water dynamic prediction method based on the improved Green-Ampt redistribution model for layered soils and the Hydrus-1D model. Among them, (a) is a comparison schematic diagram of the solutions of the method of the present invention and the HYDRUS-1D model in the first experimental scenario, (b) is a comparison schematic diagram of the solutions of the method of the present invention and the HYDRUS-1D model in the second experimental scenario, and (c) is a comparison schematic diagram of the solutions of the method of the present invention and the HYDRUS-1D model in the third experimental scenario. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0074] The following further elaborates on the present invention in detail with reference to the drawings and specific embodiments. The following described specific embodiments further elaborate on the purpose, technical solutions, and beneficial effects of the present invention. It should be understood that the following is only the specific embodiments of the present invention and is not used to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

[0075] As Figure 1 shown, a water dynamic prediction method based on an improved Green-Ampt redistribution model for layered soils according to the present invention includes the following steps:

[0076] S1. Obtain the layered soil parameters and rainfall data of the research object;

[0077] First, perform layered soil structure division on the soil. Vertically divide the soil into N layers, and independently define the hydraulic parameters for each layer of soil (the subscript j represents the j-th layer), including the hydraulic conductivity K representing the j-th layer j , the saturated volumetric water content θ when the j-th layer of soil is completely saturated s,j , simulate the initial volumetric water content θ of the j-th layer at the initial moment i,j , the residual volumetric water content θ of the j-th layer of soil r,j , and the saturated hydraulic conductivity parameter K of the maximum infiltration capacity of the i-th layer s,i , the thickness T of the i-th layer of soil i and the rainfall rate R.

[0078] In this embodiment, three rainfall scenarios are set using idealized data to simulate three experiments respectively, in order to compare the performance of the moisture dynamic prediction method model method based on the improved Green-Ampt redistribution model for layered soil and the HYDRUS-1D numerical simulation method. The rainfall dataset lasts for 12 hours. A short rainstorm lasting for 25 minutes occurs 30 minutes after the start of the simulation. Subsequently, a long rainstorm lasting for 5 hours and 50 minutes starts when the simulation proceeds to 4 hours and 55 minutes. And the rainfall intensities of the two rainfall pulses are the same. All can completely infiltrate into the soil during the first shorter rainfall pulse, while exceeding the potential infiltration rate during this period in the second longer rainfall pulse, and can also cause the wetting front depth generated by the second rainfall pulse to exceed the wetting front depth created by the first rainfall pulse and then merge. In the embodiment, it is assumed that the evapotranspiration is zero, and the simulated rainfall intensities in Experiments 1, 2, and 3 are 32 mm / h, 7 mm / h, and 15 mm / h respectively.

[0079] In this study, the soil is vertically divided into 3 layers under three experimental conditions, and the hydraulic parameters are independently defined for each layer. Table 1 shows the detailed soil parameters of the three experiments:

[0080] Table 1 Soil parameters used under three experimental conditions

[0081]

[0082] The rainfall intensity in Experiment 1 was 32 mm / h. Under the initial capillary head condition of ψ = -10000 mm, it consisted of three types of soils, namely loam, clay loam, and silty clay loam with given parameters from top to bottom, with thicknesses of 10 cm, 30 cm, and 30 cm respectively; the rainfall intensity in Experiment 2 was 7 mm / h. Under the initial capillary head condition of ψ = -1000 mm, it consisted of three types of soils, namely silt, silty clay loam, and silty clay with given parameters from top to bottom, with thicknesses of 10 cm, 30 cm, and 30 cm respectively; the rainfall intensity in Experiment 3 was 15 mm / h. Under the initial capillary head condition of ψ = -1000 mm, it consisted of three types of soils, namely sandy clay loam, sandy clay, and silt loam with given parameters from top to bottom, with thicknesses of 10 cm, 30 cm, and 30 cm respectively.

[0083] S2. The rainfall duration was segmented into multiple calculation periods of length △t. For each calculation period [t, t + △t], based on the composite multi-layer soil parameterization method and rainfall data, the starting water accumulation time t at the top of the soil within the calculation time step was determined. p , where t is the starting moment of this calculation period, and △t is the calculation time step;

[0084] The calculation of the composite saturated hydraulic conductivity dynamically calculates the composite saturated hydraulic conductivity K according to the current soil layer number n where the wetting front is located and the thickness of each layer. s,c reflects the comprehensive water conduction ability of the multi-layer soil, K s,c The calculation formula is:

[0085]

[0086] where Z is the depth reached by the wetting front, K s,i is the saturated hydraulic conductivity parameter of the maximum infiltration capacity of the i-th layer, K s,n is the saturated hydraulic conductivity parameter of the maximum infiltration capacity of the n-th layer, T i is the thickness of the i-th layer of soil.

[0087] For example, when the wetting front penetrates the first layer (T1 = 10 cm), K s,c = 31.2 cm·h -1 ; if it penetrates to the second layer (T2 = 40 cm), then

[0088] In the present invention, the composite saturated hydraulic conductivity K used within each calculation time step △t s,c will change due to different values of the wetting front depth Z. Therefore, for each calculation time step △t, it is necessary to judge the relationship between the soil infiltration capacity f r , the composite saturated hydraulic conductivity K s,c and the rainfall rate R.

[0089] Calculate the comprehensive capillary driving force G(θ b , θ) with the formula:

[0090]

[0091] where θ represents the soil volumetric water content at the wetting front, and θ b represents the soil volumetric water content below the wetting front. The values of ψ and ψ b correspond to the capillary pressures at the soil volumetric water contents θ and θ b respectively. K(ψ) is the hydraulic conductivity of the capillary pressure ψ, and K s is the saturated hydraulic conductivity.

[0092] For a certain position Z where the wetting front is located in the nth layer, the soil infiltration capacity f r can be expressed as:

[0093]

[0094] where I is the cumulative infiltration amount up to the calculation time, θ s,n is the saturated soil volumetric water content of the nth layer, θ i,n is the initial soil volumetric water content of the nth layer, and G(θ i,n , θ s,n ) is the comprehensive capillary driving force corresponding to the saturated volumetric water content at the wetting front.

[0095] If the rainfall rate R is less than the composite multi-layer saturated hydraulic conductivity K s,c , there will be no ponding within the calculation time step. At this time, the ponding time t p is greater than the end time t + △t of the calculation time step;

[0096] If the rainfall rate R is greater than the composite multi-layer saturated hydraulic conductivity K s,c but less than the soil infiltration capacity f r , there may be ponding within the calculation time step [t, t + △t]. At this time, first, the corresponding ponding occurrence time t s,c needs to be calculated according to the following formula and the composite multi-layer saturated hydraulic conductivity K p,n at this time step:

[0097]

[0098] where I t is the cumulative infiltration amount up to the calculation start time t.

[0099] Finally, if the rainfall rate R is greater than the soil infiltration capacity f r , ponding will directly occur at the beginning stage within the calculation time step [t, t + △t].

[0100] S3. Compare the start time of ponding at the soil surface \(t\) p with the end time \(t +\Delta t\) of the computational time step. If the end time of the computational time step is not greater than the start time of ponding at the soil surface, go to S4; otherwise, go to S5.

[0101] S4. No ponding occurs at the soil surface. Calculate the increment of soil volumetric water content according to the rainfall intensity and the computational time step \(\Delta t\), and then go to S6.

[0102] Calculating the increment of soil volumetric water content according to the rainfall intensity and the computational time step \(\Delta t\) specifically includes:

[0103] If the start time of ponding at the soil surface \(t\) p is greater than the end time \(t+\Delta t\) of the computational time step, no ponding will occur within the computational time step. At this time, the increment of cumulative infiltration \(\Delta I\) is calculated by the following formula:

[0104] \(\Delta I = R\cdot\Delta t\)

[0105] S5. Ponding occurs at the soil surface during the computational period. Calculate the increment of soil volumetric water content according to the rainfall intensity and the computational time step \(t - t\) within the range \([t, t\) p from the start time of the calculation to the time of ponding occurrence, and calculate the surface runoff and the increment of soil volumetric water content according to the rainfall intensity, the computational time step \(t+\Delta t - t\) p and the soil volumetric water content within the range \([t\) p , \(t+\Delta t]\) after the ponding occurrence until the end of the computational period, and then go to S6. p Calculating the increment of soil volumetric water content according to the rainfall intensity, the computational time step \(t+\Delta t - t\)

[0106] and the soil volumetric water content specifically includes: p No ponding will occur within the range \([t, t\)

[0107] from the start time of the calculation to the time of ponding occurrence. At this time, the increment of cumulative infiltration \(\Delta I\) is calculated by the following formula: p

[0108] \(\Delta I = R\cdot(t\) p \(-t)\)

[0109] For the computational period \([t\) p , \(t+\Delta t]\), when the rainfall rate \(R\) is greater than the soil infiltration capacity \(f\) r , then:

[0110]

[0111] Integrating gives:

[0112]

[0113] where t pp is the equivalent time required for infiltration into I under the condition of the time step [t, t + Δt]. t

[0114] In the method for predicting moisture dynamics based on the improved Green-Ampt redistribution model for layered soil, the Newton-Raphson iteration method is used to solve the above implicit equation as follows: First, the implicit equation is reconstructed into a functional form f(I) with the cumulative infiltration amount I as the unknown:

[0115]

[0116] During iterative solution, for the initial guess value I t+△t of the cumulative infiltration amount at the set time t + Δt, subsequently, by calculating the function value f(I t+△t ) and its derivative f'(I t+△t ), according to It +△ t,k +1 = It +△ t,k - f(It +△ t,k) / f'(It +△ t,k) gradually approaches the true solution until the residual is less than the preset tolerance, where I t+△t,k+1 is the (k + 1)-th iteration value of the cumulative infiltration amount at time t + Δt, I t+△t,k is the k-th iteration value of the cumulative infiltration amount at time t + Δt, f(It +△ t,k) is the calculation result of the function f(I) when I = I t+△t,k , and f'(It +△ t,k) is the first derivative value of the function f(I) when I = I t+△t,k ; during this process, if the wetting front crosses the soil layer interface resulting in a parameter mutation, the composite multi-layer saturated hydraulic conductivity K s,c and the comprehensive capillary driving force G(θ i,n , θ s,n ) corresponding to the saturated volumetric water content at the wetting front need to be updated to ensure cross-layer physical continuity.

[0117] S6. Update the position of the wetting front based on the improved Green-Ampt redistribution model for layered soil to realize the dynamic simulation of the infiltration and redistribution processes in multi-layer soil during the calculation period [t, t + Δt];

[0118] The composite multi-layer unsaturated hydraulic conductivity K c (ψ) of the wetting front located in the n-th layer of soil is calculated by the following formula:

[0119] ​

[0120] Among them, K j (ψ) is the hydraulic conductivity of the j-th layer when it has the same ψ value as the wetting front, K n (ψ) is the hydraulic conductivity of the n-th layer where the wetting front with ψ value is located, T j is the soil thickness of the j-th layer, and Z is the depth where the wetting front is located.

[0121] Then, the redistribution control equation of the volumetric water content θ at the wetting front can be expressed as an ordinary differential equation through the cumulative infiltration depth:

[0122]

[0123] Among them, θ o,n is the volumetric water content at the surface of the n-th layer of soil, R is the rainfall during the calculation period, K c (ψ) is the composite multi-layer hydraulic conductivity, p is an empirical parameter reflecting the shape characteristics of the wetting front. Among them, for the value of the parameter p, considering the rectangular-shaped wetting front selected in the present invention, p = 1.0 is used when the rainfall data R > 0, and p = 1.7 is used when the rainfall data R = 0, G(θ i,n , θ o,n ) is the comprehensive capillary driving force corresponding to the volumetric water content at the surface of the wetting front, I t+△t is the cumulative infiltration amount up to the start time t + △t of the calculation.

[0124] The adaptive step-size Runge - Kutta algorithm efficiently solves the change of the volumetric water content θ at the wetting front by dynamically adjusting the time step; the process is as follows: when the model enters the redistribution stage of the volumetric water content at the wetting front, the change in soil water content is described by the ordinary differential equation dθdt; the fourth-order Runge - Kutta method is used for single-step prediction, and the step-size is dynamically adjusted by comparing the errors of two half-step integrations and a full-step integration. That is, when the relative error exceeds the allowable threshold, the step-size is reduced according to an exponential law combined with a safety factor of 0.7. After the error meets the standard, the step-size is expanded by a factor of 1.3. At the same time, an error correction factor of 1 / 15 is introduced to improve the accuracy; the algorithm sets a maximum number of iteration steps of 100 and a minimum value protection of 0.001 during the integration process, and achieves the optimal balance between calculation efficiency and numerical stability during the evolution of the wetting front position through iterative loops.

[0125] At the same time, the present invention uses an empirical formula to calculate the redistribution coefficient Γ to further correct the value of the volumetric water content θ at the wetting front:

[0126]

[0127] Among them, T R is the redistribution time; N Ris the number of redistribution times.

[0128] After the update obtained according to the above calculation, the volumetric water content θ of the wetting front new is valued as

[0129]

[0130] where θ s,n is the saturated volumetric water content of the soil in the nth layer where the wetting front is located, because the volumetric water content θ of the wetting front cannot exceed the saturated volumetric water content of the soil in the layer where it is located.

[0131] The new wetting front depth Z needs to be calculated according to the mass conservation by the following formula new :

[0132]

[0133] where θ o,j is the volumetric water content on the surface of the jth layer of soil; θ i,j is the initial volumetric water content of the jth layer of soil; T j is the thickness of the jth layer of soil; I is the cumulative infiltration amount up to the calculation time;

[0134] The specific process is as Figure 2 shown, which shows a schematic diagram of the volumetric water content profile across the soil layer interface. The green area represents the initial volumetric water content of the soil, the purple area represents the volumetric water content profile of the wetting front at the start of the time interval, the yellow area represents the change in the volumetric water content profile of the wetting front at the end of the time interval, the blue line represents the initial soil volumetric water content at different depths, which varies between different soil layers but has the same capillary head, and the red line represents the saturated volumetric water content of different soil layers.

[0135] If the wetting front depth Z obtained from the above calculation exceeds the lower boundary of the nth layer of soil, it is necessary to repeat the calculation of the new wetting front depth Z by mass conservation under the condition that the capillary heads ψ of the soil volumetric water content above and below any layer interface penetrated by the wetting front are equal new of the formula until or the wetting front reaches the lowest layer boundary of the calculation area, where k is the serial number of the soil layer;

[0136] S7. Judge whether the end time t + △t is greater than the rainfall end time. If it is greater, enter S8; otherwise, return to S2;

[0137] S8. Statistically analyze the runoff and cumulative infiltration amount data of the rainfall data under the determined soil parameters, and the calculation ends.

[0138] As shown in Table 2 and Figure 3As shown in (a) to (c), in the three experiments, the performance indicators of infiltration and runoff showed high consistency and reliability. First, the KGE values exceeded 0.96 in all experiments, indicating good performance of the model in simulating infiltration and runoff. In particular, the KGE value in Experiment 1 reached 0.981, suggesting that the method for predicting soil moisture dynamics based on the improved Green-Ampt redistribution model for layered soil was very accurate under this condition. The NSE values also performed excellently. The NSE in Experiment 1 was 0.995, in Experiment 2 was 0.969, and in Experiment 3 was 0.982, all close to 1, further verifying the accuracy of the method. In terms of bias, the PBIAS values varied in different experiments. The PBIAS values were within 1% bias in all experiments, indicating that the method for predicting soil moisture dynamics based on the improved Green-Ampt redistribution model for layered soil had a prediction accuracy that met expectations in this case. The root mean square error (RMSE) values were low in each experiment. The RMSE in Experiment 1 and Experiment 2 were 0.05 and 0.045, 0.029 and 0.029 respectively, and the RMSE in Experiment 3 was 0.048 and 0.044, suggesting that the prediction error of the model was small and had high reliability.

[0139] In terms of cumulative flux, the results of the method for predicting soil moisture dynamics based on the improved Green-Ampt redistribution model for layered soil and the HYDRUS-1D model showed similar trends. The results of the method for predicting soil moisture dynamics based on the improved Green-Ampt redistribution model for layered soil and HYDRUS-1D simulation in Experiment 1 are as Figure 3 (a) shown, where the cumulative infiltration value of the patented method was 64.845 mm and that of HYDRUS-1D was 65.055 mm, and the difference between them was only -0.210 mm, indicating that the prediction results of the two methods were very close under this condition. The simulation results of Experiment 2 are as Figure 3 (b) shown, the method for predicting soil moisture dynamics based on the improved Green-Ampt redistribution model for layered soil was 16.969 mm and HYDRUS-1D was 17.043 mm, with a difference of -0.74 mm. The simulation results of Experiment 3 are as Figure 3 (c) shown, for the method for predicting soil moisture dynamics based on the improved Green-Ampt redistribution model for layered soil, its cumulative infiltration value was 30.447 mm, while the value of HYDRUS-1D was 30.145 mm, showing a difference of 0.302 mm between the two, still within an acceptable range. Overall, these data demonstrate the soil moisture dynamics changes under different experimental conditions for the two methods, verifying the effectiveness and reliability of the method for predicting soil moisture dynamics based on the improved Green-Ampt redistribution model for layered soil in hydrological process simulation, and providing an important reference basis for the function and practicality of the method.

[0140] Table 2 Error metrics describing the differences between the patented method and HYDRUS-1D simulations

[0141]

[0142]

[0143] The moisture dynamics prediction system based on the improved Green-Ampt redistribution model for layered soils according to the present invention comprises:

[0144] A data acquisition unit for acquiring the layered soil parameters and rainfall data of the research object;

[0145] A starting ponding time calculation unit for performing rainfall segmentation processing, and determining the starting ponding time at the top of the soil within the calculation time step based on the composite multi-layer soil parameterization method and precipitation data for each calculation period;

[0146] A dynamic simulation unit for comparing the starting ponding time t p at the top of the soil with the end time t+Δt of the calculation time step. If the end time of the calculation time step is not greater than the starting ponding time at the top of the soil, there is no ponding at the top of the soil, and the increment of the soil volume water content is calculated according to the rainfall magnitude and the calculation time step Δt; otherwise, there is ponding at the top of the soil within the calculation period, and within the range [t,t p from the calculation start time to the ponding occurrence time, the increment of the soil volume water content is calculated according to the rainfall magnitude and the calculation time step t p -t. Within the range [t p ,t+Δt] from the ponding occurrence time to the end of the calculation period, the surface runoff and the increment of the soil volume water content are calculated according to the rainfall magnitude, the calculation time step t+Δt-t p and the soil volume water content; the wetting front position is updated based on the improved Green-Ampt redistribution model for layered soils to realize the dynamic simulation of the infiltration and redistribution processes in the multi-layer soil within the calculation period [t,t+Δt];

[0147] A statistics unit for determining whether the end time t+Δt is greater than the rainfall end time. If it is greater, the runoff and cumulative infiltration data of the rainfall data under the determined soil parameters are statistically calculated, and the calculation ends; otherwise, it returns to the starting ponding time calculation unit.

[0148] An apparatus of the present invention includes a memory and a processor, wherein:

[0149] The memory is used for storing a computer program that can run on the processor;

[0150] A processor, configured to execute the steps of the above-described method for predicting moisture dynamics based on the improved Green-Ampt redistribution model for layered soil when running the computer program.

[0151] A storage medium of the present invention, on which a computer program is stored. When the computer program is executed by at least one processor, it implements the steps of the above-described method for predicting moisture dynamics based on the improved Green-Ampt redistribution model for layered soil, and can achieve the same technical effects as the above method.

Claims

1. A method for predicting moisture dynamics based on an improved Green-Ampt redistribution model for layered soil, characterized in that, It includes the following steps: S1. Obtain the layered soil parameters and rainfall data of the research object; S2. Process the rainfall in segments. For each calculation period, based on the composite multi-layer soil parameterization method and precipitation data, determine the time when water starts to accumulate at the top of the soil within the calculation time step; S3. Compare the time \(t\) when water starts to accumulate at the top of the soil p with the end time \(t +\Delta t\) of the calculation time step. If the end time of the calculation time step is not greater than the time when water starts to accumulate at the top of the soil, proceed to S4; otherwise, proceed to S5. S4. When there is no water accumulation at the top of the soil, calculate the increment of soil volumetric water content according to the rainfall magnitude and the calculation time step △t, and enter S6; S5. Waterlogging occurs at the top of the soil during the calculation period. In the range [t, t p from the start time of the calculation to the time when waterlogging appears, calculate the increment of soil volumetric water content according to the rainfall magnitude and the calculation time step t p -t. After waterlogging appears, in the range [t p , t + △t] from the time when waterlogging appears to the end of the calculation period, calculate the surface runoff and the increment of soil volumetric water content according to the rainfall magnitude, the calculation time step t + △t - t p and the soil volumetric water content, and enter S6; S6. Update the wetting front position based on the improved Green-Ampt redistribution model for layered soil, and realize the dynamic simulation of the infiltration and redistribution processes in the multi-layer soil within the calculation period [t, t+△t]; S7. Judge whether the end time t+△t is greater than the rainfall end time. If it is greater, enter S8; otherwise, return to S2; S8. Statistically calculate the runoff and cumulative infiltration data of the rainfall data under the determined soil parameters, and the calculation ends.

2. The method for predicting moisture dynamics based on the improved Green-Ampt redistribution model for layered soil according to claim 1, characterized in that, In step S1, the soil is vertically divided into N layers, and the hydraulic parameters are independently defined for each layer, including hydraulic conductivity, saturated hydraulic conductivity parameter with the maximum infiltration capacity, saturated soil volumetric water content when the soil is completely saturated, initial soil volumetric water content at the initial simulation moment, residual soil volumetric water content of the soil, soil thickness, and rainfall rate.

3. The moisture dynamics prediction method based on the improved Green-Ampt redistribution model for layered soil according to claim 1, characterized in that The composite multi-layer soil parameterization method in step S2 is: Dynamically calculate the composite saturated hydraulic conductivity K according to the current soil layer number n where the wetting front is located, as well as the thickness of each layer s,c K reflects the comprehensive water conductivity of multi-layer soil s,c The calculation formula is as follows: where Z is the depth reached by the wetting front, K s,i is the saturated hydraulic conductivity parameter of the maximum infiltration capacity of the i-th layer, K s,n is the saturated hydraulic conductivity parameter of the maximum infiltration capacity of the n-th layer, T i is the soil thickness of the i-th layer.

4. The moisture dynamic prediction method based on the improved Green-Ampt redistribution model for layered soil according to claim 1, wherein, Determine the time t when water starts to accumulate at the top of the soil within the calculated time step in step S2 p , including: First, calculate the comprehensive capillary driving force G(θ b , θ) with the formula as follows: where θ represents the soil volumetric water content at the wetting front, and θ b represents the soil volumetric water content below the wetting front, and ψ and ψ b correspond to the capillary pressures at the soil volumetric water contents θ and θ b respectively, K(ψ) is the hydraulic conductivity of the capillary pressure ψ, and K s is the saturated hydraulic conductivity; When the wetting front is at a certain position Z in the nth layer, the soil infiltration capacity f r is expressed as: Among them, K s,n is the saturated hydraulic conductivity parameter of the maximum infiltration capacity of the nth layer, θ s,n is the saturated soil volumetric water content of the nth layer, θ i,n is the initial soil volumetric water content of the nth layer, G(θ i,n , θ s,n ) is the comprehensive capillary driving force corresponding to the saturated volumetric water content at the wetting front, I is the cumulative infiltration amount up to the calculation time, K s,c is the composite saturated hydraulic conductivity; If the rainfall rate R is less than the composite saturated hydraulic conductivity K s,c , there will be no ponding within the calculated time step, and in this case, the ponding time t p is greater than the end time t + △t of the calculated time step; If the rainfall rate R is greater than the composite saturated hydraulic conductivity K s,c but less than the soil infiltration capacity f r , then ponding may occur within the calculation time step [t, t+△t]. At this time, it is first necessary to calculate the corresponding ponding occurrence time t according to the following formula and the composite saturated hydraulic conductivity K at this time step s,c : p,n : Among them, I t is the cumulative infiltration volume up to the start time t of the calculation; Finally, if the rainfall rate R is greater than the soil infiltration capacity f r , then ponding will directly occur at the beginning stage within the time step [t, t + △t].

5. The moisture dynamic prediction method based on the improved Green-Ampt redistribution model for layered soil according to claim 1, wherein, In step S4, calculating the increment of soil volumetric water content according to the rainfall magnitude and the calculation time step △t includes: If the time t when water starts to accumulate at the top of the soil p is greater than the end time t + △t of the calculated time step, no water accumulation will occur within the calculated time step. In this case, the increment △I of the cumulative infiltration amount is calculated by the following formula: △I = R·△t where R is the rainfall rate.

6. The moisture dynamic prediction method based on the improved Green-Ampt redistribution model for layered soil according to claim 1, wherein, In step S5, calculate the increment of soil volume water content based on rainfall magnitude, calculation time step t+△t-t p and soil volume water content, including: There will be no ponding within the range [t, t p from the start time of calculation to the ponding occurrence time. At this time, the increment △I of the cumulative infiltration amount is calculated by the following formula: △I = R·(t p - t) where R is the rainfall rate; For the calculation period [t p , t + △t], when the rainfall rate R is greater than the soil infiltration capacity f r , then: where I is the cumulative infiltration amount up to the calculation time, K s,n is the saturated hydraulic conductivity parameter of the maximum infiltration capacity of the nth layer, θ s,n is the saturated soil volumetric water content of the nth layer, θ i,n is the initial soil volumetric water content of the nth layer, G(θ i,n , θ s,n ) is the comprehensive capillary driving force at the wetting front of the soil, K s,c is the composite saturated hydraulic conductivity; Integrating to obtain: Among them, t pp is the equivalent time required for infiltration into I under the parameter conditions at the time step [t, t + △t], t where I t is the cumulative infiltration amount up to the start time t of the calculation.

7. The method for predicting water dynamics based on the improved Green-Ampt redistribution model for layered soil according to claim 1, characterized in that In step S6, realizing the dynamic simulation of the infiltration and redistribution processes in the multi-layer soil within the calculation period [t, t+△t] includes: The composite unsaturated hydraulic conductivity K c (ψ) at the wetting front in the n-th layer of soil is calculated by the following formula: where K j (ψ) is the hydraulic conductivity of the j-th layer when it has the same ψ value as the wetting front, K n (ψ) is the hydraulic conductivity of the n-th layer where the wetting front with ψ value is located, T j is the soil thickness of the j-th layer, and Z is the depth where the wetting front is located; Then the redistribution control equation of the wetting front soil volumetric water content θ is expressed as an ordinary differential equation through the cumulative infiltration depth: where θ o,n is the volumetric water content at the surface of the n-th soil layer, θ i,n is the initial volumetric water content of the n-th soil layer, R is the rainfall during the calculation period, K c (ψ) is the composite multi-layer hydraulic conductivity, p is an empirical parameter reflecting the shape characteristics of the wetting front, K s,n is the saturated hydraulic conductivity parameter of the maximum infiltration capacity of the n-th layer, G(θ i,n , θ o,n ) is the comprehensive capillary driving force corresponding to the volumetric water content at the surface of the wetting front, I t+△t is the cumulative infiltration amount up to the calculation start time t + △t; Meanwhile, further correct the volumetric water content θ value at the wetting front to obtain the updated volumetric water content θ at the wetting front new The value is: where θ s,n is the saturated volumetric water content of the nth layer where the wetting front is located, and Γ is the redistribution coefficient; The new wetting front depth Z needs to be calculated based on the mass conservation according to the following formula new :[[]]END]] Among them, T j is the thickness of the j-th layer of soil, I is the cumulative infiltration amount up to the calculation time, and θ o,j is the volumetric water content at the surface of the j-th layer of soil, and θ i,j is the initial volumetric water content of the j-th layer of soil; If the wetting front depth Z calculated above exceeds the lower boundary of the nth soil layer, it is necessary to repeat the calculation of the new wetting front depth Z through mass conservation under the condition that the capillary heads ψ of the soil volumetric water contents above and below any layer interface penetrated by the wetting front are equal new the formula of until either the wetting front reaches the lowest layer boundary of the calculation area, where k is the serial number of the soil layer 8. A moisture dynamic prediction system based on the improved Green-Ampt redistribution model for layered soil, characterized in that, It includes: A data acquisition unit for obtaining the layered soil parameters and rainfall data of the research object; A start water accumulation time calculation unit for processing the rainfall in segments. For each calculation period, based on the composite multi-layer soil parameterization method and precipitation data, determine the time when water starts to accumulate at the top of the soil within the calculation time step; A dynamic simulation unit is used to calculate the time t when water starts to accumulate at the top of the soil p and compare it with the end time t+△t of the calculation time step. If the end time of the calculation time step is not greater than the time when water starts to accumulate at the top of the soil, there is no water accumulation at the top of the soil, and the increment of soil volume water content is calculated according to the rainfall magnitude and the calculation time step △t; otherwise, there is water accumulation at the top of the soil during the calculation period. In the range [t, t p from the start time of the calculation to the time when water accumulation occurs, the increment of soil volume water content is calculated according to the rainfall magnitude and the calculation time step t p -t. In the range [t p , t+△t] from the time when water accumulation occurs to the end of the calculation period, the surface runoff and the increment of soil volume water content are calculated according to the rainfall magnitude, the calculation time step t+△t-t p and the soil volume water content; Update the wetting front position based on the improved Green-Ampt redistribution model for layered soil, and realize the dynamic simulation of the infiltration and redistribution processes in the multi-layer soil within the calculation period [t, t+△t]; A statistical unit for judging whether the end time t+△t is greater than the rainfall end time. If it is greater, statistically calculate the runoff and cumulative infiltration data of the rainfall data under the determined soil parameters, and the calculation ends; otherwise, return to the start water accumulation time calculation unit.

9. A device, characterized in that, It includes a memory and a processor, where: The memory is used to store a computer program that can run on the processor; The processor is used to execute the steps of the moisture dynamic prediction method based on the improved Green-Ampt redistribution model for layered soil as described in any one of claims 1-7 when running the computer program.

10. A storage medium, characterized in that, The computer program is stored on the storage medium, and when the computer program is executed by at least one processor, it realizes the steps of the moisture dynamic prediction method based on the improved Green-Ampt redistribution model for layered soil as described in any one of claims 1-7.

Images