Estimating Near-Surface Soil Hydraulic Properties through Sensor-Based Soil Infiltrability Measurements and Inverse Modeling

Abstract

Near-surface soil hydraulic properties (SHPs) are fundamental for describing and predicting water and energy exchange processes, particularly at the soil–atmosphere interface, and regulating evapotranspiration, infiltration, and runoff in different ecosystems. In this study, a new method was proposed to estimate near-surface SHPs by combining sensor-based soil infiltrability measurements with inverse modeling using HYDRUS-2D. The infiltration rate (IR) was estimated by combining the linear source inflow method with image processing, and K_s was estimated from the near-surface steady-state IR (NSIRM). The SWRC parameters described by the van Genuchten model were estimated using the inverse modeling method of HYDRUS-2D for the fitting of sensor-measured infiltration data. Subsequently, the parameters of the van Genuchten model, including $\alpha$, $n$, and $l$, were inversely estimated. Three undisturbed soils, including two stand humus samples from cork oak (Quercus suber L.) and oleander (Pinus tabuliformis L.) stands and one sandy loam from a farmland, were sampled near the soil surface to validate the proposed method. The estimated K_s was evaluated by the constant head method (CHM). The estimated parameters of the SWRC were validated by those determined through the simultaneous measurement of the soil moisture content and water potential using sensor techniques. The results showed that the K_s estimated from the NSIRM for each soil sample were 23.40 ± 1.21, 23.86 ± 1.83, and 22.99 ± 2.26 mm h⁻¹, respectively. In comparison, the K_s determined by the CHM were 24.41 ± 1.53, 24.26 ± 0.37, and 23.81 ± 0.10 mm/h, respectively. The relative errors of the proposed method were 4.14%, 1.64%, and 3.42%, respectively. For the SWRC estimation, the normalized root mean square errors (NRMSEs) between the measurements and the estimates for each soil sample were 0.1724, 0.1454, and 0.0606, respectively. Based on this, the AWC was obtained, and K_u was deduced from the estimated K_s and SWRC parameters for each soil sample. In general, the proposed method successfully estimates near-surface SHPs, simplifies the measurement device, and provides a new perspective for the in situ determination of near-surface SHPs under field conditions in the near future.

1. Introduction

Knowledge of soil hydraulic properties (SHPs) is essential to assess the water cycle in ecosystems [1]. This includes knowledge regarding the infiltration rate (IR), soil water retention curve (SWRC), available water capacity (AWC), and hydraulic conductivity (K), which are all important parameters for determining the quality of the soil environment and the capability of the soil to serve the ecosystem. In contrast to other land uses, studies reveal that forest soils usually show higher hydraulic conductivity [1,2,3]. Moreover, the infiltration capacity and soil water retention in forested soils control the formation of surface runoff, acting as a natural flood regulator [3,4]. Compared to deeper soil layers, the SHPs near the soil surface are subject to rapid changes in response to changes in soil properties, interactions with water (rainfall, infiltration, evaporation), and biopores created by the decomposition of the root. Thus, near-surface SHPs are fundamental for describing and predicting water and energy exchange processes, particularly at the soil–atmosphere interface, regulating evapotranspiration, infiltration, and runoff.

In recent years, studies have witnessed remarkable progress in determining the saturated hydraulic conductivity (K_s), unsaturated hydraulic conductivity (K_u), soil water retention curve (SWRC), and available water capacity (AWC) of SHPs; however, the infiltration rate (IR) has rarely been investigated. In other words, it is difficult to obtain all SHPs from a single method based on existing studies.

Diverse commercial sensors, devices, and methods have been developed and applied under both laboratory and field conditions to determine K_u directly or indirectly. For instance, previous studies have applied a method to determine K_u by measuring soil water potential (ψ) at two soil depths using the Buckingham–Darcy equation during an evaporation process [5,6,7,8,9]. However, this method requires the simultaneous application of multiple sensors, which not only increases the cost, but also decreases the measurement accuracy, owing to sensor-to-sensor uncertainties [10,11]. Consequently, the indirect method has been applied to determine K_u using the measured K_s [12,13,14,15,16,17], and the fitted parameters of SWRC based on hydraulic models [14,18].

