A robust calculation method of soil hydrological parameters based on one-dimensional upward infiltration test

Abstract

The application discloses a kind of soil hydraulics parameter robust calculation method based on one-dimensional upward infiltration test.The application constructs piston type water profile model (PASCP), and the innovative method is converted into equivalent homogeneous region by the inhomogeneous water content distribution of unsaturated zone, effectively expand the applicable scope of model.Based on the model, the application establishes analytical infiltration model considering saturation time, and introduces interval constraint optimization objective function, significantly improves the inversion accuracy and robustness of key hydraulic parameters of Brooks-Corey model.Experiments prove that the method shows excellent applicability under different soil texture conditions, especially under the condition of measurement error, it shows higher result reliability.The parameter estimation method based on standard cutting ring provided by the application provides an efficient and practical technical means for the rapid and accurate characterization of soil hydraulic properties, and has wide application prospect in the fields of agricultural irrigation, hydrological simulation, etc.

Description

Technical Field

[0001] This invention relates to a robust calculation method for soil hydraulic parameters based on a one-dimensional upward infiltration test, belonging to the fields of soil physics and soil hydrology. Background Technology

[0002] Soil hydraulic properties, including the soil moisture characteristic curve θ(h) and the unsaturated hydraulic conductivity curve K(θ), are key input parameters for hydrological models based on the Richards equations. Together, they determine the transport characteristics of soil moisture in the unsaturated zone. Accurate measurement of these soil hydraulic parameters is crucial not only for field-scale studies of water and solute transport but also for global-scale hydrological and energy cycle simulations. In the laboratory, traditional methods for measuring θ(h) and K(θ) are the pressure film method and the transient profile method, respectively. While theoretically simple, these methods are cumbersome and time-consuming in practice. To improve efficiency, numerical inversion techniques based on transient flow experiments to simultaneously estimate soil θ(h) and K(θ) have rapidly developed, including single / multi-step drainage experiments, evaporation methods, and upward infiltration methods. Among these, upward infiltration inversion techniques based on ring sample soil samples have achieved significant breakthroughs in the last decade. This method utilizes capillary-driven anti-gravity infiltration, resulting in slower water movement. It is particularly suitable for accurately measuring water flow through a 5cm high ring cutter, thus providing an effective approach for large-scale soil hydraulic surveys. However, existing numerical inversion techniques still face a fundamental challenge: the lack of a universal and efficient solution algorithm for the highly nonlinear Richards equations. Current numerical schemes often require customized spatial discretization strategies based on different soil types and boundary conditions. This strong scenario dependence severely restricts their wider practical application.

[0003] To address the limitations of generality and complexity in the aforementioned numerical methods, researchers have gradually shifted their focus to analytical and semi-analytical physical models in order to more efficiently predict soil hydraulic parameters. Physical models have also been applied to parameter inversion under upward infiltration conditions. Moret-Fernandez and Latorre (2016), building upon the Haverkamp model, combined it with the saturated hydraulic conductivity K determined based on Darcy's law. s The permeability S and shape parameter β were successfully derived, and θ(h) was further estimated. Chinese Patent Publication No. CN114397427A discloses a method for predicting soil hydraulic properties based on the infiltration process of ring-sampling soil samples, which provides a new algorithm that can simultaneously estimate θ(h) and permeability S. However, both of the above methods require additional determination of the soil's saturated hydraulic conductivity K. sAs an input parameter, θ(h) is estimated. And as van Genuchten (1991), the proposer of the classical soil hydraulics model (van Genuchten model), pointed out, K... s Parameters are often difficult to measure accurately. To address this issue, Wu et al., in a 2022 paper published in the *Journal of Hydrology*, proposed a constant pressure upseeking physical model (ASCP) based on patent CN114397427A. This model can simultaneously estimate key parameters describing θ(h) and K(θ), including the pore distribution index n and the reciprocal of the inlet suction force 1 / h. d and saturated hydraulic conductivity K s However, accurate estimation of parameters using the ASCP method relies on strict equality constraints, namely, the need to accurately determine the cumulative infiltration at a single point when the wetting front first reaches the soil surface (I0). * ) and the corresponding height of the soil column (L) * However, in actual laboratory measurements, due to the limited accuracy of barometric pressure sensors and the interference of air bubbles in the Marble flask on the indoor pressure (Latorre et al., 2015), the pressure at which the wetting front just reaches the soil surface is significantly reduced. * Precise measurement is often difficult, leading to the inability to meet strict equality constraints and thus weakening the performance of the ASCP method in real-world scenarios. Therefore, how to reduce the impact of measurement error on I... * The impact of wetting fronts and the need to enhance the robustness of methods have become key issues that urgently need to be addressed. Furthermore, experimental evidence from Mort-Fernandez and Latorre (2019) indicates that soil water absorption continues slowly after the wetting front reaches the soil surface until a final steady state is reached. Considering the significant differences in the time required for different soil textures to reach a steady state, infiltration data at this stage may contain more crucial information about soil hydraulic properties. Integrating this information into the ASCP method could further improve the accuracy of hydraulic parameter estimation. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of the prior art and provide a robust calculation method for soil hydraulic parameters based on a one-dimensional upward infiltration test, thereby further improving the accuracy of hydraulic parameter estimation.