K_s can be determined using a single-ring permeameter [19], a Guelph permeameter [20], or a dual-ring permeameter in the field, whereas the falling method or the CHM is employed in the laboratory [21,22,23,24]. The SWRC can be measured using the ceramic pressure plate extractor [25]. However, due to the difference in contact area between the soil and the plate, the measurements of the pressure plate lack reliability in the dry range [26]. Other existing methods for SWRC estimation include fractal geometry, pedo-transfer functions (PTFs) [27], soil morphology, and empirical models.

Mathematical modeling plays a crucial role in estimating SHPs, thus providing valuable insights into water movement in the subsurface. Various methods and tools have been developed to capture the complex relationships between soil properties and hydraulic behavior. Here, we delve into different approaches employed in mathematical modeling for hydraulic property estimation. Numerical simulation models, such as finite element or finite difference methods, allow for a more detailed representation of the soil profile, like HYDRUS, SWAP, and MODFLOW. Numerical models demand computational resources and data-intensive input, making them challenging for certain applications [28]. Pedo-transfer functions (PTFs) are empirical equations that estimate SHPs based on easily measurable soil characteristics. These functions leverage statistical relationships derived from extensive soil databases. PTFs provide a quick and practical way to estimate SHPs, but may have limitations in capturing site-specific variations like soils under conservation or conventional agriculture [29,30]. Recent advances in machine learning (ML) have introduced innovative techniques for estimating hydraulic properties. ML algorithms, including support vector machines, neural networks, and random forests, can learn complex patterns from large datasets. These data-driven approaches offer flexibility and adaptability to diverse soil conditions. But ML approaches may be considered “black-box” models, making it challenging to interpret their internal workings.

In fact, most studies have only measured K_s, K_u, and SWRC, and only a few have measured the IR and AWC. AWC is defined as the ability of soil to store water between the field capacity (FC) and the permanent wilting point (PWP). FC refers to the amount of water retained in the soil after any excess water has drained out, leaving the soil saturated but not waterlogged. It represents the maximum amount of water the soil can hold against the force of gravity and is available for plant use. PWP is the soil moisture content at which plants are unable to extract more water, leading to irreversible wilting. It is the lower limit of plant-available water in the soil, indicating the point at which the soil water potential becomes too low for plants to sustain normal growth [31,32]. Veihmeyer and Hendrickson [33] used sensor measurements, in situ moisture content monitoring, and the laboratory measurements of soil cores in pressure chambers to directly acquire the AWC. It is worth mentioning that the FC is an empirical concept whose physical definition remains under discussion [34]. Th FC is generally defined as the soil moisture content (θ) retained at ψ levels approximately between −10 and −33 kPa, whereas the PWP represents the SWRC at −1500 kPa [35]. Hence, the FC and PWP are defined as the θ at the −33 kPa and −1500 kPa, respectively, based on the estimated SWRC in this study. The above-mentioned parameters are associated, except for the IR, which can be estimated using the linear source method [36] and point source method [37].

This study proposes a new method that uses a simplified device for estimating near-surface SHPs. Assuming soil samples are homogeneous, a steady flow rate is provided to the soil during the experimental process. The soil wetting area is collected via the use of image sensors, while the soil moisture information is obtained using moisture sensors. Both multi-sensor measurements and inverse modeling methods were used to estimate the IR, K_s, SWRC, K_u, and AWC. The SHPs estimated by the proposed method are validated against those estimated by traditional methods. Upon ensuring the accuracy of this method, it is demonstrated that only two sets of data, namely the wetting area and soil moisture information, combined with an inversion module, are required to obtain SHPs.

2. Materials and Methods

2.1. Mathematical Models for the Estimation of the Near-Surface Soil Hydraulic Properties

2.1.1. Linear Source Inflow Method for Estimating Soil Infiltration Rate and Saturated Soil Hydraulic Conductivity

Assuming the soil is uniform and non-layered, when the soil moisture content near the soil surface saturates during the water infiltration process, the soil IR approaches a constant value (steady-state IR), which is commonly considered to be the K_s value. To determine the soil IR, a linear source inflow method proposed by Mao et al. [36] was combined with the camera image processing method to estimate K_s. Assuming a near-uniform soil texture and a sufficient water supply, the soil water infiltration process can be calculated using the relationship between the soil-wetted area and time (Figure 1a). Initially, the wetting area is near zero, but has the highest IR. Subsequently, the IR and the increase in the wetting area gradually decrease. When the wetted area and the IR reached the steady-state value, K_s was estimated from the steady-state IR.

However, during the process of surface water flow, the starting time of the IR at different points on the surface is different. Figure 1b illustrates the distribution of the IR at different spatial locations and times. At time t₁, the water has just reached spatial point A₁, where the soil IR is the maximum initial value, denoted as i₁. Subsequently, at time t₂, the water has just reached the spatial point A₂, where the IR remains at the initial value of i₁. However, at this time, the IR performance at the spatial point A₁ has decreased, as indicated by the curve in Figure 1b, which has decreased to i₂. As time progresses to t₃, the IR at points A₁ and A₂ have decreased to i₃ and i₂, respectively, following the trend shown in Figure 1b.

The relationship among the soil IR, the wetted area in the horizontal direction, and the time during the infiltration process can be obtained using the following equation [36]:

$$q={\int }_0^{A}i\left(A,t\right)dA$$

(1)

where $q$ is the in-flow rate ($\mathsf{\mu }\mathrm{L} {\mathrm{h}}^{-1}$)$, i$ is the infiltration rate ($\mathrm{m}\mathrm{m}{\mathrm{h}}^{-1}$), and A is the wetted area in the horizontal direction (${\mathrm{m}\mathrm{m}}^2$). In fact, directly obtaining the analytical solution of Equation (1) is challenging. Therefore, we divide time into infinitesimally small intervals, and then select the corresponding soil wetting area increment steps (∆A) for each interval. Within each time interval and ∆A, we consider the IR as the average value for that segment. At time t₂ (Figure 1b), with a wetting area increment of ΔA₂, the infiltration rate is i₁, consistent with that observed on ΔA₁ at time t₁. However, the IR on ΔA₁ decreases to i₂ at this moment. Similar situations occur in other time intervals. Eventually, we iteratively derive an approximate estimate of i_n. The specific computational process is described as follows:

According to the principle of water flow balance at time $t_1$ [36],

$$q_1=i_1\cdot ∆A_1$$

(2)

at time $t_2$,

$$q_2=i_2\cdot {∆A}_1+i_1\cdot ∆A_2$$

(3)

at time $t_n$,

$$q_n=i_n\cdot ∆A_1+i_{n-1}\cdot ∆A_2+\cdots +i_2\cdot ∆A_{n-1}+i_1\cdot ∆A_n$$

(4)

From Equations (2)–(4), the $i_n$ at different times is given by the following equation:

$$i_n=\frac{q_n-{\sum }_{j=1}^{n-1}{i}_j\cdot ∆A_{n-j+1}}{∆A_1}\left(n=\mathrm{1,2},3\dots \right)$$

(5)

Finally, K_s was estimated from the steady-state $i_n$.

2.1.2. Inverse Estimation of the Soil Water Retention Curve Using HYDRUS-2D

A two-dimensional soil hydrodynamic model was used to describe the soil water vertical infiltration and horizontal diffusion processes. Assuming that the soil properties and texture are identical near the soil surface and that the influence of the roots can be ignored without vegetation cover, the model can be described by Richards’ equation:

$$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial x}\left(K_u\frac{\partial \psi }{\partial x}\right)+\frac{\partial }{\partial z}\left(K_u\frac{\partial \psi }{\partial z}\right)+\frac{\partial K_u}{\partial z}$$

(6)

where $\theta$ is the soil moisture content $({\mathrm{m}}^3 {\mathrm{m}}^{-3})$, $t$ is the infiltration time $\left(\mathrm{m}\mathrm{i}\mathrm{n}\right),$ $\psi$ is the soil matric potential ($\mathrm{k}\mathrm{P}\mathrm{a}$), $x$ and $z$ are the horizontal and vertical diffusion (infiltration) distances ($\mathrm{m}$), respectively, and $K_u$ is the unsaturated soil hydraulic conductivity.

The relationship between $\theta$, $\psi$, and $K_u$ in Richards’ equation are given by the van Genuchten [10] (VG) model:

$$\theta \left(\psi \right)=\left\{\begin{array}{c}{\theta }_r+\frac{{\theta }_s-{\theta }_r}{{\left[1+{\left|\alpha \cdot h\right|}^n\right]}^m} \mathrm{h}<0\hfill \\ {\theta }_s \mathrm{h}\ge 0\hfill \end{array}\right.$$

(7)

$$K_u=K_s\cdot S_e^l\cdot [1-{\left(1-S_e^{1/m}\right)}^m{]}^2$$

(8)

$$S_e=\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}=\frac{1}{{\left(1+{\left|\alpha \cdot h\right|}^n\right)}^m}$$

(9)

where ${\theta }_r$ and ${\theta }_s$ denote the residual and saturated moisture content (${\mathrm{m}}^3{\mathrm{m}}^{-3}$), respectively, $S_e$ is the effective saturation, $\alpha \left({\mathrm{m}}^{-1}\right)$ is the inverse of the air-entry value (or bubbling pressure) relating to pore-size distribution, $l$ is a pore-connectivity parameter assumed to be 0.5 as an average for many soils [13], and $m=1-n^{-1}$ [14].

To estimate the SWRC, the inversion module of HYDRUS-2D was utilized, with the Levenberg–Marquardt (LM) algorithm being selected to optimize the parameters during the inversion process. The optimization principle of the LM algorithm is to minimize the difference between the measured and simulated variables (e.g., ψ and θ). The sum of the squares of these differences is represented by the objective function $\Phi$, which is expressed as

$$\Phi \left(\theta ,\alpha ,n\right)=\sum _{i=1}^N{\left({\theta }^{\mathrm{*}}\left(t_i\right)-\theta \left(t_i,\alpha ,n\right)\right)}^2$$

(10)

where $\alpha$ and $n$ are the empirical parameters of the VG model, ${\theta }^{\mathrm{*}}$ is the measured value of the moisture content of the soil vertical profile $\left({\mathrm{m}}^3 {\mathrm{m}}^{-3}\right)$, $\theta$ is the soil moisture content when the parameter $\alpha$, and $n$ takes different values $({\mathrm{m}}^3 {\mathrm{m}}^{-3})$.

2.2. System Design

The system design (Figure 2) included the water supply part (Mariotte bottle, linear source water distribution device), soil moisture measurement device, image recording part (camera), a rectangular acrylic container, and a laptop. This camera model is SD-200 (Shenzhen Micro Technology Co., Ltd., Shenzhen, China), with an effective resolution of 1920 H*1080 V (2 megapixels). The outlet of the Mariotte bottle was connected to a linear-source water distribution device in order to supply water at a constant rate. The linear source water distribution device uses a copper adjustable diverter valve, placed directly above the soil surface to ensure that water linearly flows into the soil. The soil sample was filled into an acrylic container (1 m × 0.3 m × 0.2 m). The water supply rates were 27.33 ± 6.17 mL min⁻¹, 31.33 ± 11.67 mL min⁻¹, and 33.833 ± 8.0 mL min⁻¹ for the soil from the cork oak (Quercus suber L.) stand, the soil from the oleander (Pinus tabuliformis L.) stand, and the sandy loam from a farmland, respectively. Commercial soil water potential sensors Teros 21 (METER Group, Inc., Pullman, WA, USA) were used to validate the estimated SWRC.

Based on the study conducted by Yu et al. [38], we designed the moisture sensor, and the diagram and photograph are shown in Figure 3. To calibrate the sensor, we used three types of soil samples which were required for the experiment: the humous soil from both the cork oak and oleander stands, and the sandy loam from the farmland. All soil samples were dried at 105 °C for 48 h, and then moistened with different water contents. Under laboratory conditions at room temperature (25 ± 1 °C), the mixed soil samples were filled into cylindrical PVC containers. The sensor surrounded by the soil samples was placed at the center of the PVC container.

Table 1. The relations and R² values in the different soil samples: the humous soil from both the cork oak (Quercus suber L.) and oleander (Pinus tabuliformis L.) stands, and the sandy loam from the farmland. In which “y” represents the soil moisture content and “x” represents the output voltage value of the moisture sensor.