[0005] The technical solution adopted in this invention is as follows:

[0006] A robust method for calculating soil hydraulic parameters based on a one-dimensional upward infiltration test includes the following steps:

[0007] Step A: Apply a constant lower boundary water head h to the ring sampler soil sample. pOne-dimensional upward infiltration tests were conducted on homogeneous soil. The cumulative infiltration I of the ring sample was recorded as a function of time t, the time t* when the wetting front first reached the surface of the ring sample, and the cumulative infiltration I when the infiltration curve tended to stabilize. end ;

[0008] Step B: Based on the initial volumetric water content θ of the ring sample soil sample i The height of the ring cutter L* and the cumulative infiltration rate I when the soil reaches full saturation. end The saturated volumetric water content θ of the soil was calculated. s The calculation formula is:

[0009]

[0010] Step C: Based on the relationship between the cumulative infiltration amount I and time t obtained in Step A, determine the single-point cumulative infiltration amount I* corresponding to the moment t* when the wetting front first reaches the surface of the ring sample soil sample;

[0011] Step D: Combine the saturated volumetric water content θ of the ring sample. s Initial volumetric water content θ i Residual volumetric water content θ r The cumulative infiltration rate I* at a single point when the wetting front first reaches the soil surface, and the observed data I of the cumulative infiltration rate changing over time. obs (t i The parameters n and 1 / h to be estimated are obtained by using a parameter optimization inversion method based on interval constraints. d and K s .

[0012] Preferably, in step D, the objective function of the interval-constrained parameter optimization inversion method is constructed according to formula (14):

[0013]

[0014] Preferably, the first term of formula (14), I obs (t i The measured value obtained in step A represents the observation time point t. i The cumulative infiltration amount I at time t, I sim (t i ;n,1 / h d ,K s ) represents the observation time point t i The simulated cumulative infiltration amount is calculated using the following formula:

[0015] ① Time to soil saturation (t) s The function formula for ) is:

[0016]

[0017] In the formula z s * It is the equivalent wetting front z fe The function expression, when the equivalent wetting front z fe Reaching the soil surface (i.e., z) fe =L * When z reaches saturation, the soil is saturated. s * This can be explicitly stated as:

[0018]

[0019] ② Regarding the infiltration time t <t s At that time, the cumulative soil infiltration I is a function of time t as follows:

[0020]

[0021] ③ When t≥t s When the soil reaches full saturation, the cumulative infiltration remains constant, and the calculation formula is as follows:

[0022] I=(θ s -θ i )L * (10);

[0023] The parameters to be estimated, n and 1 / h d and K s Substituting this into the equation, we can obtain the measured value I. obs (t i ) and simulated value I sim (t i ;n,1 / h d ,K s The root mean square error between them tends to be minimized;

[0024] The second term of formula (14), I sim (L * ;n,1 / h d ,K s The cumulative infiltration amount is calculated using the following formula:

[0025]

[0026] Where z f =L*, where σ is the penalty factor [–], and ∈ is the constraint threshold [mm], the value of which depends on the upper limit of the maximum measurement error in actual observation.

[0027] The remaining parameters are expressed according to the following relationship:

[0028]

[0029] Preferably, formula (14) is solved using a genetic algorithm (GA package) in R language: the algorithm generates 1000 sets of parameters to be estimated n and 1 / h in each generation. d and K s The algorithm stops running when the value of φ in formula (14) no longer increases within 500 consecutive iterations, or when the cumulative number of iterations reaches 5000; a local search strategy is enabled; the parameter search range is K. s =[10 -6 10 -2 ]cm·s -1 1 / h d =[0.01,1]cm -1 , n = [0.1, 0.7], for K s With 1 / h d Perform logarithmic transformation; set the constraint threshold ε to 0.4 mm; set σ to... Where .Machine$double.xmax represents the maximum double-precision floating-point number that can be represented in R.

[0030] This invention proposes a robust calculation method for soil hydraulic parameters based on a one-dimensional upward infiltration experiment (PASCP method). Unlike the ASCP method, which can only describe the infiltration process immediately before the wetting front reaches the soil surface, PASCP further expands its applicability by incorporating soil hydraulic characteristics from the steady-state stage (the infiltration stage that continues after the wetting front reaches the soil surface) into the parameter estimation. This achieves a complete simulation of the entire upward infiltration process from the initial wetting state to full saturation, significantly improving the accuracy of hydraulic parameter estimation. Furthermore, to reduce the impact of measurement errors on the single-point cumulative infiltration (Ig) at the first arrival of the wetting front at the soil surface... * To mitigate the impact of measured values, a parameter optimization method based on interval constraints is proposed using the PASCP method to robustly estimate soil hydraulic parameters (n, 1 / h) of the Brooks-Corey model simultaneously. d and K s ).

[0031] The beneficial effects of this invention are as follows:

[0032] 1. The PASCP method effectively solves the problem of limited applicability of the ASCP method after the wetting front reaches the soil surface, and provides an important new approach for robust estimation of soil hydraulic parameters;

[0033] 2. By introducing a parameter optimization inversion method with interval constraints, the impact of measurement errors on the single-point cumulative infiltration (Ig) at the first arrival of the wetting front at the soil surface is reduced. * The influence of the measured values ​​significantly improved the robustness of soil hydraulic parameter estimation;

[0034] 3. The parameter estimation method based on standard ring sampler provides an effective technical means for rapidly and reliably characterizing soil hydraulic properties. Attached Figure Description

[0035] Figure 1 This is a schematic diagram of the evolution of the water profile during soil infiltration. (a) Profile shape when the wetting front first reaches the soil surface; (b) Profile changes as the infiltration process continues after the wetting front reaches the surface; (c) Profile shape when the soil reaches a fully saturated steady state.

[0036] Figure 2 This is a flowchart of the calculation process of the present invention.

[0037] Figure 3 This is an example of a parameter inversion experiment for clay soil samples. (a) Schematic diagram of the one-dimensional infiltration experimental setup and measurement process for laboratory ring sample soil; (b) Measured cumulative infiltration curve and key time points, including the moment when the wetting front first reaches the soil surface (t). * ) and the time when the soil reaches complete saturation (t) s (c) Iteration trajectory of PASCP and ASCP methods in the parameter estimation process; (d) Comparison of soil moisture characteristic curves obtained by PASCP and ASCP methods with measured values.

[0038] Figure 4 This is a parameter estimation diagram of five soils with significant differences in texture, determined in the laboratory.

[0039] Figure 5 This is a parameter estimation map based on numerical simulation of seven soils with significant differences in texture. (a) The PASCP method includes infiltration information during the steady-state stage; (b) The PASCP method excludes infiltration information during the steady-state stage.

[0040] The present invention will be further described below with reference to specific embodiments. Detailed Implementation

[0041] The following description, combined with soil infiltration tests, further illustrates the content of this invention, but should not be construed as limiting the invention. Modifications and substitutions made to the methods, steps, or conditions of this invention without departing from its spirit and essence are all within the scope of this invention. Unless otherwise specified, the technical means used in the following embodiments are conventional means well known to those skilled in the art.

[0042] Example 1: Derivation of the PASCP Method

[0043] 1. Piston-type moisture profile construction method

[0044] The soil hydraulic parameter inversion method based on the piston-type moisture profile model (PASCP) of this invention is, specifically, an improvement and extension of the ASCP method previously proposed by the applicant (Wu et al., 2022, Journal of Hydrology). In the ASCP method previously proposed by the applicant, the soil moisture profile consists of saturated and unsaturated regions, where the moisture distribution is uniform in the saturated region and exhibits a significant gradient change in the unsaturated region. Figure 1 a) Its cross-sectional mathematical expression is:

[0045]

[0046] Where z represents soil depth [cm], with its zero point located at the bottom of the soil and the coordinate axis pointing upwards; z f The distance the moist front advances [cm]; z s S represents the length of the saturation region [cm]. ei To correspond to the initial water content θ i The effective saturation [–]. n is the porosity distribution index. Parameters a and b are defined as follows:

[0047]

[0048] However, the ASCP method can only describe the infiltration process just before the wetting front reaches the soil surface. Once the wetting front reaches the surface, the unsaturated area gradually becomes saturated. Figure 1 (b) The water profile deviates from the description in the ASCP formula (Equation 1). To address this limitation, this invention proposes a novel piston-type water profile construction method (Modified ASCP with piston-type water profile, or PASCP for short) based on Green and Ampt's (1911) piston flow assumption. Its innovation lies in avoiding explicit simulation of the gradual saturation process in the unsaturated region while accurately reflecting the cumulative infiltration process at this stage, effectively expanding the applicability of the ASCP method. Specifically, the heterogeneous water content distribution in the unsaturated region is transformed into an equivalent homogeneous region using an integral method, and the equivalent length z of the unsaturated region is defined. ufe The calculation formula is as follows:

[0049]

[0050] Among them, z ufe z represents the equivalent length of the unsaturated region. s Let z be the length of the saturation region. s With z ufe By adding them together, the equivalent piston-type wetting front depth z can be constructed. fe Its expression is:

[0051] zfe =z s +z ufe (5)

[0052] Using the aforementioned piston-shaped moisture profile model, the PASCP method can effectively simulate the upward movement of a wetting front in a piston-shaped moisture profile. When the equivalent wetting front z... fe When it reaches the soil surface (e.g.) Figure 1 As shown in c), the soil column enters a fully saturated state and reaches a steady state; this moment is defined as the soil saturation time (t). s [min].

[0053] 2. Analytical Model Derivation

[0054] Using the aforementioned piston-shaped water profile, a one-dimensional analytical model for water infiltration applicable to constant head boundary conditions is further derived based on the principle of mass conservation.

[0055] ④ Time to soil saturation (t) s The function formula for ) is:

[0056]

[0057] In the formula z s * It is the equivalent wetting front z fe The function expression. When the equivalent wetting front z fe Reaching the soil surface (i.e., z) fe =L * When z reaches saturation, the soil is saturated. s * This can be explicitly stated as:

[0058]

[0059] It is worth noting that t s Depends on parameters n and 1 / h d and K s This demonstrates its application value in parameter inversion.

[0060] ⑤ Regarding the infiltration time t <t s At that time, the cumulative soil infiltration I is a function of time t as follows:

[0061]

[0062] ⑥ When t≥t s When the soil reaches full saturation, the cumulative infiltration remains constant, and the calculation formula is as follows:

[0063] I=(θ s -θ i )L* (10)

[0064] ⑦ Cumulative infiltration rate I and the advancing distance of the wetting front z f The functional relationship is as follows:

[0065]

[0066] The unknown parameter in the above formula is defined as follows: θ i The initial volumetric water content of the soil [cm] -3 cm -3 ], θ s Soil saturated volumetric water content [cm] -3 cm -3 ], h p Water head [cm] for the lower boundary of the soil, L * The soil height is [cm]. The other parameters are expressed according to the following relationship:

[0067]

[0068] Soil residual volumetric water content θ r The determination method is as follows: The sieved soil is stored in a dry air environment with a relative humidity of less than 15% for several months until its mass remains constant. The water content of the soil at this point is the residual mass water content, which is then multiplied by the bulk density of the ring sample to obtain the residual volumetric water content. The initial volumetric water content θ is also measured. i The soil saturated volumetric water content θ was determined using the oven-drying method. s The determination method is as follows: it is obtained by inverting formula (10), and its expression is:

[0069]

[0070] Among them, I end The value [cm] represents the stable value when the cumulative infiltration curve reaches a plateau, and represents the total amount of water absorbed by the soil sample.

[0071] 3. Inversion method for parameter optimization with interval constraints

[0072] This invention further proposes an interval-constrained optimization inversion method based on the PASCP method to improve soil hydraulic parameters (n, h) d and K s The accuracy and robustness of the estimation. The objective function expression of the specific algorithm is:

[0073]