Soil Samples		Relations	R2
Humus soil	Cork oak stand	$y=22.692x^3-27.274x^2+11.809x-1.7261$	0.9898
	Oleander stand	$y=-7.5084x^3+10.08x^2-2.8212x+0.1402$	0.9986
Sandy loam	Farmland	$y=-186.67x^4+383.46x^3-2877.77x^2+94.603x-11.533$	0.9887

displays the calibration results of the moisture sensor in the three soil samples, with the fitted correlation R² between the sensor output and the soil moisture content ranging from 0.9887 to 0.9986.

2.3. Experimental Preparation and Procedure

2.3.1. Experimental Materials

The experimental materials used were undisturbed soil samples that were collected from two forest stands and a farmland. The humus soil samples were obtained from the surface high-organic-matter soil in the cork oak stand area and the oleander stand area of the Jiufeng National Stand Park (39°54′ N, 116°28′ E), Haidian District, Beijing, with initial moisture contents of 1.42%, and 22.93%, respectively. The sandy loam sample was taken from the 0–10 cm soil layer of San Qingyuan (40°0′ N, 116°20′ E) from the Beijing Forestry University, with an initial moisture content of 2.83%. The soil colors were determined according to the Munsell chart, with the cork oak stand soil being 5GY6/1, the oleander stand soil being 7.5Y8/4, and the sandy loam being 10Y7/2. The experiments for each sample were repeated three times. The physical properties of the soil samples are listed in

Table 2. Physical properties of the experimental soil samples.

Soil Texture		Depth (cm)	Initial Moisture Content (%)	Sand (%)	Silt (%)	Clay (%)	Organic Matter Content (%)	Bulk Density (g cm3)	Porosity (%)
Humus soil	Cork oak stand	0–10	1.42	/	/	/	8.35	1.18	55.47
	Oleander stand		22.93				14.95	1.071	59.58
Sandy loam	Farmland		2.83	39.9	46.6	13.5	1.85	1.33	50.19

.

2.3.2. Experimental Procedure

In the measurement system, a hole was drilled in the undisturbed soil samples for the in situ installation of the PVC pipe. At the beginning of the measurement, the sensor probe was controlled to descend to different depths, measuring the soil moisture content and time.

Using the linear source inflow process, images of water diffusion at the soil surface at different times were captured using a camera. Image processing technology was employed to calculate the wetting areas (A) from the pictures using OpenCV. It includes three steps: distortion correction, binarization, and counting the number of pixels. Distortion correction refers to the phenomenon in which the actual and ideal images differ because of the lens itself. The captured images were cropped to obtain the mesh corner point information to calculate the camera correction coefficient. Binarization was used to set the gray values of the pixels to 0 or 255 for each image. After binarizing the images, the number of pixels in the mesh image was determined using OpenCV. Finally, the wetting area was calculated according to the number of mesh pixels and individual mesh areas.

The inversion module of the HYDRUS-2D model was used to estimate the SWRC. In managing the experiment in 2D, we considered both the vertical (z) and horizontal (x) directions. Despite the horizontal wetting area extending in both the x and y directions, we opted for the linear source infiltration method and assumed homogeneity within the soil samples. Hence, for simplicity, we focused solely on the x direction. Consequently, in HYDRUS, we constructed a soil model with dimensions of 30 cm in the x direction and 20 cm in the z direction, divided into 1200 triangular grids. Furthermore, the inversion module was selected to simulate the water flow dynamics of the model within 120 min, with the estimated results (i.e., output data) set as the soil hydraulic parameters. The soil moisture content was used as the initial condition for the model, with initial moisture contents of the two humus soils and the sandy loam being 0.015 cm³ cm⁻³, 0.24 cm³ cm⁻³, and 0.17 cm³ cm⁻³. Observation points were selected at depths ranging from 2 to 5 cm, with a depth interval of 1 cm, and the soil moisture sensor was used to measure the soil moisture changes at four depth points, providing input data. Additionally, the saturation moisture content (θ_s) in the soil hydraulic parameters was measured from the sensor in saturated soil conditions, while the volumetric moisture content (θ_r) was estimated from the soil moisture characteristic curve corresponding to −1500 kPa. The K_s was obtained using the NSIRM. The remaining soil hydraulic parameters, α and n, were the output data. The upper boundaries of the three soil samples were set at flow rates of 1.64 L h⁻¹, 1.85 L h⁻¹, and 2.06 L h⁻¹, respectively, to simulate atmospheric rainfall boundaries, while the lower boundary was set as a free drainage boundary.

The AWC was estimated using the SWRC. The FC and PWP of the three soils were taken from the moisture contents corresponding to −33 kPa and −1500 kPa in the SWRC, respectively. The K_u value can be calculated using the above parameters with Equation (8). A flowchart for estimating the near-surface SHPs is shown in Figure 4.

2.4. Method for Validation and Error Analysis

2.4.1. Validation of the Saturated Hydraulic Conductivity Using the CHM

The measurement error of K_s, estimated by the NSIRM (Figure 5), assuming the constant head measurement as the true value, can be calculated as follows:

$$K_s=\frac{Q\cdot L}{A\cdot \mathsf{\Delta }h}$$

(11)

where $Q$ is the amount of water flowing out of the test container (${\mathrm{m}\mathrm{m}}^3 {\mathrm{h}}^{-1}$), $L$ is the length of the soil column in the measurement section ($\mathrm{m}\mathrm{m}$), $Ac$ is the cross-sectional area of the soil column in the test container (${\mathrm{m}\mathrm{m}}^2$), and $\mathsf{\Delta }h$ is the height difference between the lowest point of the conduit in the Mariotte bottle and the water outlet of the test container ($\mathrm{m}\mathrm{m}$). The measurement error can be calculated using Equation (12).

$$\delta =\left|\frac{I_2-I_1}{I_2}\right|\times 100\%$$

(12)

where $I_2$ is the K_s measured using the CHM, and $I_1$ is the K_s value measured using the NSIRM.

2.4.2. Comparison of the SWRC and the Measured Results

The HYDRUS-2D software (v2.04) simulation was quantitatively validated using the normalized root mean square error (NRMSE) as follows:

$$\mathrm{RMSE}=\sqrt{\frac{{\sum }_{i=1}^N{\left(P_i-O_i\right)}^2}{N}}$$

(13)

$$\mathrm{NRMSE}=\frac{\mathrm{RMSE}}{O_{max}-O_{min}}$$

(14)

where $N$ is the total number of observation points of θ or ψ, $P_i$ and $O_i$ are the simulated and measured values of stand θ or ψ of the ith observation point, respectively, and $O_{max}$ and $O_{min}$ are the measured maximum and minimum values of θ or ψ at the observation point, respectively. The root mean square difference represents the fitting degree between the simulated and actual values. NRMSE is the normalized RMSE, used to facilitate comparisons.

3. Results

3.1. Estimation Results of the Soil Wetting Area Using Image Processing

After the experiment begins, the camera captures images of the soil wetting situation every 2 min. Here, we have chosen one example to explain the image processing process. Figure 6 shows the image processing results. A checkerboard (Figure 6a) was used to calibrate and correct the camera images, and 54 corners were extracted in order to establish the corner point coordinates (Figure 6b). According to the relationship between the corner point, world, and image coordinates, the internal and external parameter matrices and distortion coefficients of the camera were obtained, and the calibration result (Figure 6c) of the camera was 0.0766. The area of one checkerboard is 4.84 cm² with 3600 pixels; thus, one pixel corresponds to an area of 0.134 mm². The internal and external parameters of the camera and the distortion coefficient are shown in Figure 6d, and the corrected results are shown in Figure 6e. Following binarization (Figure 6f), the soil wetting area was calculated according to the correspondence between the pixels and the area.

We recorded the soil wetting area calculated from image processing every 2 min, and fitted it into a curve representing the change over time (Figure 7). In the initial stage, as infiltration begins, the soil wetting area gradually expands as water penetrates the soil and spreads outward. At the onset of the infiltration experiment, the rate at which the increase in the spread of the wetting area is rapid, especially when there are no significant obstacles on the soil surface. As time progresses, the rate of increase in the wetting area starts to decline because the soil has begun to saturate, and any further infiltration of water is limited or slows down. After a period of time, when the infiltration experiment ceases, the soil wetting area reaches a stable state.

The regression equations of the three soil samples according to the relationship between the wetting area and time are shown in