[0074] The first item represents the measured data I. obs (t i ) and simulation data Isim (t i ;n,1 / h d ,K s The root mean square error (RMSE) between the simulated values ​​I and I. sim (t i ;n,1 / h d ,K s ) is calculated using formulas (6)-(10), where t i For the observation time point, I sim This represents the cumulative soil infiltration at the corresponding time point. The second term indicates the interval constraint introduced through the penalty function to ensure the simulated cumulative infiltration I. sim (L * ;n,1 / h d ,K s The cumulative infiltration at a single point when the wetting front first reaches the soil surface (I) * The deviation between them is maintained within ±ε to account for the impact of experimental measurement errors on the inversion accuracy. Where I sim (L * ;n,1 / h d ,K s The result is obtained by formula (11), i.e., z. f =L*); σ is the penalty factor [–], ∈ is the constraint threshold [mm], the value of which depends on the upper limit of the maximum measurement error in actual observation.

[0075] The above measured data (cumulative infiltration I) i With time t i The relationship of change can be obtained using a portable rapid soil hydraulic property measuring device disclosed by the applicant (patent publication number: CN216209116U), with a constant water head h applied to the bottom of the ring sampler. p The measurement frequency was once per second. A camera was used to determine the precise moment (t*) when the wetting front reached the soil surface. From the measurement data (I~t), based on t... * The cumulative infiltration rate I when the wetting front reaches the soil surface is obtained. * Record the cumulative infiltration I when the infiltration curve tends to stabilize. end Measure the height L of the ring cutter. * .

[0076] Formula (14) is solved using a genetic algorithm (GA package) in R language: the algorithm generates 1000 sets of parameters to be estimated n and 1 / h in each generation. d and K s The algorithm stops running when the value of φ in formula (14) no longer increases within 500 consecutive iterations, or when the cumulative number of iterations reaches 5000; a local search strategy is enabled; the parameter search range is K. s =[10-6 10 -2 ]cm·s -1 1 / h d =[0.01,1]cm -1 , n = [0.1, 0.7], for K s With 1 / h d Perform a logarithmic transformation. Set the constraint threshold ε to 0.4 mm. Set σ to... Where .Machine$double.xmax represents the maximum double-precision floating-point number that can be represented in R.

[0077] Example 2

[0078] This embodiment conducts a parameter inversion experiment based on the PASCP method. Figure 2 The specific steps are as follows:

[0079] Step A: Infiltration Test Measurement Stage

[0080] (1) Apply a constant water head h to the bottom of the ring sampler soil sample. p ;

[0081] (2) Record the relationship between the cumulative infiltration volume I of the soil sample and time t;

[0082] (3) Use a camera to record the time t when the wetting front first reaches the soil sample surface. * .

[0083] Step B: Experimental Parameter Acquisition Stage

[0084] (1) From the measurement data, according to t * The cumulative infiltration rate I when the wetting front reaches the soil surface is obtained. * ;

[0085] (2) Determine the initial moisture content θ of the soil sample i and residual moisture content θ r ;

[0086] (3) Record the cumulative infiltration amount I when the infiltration curve tends to stabilize. end Calculate the saturated water content θ s ;

[0087] (4) Measure the ring cutter height L * .

[0088] Step C: Model parameter inversion stage:

[0089] (1) Construct an optimization objective function based on the genetic algorithm (GA) in R language, including RMSE term and interval constraint term;

[0090] (2) Set the population size to 1000, the maximum number of iterations to 5000 (or 500 generations), and enable the local search strategy;

[0091] (3) The parameter search range is K s =[10 -6 10 -2 ]cm·s -1 1 / h d =[0.01,1]cm -1 , n = [0.1, 0.7], for K s With 1 / h d Perform a logarithmic transformation;

[0092] (4) Soil hydraulic parameters n and 1 / h were obtained through optimization calculations. d and K s The best estimate.

[0093] The measured data in step A (the relationship between cumulative infiltration I and cumulative time t) obs (t i The results were obtained using a portable rapid soil hydraulic properties measuring device disclosed by the applicant (patent publication number: CN216209116U), with a measurement frequency of once per second.

[0094] The residual soil volumetric water content θ in step B r The determination method is as follows: The sieved soil is stored in a dry air environment with a relative humidity of less than 15% for several months until its mass remains constant. The water content of the soil at this point is the residual mass water content, which is then multiplied by the bulk density of the ring sample to obtain the residual volumetric water content. The initial volumetric water content θ is also measured. i The soil saturated volumetric water content θ was determined using the oven-drying method. s The determination was performed using formula (13).

[0095] The objective function in step C is shown in formula (14), and is solved using a genetic algorithm (GA) based on R language. The population size is set to 1000, and the algorithm is stopped after a maximum of 5000 iterations or 500 generations, provided that the fitness is not improved. A local search strategy is enabled, and the parameter search range is K. s =[10 -6 10 -2 ]cm·s -1 1 / h d =[0.01,1]cm -1 , n = [0.1, 0.7], for K s With 1 / h d Perform a logarithmic transformation. Set the constraint threshold ε to 0.4 mm. Set σ to... Where .Machine$double.xmax represents the largest double-precision floating-point number that can be represented in R. The simulated infiltration rate I in formula (14) sim (L * ;n,1 / h d ,K s ) Calculated using formula (11); Simulated data I sim (t i ;n,1 / h d ,K s The formulas (6)-(10) are used for calculation.

[0096] Example 3

[0097] This embodiment aims to verify the ability of the PASCP method to handle measurement errors under actual measurement conditions, and to demonstrate the improvement in parameter estimation robustness compared to the ASCP method. Following the steps of Embodiment 2, the following specific experiments were performed:

[0098] The collected soil samples were first air-dried and then sieved through a 2mm sieve. The sieved soil was then stored in a dry air environment with a relative humidity below 15% for several months until its mass remained constant. The residual water content was then determined by gravimetric analysis, and the residual water content θ was calculated based on a predetermined bulk density of 1.4. r =0.025cm 3 cm -3 Subsequently, soil samples were slowly and layered into a ring cutter with an inner diameter of 50.46 mm and a height of 50 mm, with a nylon mesh at the bottom. During the filling process, the ring cutter was gently tapped and the layers were roughened to achieve a uniform bulk density. After the ring cutter filling was completed, the initial moisture content θ of the soil sample was... i =0.025cm 3 cm -3 .

[0099] For the completed ring cutter, the cumulative infiltration curve (I~t) is measured using a rapid soil hydraulic property measurement device, specifically using a constant water head h. p =0cm. Infiltration data per second was collected using a portable infiltration device, and the time t for the wetting front to reach the surface was recorded using a camera. * =3356s, corresponding to the cumulative infiltration I when the wetting front reaches the soil surface. * = 1.848 cm. Record the cumulative infiltration I when the infiltration curve tends to stabilize. end =1.984cm, calculate the saturated water content θ s =0.422cm 3 cm -3 Measure the height L of the ring cutter. * =5cm.

[0100] The objective function (formula (14)) is constructed using a genetic algorithm (GA) based on R language. Specifically, the population size is set to 1000, the maximum number of iterations is 5000 (or 500 generations), and a local search strategy is enabled; the parameter search range is K. s =[10 -6 10 -2 ]cm·s -1 1 / h d =[0.01,1]cm -1 , n = [0.1, 0.7], for K s With 1 / h d Perform a logarithmic transformation; transform θ s θ i θ r I * L * and h p As a known parameter, I in formula (14) is calculated by combining formulas (6)-(11). sim (L * ;n,1 / h d ,K s ) and I sim (t i ;n,1 / h d ,K s The constraint threshold ε was set to 0.4 mm. Finally, the soil hydraulic parameters n and 1 / h were obtained by optimizing the calculation formula (14). d and K s The best estimate.

[0101] Figure 3 The parameter iterative estimation process of PASCP and ASCP methods is shown respectively. ASCP method is used as a comparison reference to evaluate the improvement of PASCP method in parameter estimation robustness. Figure 3 The study demonstrates that the PASCP method has a wider parameter search range, while the ASCP method, due to its strict equality constraints, has a limited search space. Interestingly, the parameter values ​​ultimately determined by the ASCP method overlap with some parameter values ​​searched during the iteration process of the PASCP method. This phenomenon indicates that the parameter search space of the PASCP method encompasses the search range of the ASCP method, highlighting the robustness of the PASCP method in the presence of large measurement errors. Figure 3 b further shows that the SWRC predicted by the PASCP method is in better agreement with the measured value (RETC), while the prediction by the ASCP method shows a significant deviation (RMSE_ASCP=0.030, RMSE_PASCP=0.015, RMSE_RETC=0.010).

[0102] Example 4

[0103] This embodiment aims to broadly verify the accuracy and parameter identification ability of the present invention in soil hydraulic parameter estimation. Five soils with significant differences in texture were used as research objects (as shown in Table 1, the properties of the five soils include: soil texture composition (percentage content of sand, silt and clay), soil organic carbon content (SOC), bulk density (ρ). b ), soil residual water content (θ) r ), soil saturated water content (θ) s ), initial soil moisture content (θ) i ) and the cumulative infiltration at a single point when the wetting front reaches the soil surface (I * The data in the table are measured using standard methods (Klute, 1986), and are compared with the ASCP method to highlight the advantages of this invention. Since the specific measurement steps for the five soil types are the same as in the aforementioned embodiments, only the final results are shown here.

[0104] Figure 4 This study demonstrates the dispersion between the estimated and measured parameters of the PASCP and ASCP methods for the five soil types mentioned above. The aim is to verify whether the PASCP method can accurately capture the parameter variation trends across different soil textures, and to show the stronger parameter estimation robustness of the PASCP method compared to the ASCP method. Regression analysis results show that the parameters estimated by the PASCP method (n, 1 / h) are... d , and log 10 (K s There is a significant linear correlation between the measured parameters and the actual parameters (R0). 2 The results (≥0.804, P<0.05) indicate that the PASCP method can effectively reflect parameter changes caused by soil texture differences. Among them, the estimated values ​​of parameters n and log10(Ks) are closer to the measured values ​​(Slope_n=0.866, Slope_log10(Ks) ≥0.804, P<0.05), indicating that the PASCP method can effectively reflect parameter changes caused by soil texture differences. Specifically, the estimated values ​​of parameters n and log10(Ks) are closer to the measured values ​​(Slope_n=0.866, Slope_log10(Ks) ≥0.804, P<0.05). 10 (K s ()=0.803); while parameter 1 / h d The estimated values ​​show some bias, which may be related to the hysteresis effect present in the measured soil moisture characteristic curve. In contrast, the ASCP method exhibits a significant decrease in parameter estimation accuracy under measurement error conditions, and the correlation and statistical significance between the estimated and measured values ​​are also significantly reduced. This further highlights the better stability and reliability of the PASCP method under measurement error interference conditions.

[0105] In summary, the PASCP method can not only effectively capture the parameter variation patterns among different soil textures, but also maintains excellent robustness and reliability even in the presence of actual measurement errors.

[0106] Table 1 Properties of five types of soil

[0107]

[0108] Example 5

[0109] This embodiment has two main objectives: first, to evaluate the soil saturation time (t) introduced in this invention. s The impact on parameter estimation accuracy; secondly, to compensate for the lag effect in the actual measured moisture characteristic curve, which leads to inaccuracies in parameter 1 / h. d The estimation bias may underestimate the advantage of this invention in terms of parameter estimation accuracy. Therefore, this embodiment uses a HYDRUS-1D-based numerical simulation, selecting seven virtual soils with significant texture differences as the research objects (see Table 2; the properties of the seven soils include: soil residual water content (θ)). r ), soil saturated water content (θ) s ), initial soil moisture content (θ) i ), Pore distribution index n, and reciprocal of intake suction 1 / h d and saturated hydraulic conductivity K s Each soil simulation was repeated ten times to assess the dispersion of the estimated parameters. The measurement procedures used in the simulations were consistent with those in the aforementioned embodiments; only the simulation results are shown here.

[0110] Figure 5 The accuracy of parameter inversion using the PASCP method on seven different textures of virtual soil is demonstrated. A saturation time (t) is introduced. s When used as additional inversion constraint information ( Figure 5 a) Parameters estimated by PASCP (n, 1 / h) d and K s The relationship between the actual value and the theoretical value exhibits a highly consistent linear relationship (R0). 2 >0.99, P<0.001), where the estimated value of parameter n almost perfectly matches the theoretical value. Although parameter 1 / h d and K s It is slightly overestimated, but its bias shows a systematic trend, and PASCP can still accurately reflect the differences in 1 / h between different soil textures. d and K s The pattern of change (R) 2 ≥0.998). Furthermore, the low standard deviations of all inversion parameters indicate that this method can uniquely estimate hydraulic parameters. In contrast, without introducing a saturation time (t... sAs an inversion constraint ( Figure 5 b) The parameter n is significantly underestimated (Slope = 0.691), the estimation accuracy decreases and the standard deviation increases significantly, while the parameter 1 / h... d and K s The prediction accuracy also decreased. This result indicates that the saturation time (t) s As inversion constraint information, it plays an important role in improving the accuracy and uniqueness of parameter estimation. In summary, the algorithm proposed in this invention can accurately and uniquely estimate soil hydraulic parameters.

[0111] Table 2 Properties of Seven Soils

[0112]

Claims

1. A robust calculation method for soil hydraulic parameters based on a one-dimensional upward infiltration test, characterized in that... The steps include: Step A: Apply a constant lower boundary water head to the ring sampler soil sample. h p One-dimensional upward infiltration tests were conducted on homogeneous soil, and the cumulative infiltration volume of the ring soil samples was recorded. I Over time t The relationship between changes, and the moment when the wetting front first reaches the surface of the ring sample. t * and the cumulative infiltration rate when the infiltration curve tends to stabilize. I end ; Step B: Based on the initial volumetric water content of the ring sampler soil. θ i Ring cutter height L * and the cumulative infiltration rate when the infiltration curve tends to stabilize I end The saturated volumetric water content of the soil was calculated. θ s The calculation formula is: Step C: Cumulative infiltration volume obtained in Step A I With time t The relationship between the changes was used to determine the moment when the wetting front first reached the surface of the ring sample. t * Corresponding cumulative infiltration at a single point I *; Step D: Combine the saturated volumetric water content of the ring sample. θ s Initial volumetric water content θ i Residual volumetric water content θ r The cumulative infiltration at a single point when the wetting front first reaches the soil surface I * and observational data on the change of cumulative infiltration over time. I obs ( t i The parameters to be estimated are obtained using a parameter optimization inversion method based on interval constraints. n 1 / h d and K s ; In step D, the objective function of the interval-constrained parameter optimization inversion method is constructed according to formula (14): ,in I obs ( t i The measured value obtained in step A represents the observation time point. t i Cumulative infiltration at time point I , I sim ( t i ; n , 1 / h d , K s (at the observation time point) t i At that time, the simulated cumulative infiltration volume was obtained. I sim ( L * ; n , 1 / h d , K s ) represents the cumulative infiltration amount calculated in the simulation.

2. The robust calculation method for soil hydraulic parameters according to claim 1, characterized in that: I sim ( t i ; n , 1 / h d , K s The following formula is used to calculate: ① Soil saturation time t s The function formula is: In the formula z s * It is an equivalent moistened front z fe The function expression, when the equivalent wetting front z fe Reaching the soil surface, i.e. z fe = L * At this point, the soil reaches saturation. z s * This can be explicitly stated as: ②Regarding infiltration time t < t s At that time, the cumulative infiltration rate of the soil I With time t The functional relationship is as follows: ③When t ≥ t s When the soil reaches complete saturation, the cumulative infiltration remains constant, and its calculation formula is as follows: I = (θ s -θ i )L * (10); The parameters to be estimated n 1 / h d and K s Substituting into it, we get the measured value I obs ( t i ) and simulated values I sim ( t i ; n , 1 / h d , K s The root mean square error between them tends to be minimized; I sim ( L * ; n , 1 / h d , K s The following formula is used to calculate: in z f = L *, where σ is the penalty factor. This is a constraint threshold, measured in mm, and its value depends on the upper limit of the maximum measurement error in actual observations. The remaining parameters are expressed according to the following relationship: ； Where parameters a and b The definition is as follows: (2) (3) S ei To correspond to the initial volumetric water content θ i Effective saturation.

3. The robust calculation method for soil hydraulic parameters according to claim 2, characterized in that: Formula (14) is solved using a genetic algorithm in R language: the algorithm generates 1000 sets of parameters to be estimated in each generation. n 1 / h d and K s The algorithm stops running when the value of φ in formula (14) no longer increases within 500 consecutive iterations, or when the cumulative number of iterations reaches 5000; a local search strategy is enabled; the parameter search range is... K s =[10 -6 , 10 -2 ] cm·s -1 , 1 / h d =[0.01, 1] cm -1 , n =[0.1,0.7], for K s With 1 / h d Perform logarithmic transformation; constrain threshold ε Set to 0.4 mm; σ is set to , where .Machine$double.xmax represents the largest double-precision floating-point number that can be represented in R.

Images