Table 3. Regression equations between the wetting area and time of the three soil samples.

Soil Samples		$\mathbf{Wetting} \mathbf{Area} ({\mathbf{m}\mathbf{m}}^2$)	R2
Humus soil	Cork oak stand	$A=5.93\times 10^2-4.848\times 10^2\cdot e^{-0.02265t}$	0.9938
	Oleander stand	$A=7.599\times 10^2-5.226\times 10^2\cdot e^{-0.02204t}$	0.9947
Sandy loam	Farmland	$A=7.146\times 10^2-5.724\times 10^2\cdot e^{-0.02438t}$	0.9988

. The coefficients of determination (R²) are 0.9938, 0.9947, and 0.9988, respectively, indicating that the regression equations can better describe the water diffusion of the soil surface over time.

3.2. Estimation Results of K_s Using the Estimated Stable IR

The IR of each soil sample (Figure 8) was calculated using Equation (5). Initially, the IR reaches its peak because the soil is not yet saturated, thus allowing water to quickly penetrate into the deeper layers of the soil. As time progresses, the IR gradually decreases, and the speed of water infiltration slows down. When the soil becomes fully saturated, the IR stabilizes, meaning the soil can no longer absorb more water. We defined a threshold value of 0.1 mm min⁻¹, and the soil IR was regarded as K_s when the slope of the IR curve was less than the threshold value. The estimated and measured K_s values are summarized in

Table 4. Results of K_s as measured using the CHM and NSIRM.

Soil Texture		Constant Head Standard (mm min−1)	Linear Source Inflow (mm min−1)	Relative Error (%)
Humus soil	Cork oak stand	24.41 ± 1.53	23.40 ± 1.21	4.14
	Oleander stand	24.26 ± 0.37	23.86 ± 1.83	1.64
Sandy loam	Farmland	23.81 ± 0.10	22.99 ± 2.26	3.42

. The K_s of the humus samples from the cork oak stand and the oleander stand were 23.40 ± 1.21 mm min⁻¹ and 23.86 ± 1.83 mm min⁻¹, respectively, and the sandy loam was 22.99 ± 2.26 mm min⁻¹. The soil texture and bulk density of the test soil column (19 cm long, 8 cm diameter) were the same as those of the drying test. The pressure head is set to 15 cm. When water flows out of the soil column at a constant rate, the increasing rate of the weight of the outgoing water as time elapsed was recorded, which was used to determine K_s. The K_s of the three soil samples calculated using Equation (11) were 24.41 ± 1.53 mm min⁻¹, 24.26 ± 0.37 mm min⁻¹, and 23.81 ± 0.10 mm min⁻¹, respectively. The relative errors of the three soil samples between the two methods were 4.14%, 1.64%, and 3.42%, indicating that the NSIRM can accurately estimate K_s.

3.3. Estimation Results of the SWRC, AWC, and K_u

The optimized parameters α and n in the SWRC for each sample were inversely yielded when Equation (10) was minimized based on the comparisons between the simulated and measured θ (Figure 9a–c). Comparisons between the inverse-estimated and measured SWRC (Figure 9d–f) and the inverse estimation errors (

Table 5. NRMSE of the simulated and measured θ and SWRC of the three soil samples.

Soil Texture		NRMSE of Soil Water Content (θ)		NRMSE of SWRC	Averaged NRMSE of SWRC
Humus soil	Cork oak stand	First	0.1958	0.1489	0.1724
		Second	0.1947	0.1597
		Third	0.2019	0.2086
	Oleander stand	First	0.1216	0.1209	0.1454
		Second	0.1471	0.1286
		Third	0.1688	0.1867
Sandy loam	Farmland	First	0.1337	0.0672	0.0606
		Second	0.1552	0.0726
		Third	0.1542	0.0421

) indicated that the proposed method is feasible for the inverse estimation of the SWRC.

The AWCs (

Table 6. AWC of the three soil samples.

Soil Texture			Field Capacity (FC)	Permanent Wilting Point (PWP)	Available Water Capacity (AWC)	Average of AWC
Humus soil	Cork oak stand	First	0.0227	0.016	0.0067	0.0058
		Second	0.0166	0.015	0.0016
		Third	0.0221	0.013	0.0091
	Oleander stand	First	0.4031	0.0676	0.3355	0.307
		Second	0.3572	0.0534	0.3038
		Third	0.3274	0.0457	0.2817
Sandy loam	Farmland	First	0.3062	0.0422	0.264	0.221
		Second	0.3062	0.113	0.1932
		Third	0.314	0.1083	0.2057

) of the three soil samples were obtained from the SWRC. Based on the estimated K_s and SWRC, the results of the parameter estimation of the three soil samples using the inverse module are listed in

Table 7. Parameter identification results of the three soil samples.

Soil Texture		θr	θs	α	n	Ks (mm min−1)	l
Humus soil	Cork oak stand	0.0147	0.460	0.0986	5.0757	23.40	0.5
	Oleander stand	0.0200	0.486	0.0357	1.6496	23.86	0.5
Sandy loam	Farmland	0.0126	0.340	0.0177	1.5127	22.99	0.5

.

4. Discussion

Among the current direct methods for estimating soil hydraulic properties, some require measuring ψ at two depths as well as soil evaporation to determine K_u [5,6,7,8,9], while others necessitate flow data measurements using a permeameter to determine K_s [19,20]. Additionally, there are methods that rely on a tensiometer to measure the soil suction and on θ to determine the SWRC [25,26]. However, these direct methods have drawbacks, such as extensive data requirements or the disruption of the soil structure. Like the underlying logic of this study, indirect methods like PTFs [27] estimate SHPs based on readily available soil information. However, research has identified up to 31 factors that can affect soil hydraulic properties (particle size distribution, hydraulic, morphological, and chemical properties) [39], indicating that this method requires numerous soil property parameters. Regardless of the method used, it only estimates a single SHP, and, if one aims to simultaneously obtain all SHPs, the required measurement data will accumulate further.

When compared to the traditional methods, the approach that combines sensors and data inversion retains the advantages of real-time monitoring, high accuracy precision, integration with automated systems, and the inherent labor savings of sensor-based measurement methods. Simultaneously, this approach harnesses the comprehensiveness and flexibility of data inversion, thus allowing the integration of diverse data sources and reducing the research costs. Specifically, this study only requires one soil infiltration experiment. Based on the easily obtainable θ and the estimated K_s, the parameters of the VG model can be inversely estimated, and then K_u can be estimated.

While this method can meet the practical measurement needs, it is still constrained by individual proficiency, and there are still some shortcomings. The collection of the soil wetting area requires manual handling, leading to a disconnection from the overall measurement process. The water supply arrangement for the linear infiltration method is also a manual operation, and the flow rate setting is cumbersome with inherent randomness. Additionally, this study is limited to laboratory measurements; in the future, point source experiments could be considered to extend the application of this method to field environments.

5. Conclusions

This study proposed a method for estimating near-surface SHPs using sensor-based soil infiltrability measurements and the inverse modeling of HYDRUS-2D. The near-surface soil infiltration rate was estimated using a linear source inflow method combined with an imaging processing method. The near-surface K_s was estimated using the stable infiltration rate. The estimated K_s values of two humus soil from the cork oak and oleander stands and one sandy loam from a farmland were 23.40 ± 1.21, 23.86 ± 1.83, and 22.99 ± 2.26 mm h⁻¹, respectively. In comparison, the K_s of the three soil samples measured using the constant water head method were 24.41 ± 1.53, 24.26 ± 0.37, and 23.81 ± 0.1 mm h⁻¹, respectively. The relative errors of the two methods were 4.14%, 1.64%, and 3.42%, respectively. The NRMSE of the SWRC for the three soil samples were 0.1724, 0.1454, and 0.0606, respectively. The AWC and K_u were also estimated based on these results. The method presented in this paper introduces a novel approach where SHPs can be obtained using just two devices, namely the image sensor and the moisture sensor, combined with the HYDRUS inversion module. In this process, data collected from the moisture sensor and the image sensor serve as input data for HYDRUS, and the output data from HYDRUS represent the estimated soil hydraulic properties. This represents a novel approach in the field, utilizing fewer devices that are integrated into a single apparatus, allowing for simultaneous data acquisition, while ensuring accuracy in estimation results. Further studies are needed to evaluate the performance of the proposed method under different soil classes.