Estimation of Stagnosol Hydraulic Properties and Water Flow Using Uni- and Bimodal Porosity Models in Erosion-Affected Hillslope Vineyard Soils

Abstract

Erosion has been reported as one of the top degradation processes that negatively affect agricultural soils. The study objective was to identify hydropedological factors controlling soil water dynamics in erosion-affected hillslope vineyard soils. The hydropedological study was conducted at identically-managed Jastrebarsko (location I), and Jazbina (II) and (III) sites with Stagnosol soils. Soil Hydraulic Properties (SHP) were estimated on intact soil cores using Evaporation and WP4C methods; soil hydraulic functions were fitted using HYPROP-FIT software. For Apg and Bg/Btg horizons, uni- and bimodal soil hydraulic models could be well fitted to data; although, the bimodal model performed better in particular cases where data indicated non-uniform pore size distribution. With these SHP estimations, a one-year (2020) water flow scenario was simulated using HYDRUS-1D to compare water balance results obtained with uni- and bimodal hydraulic functions. Simulation results revealed relatively similar flux distribution at each hillslope position between the water balance components infiltration, surface runoff, and drainage. However, at the bottom profile at Jastrebarsko, bimodality of the hydraulic functions led to increased drainage. Soil water storage was reduced, and the vertical movement increased due to modified soil water retention curve shapes. Adequate parameterization of SHP is required to capture the hydropedological response of heterogenous erosion-affected soil systems.

1. Introduction

Anthropogenic activities have a large impact on soil erosion processes of landscapes; increased erosion rates are connected with various land management practices [1]. Erosion by water, wind, and tillage has been reported and the soil erosion phenomena include a broad spectrum of processes [1,2], which come with different characteristics (form, intensity, and frequency), and are present on all continents [3,4].

In vineyards, soil redistribution by both tillage [5] and water erosion [6] contributes to relatively high soil erosion rates during shorter [7] and longer terms [8,9,10]. Sloped vineyards are one of the most intensively managed arable landscapes with complex hydrology due to the effects of erosion, tillage, crop, and compaction caused by trafficking. The combination of these conditions and measures presents a relevant basis for a close-up of particular processes governing water dynamics and redistribution in agroecosystems. Recent studies showed that non-sustainable soil erosion rates have been observed [8,10,11,12,13,14] on the vineyard hillslopes in Mediterranean countries due to the large interrow widths, low plant cover, and steep slopes [15], as well as intensive agrotechnical practices (e.g., tillage).

In arable fields, spatial soil heterogeneity can result from long-term soil management in combination with erosion processes at the sloped areas [16]. Tillage-induced erosion has a significant impact on local-scale soil hydraulic properties [17], whereas trafficking often leads to the formation of impeding subsurface layers due to soil compaction [18]. In erosion-affected soils, topsoil horizons usually vary in thickness depending on their position along the slope, and they can also be significantly altered or completely removed [19,20,21]. These modifications can then substantially affect soil water dynamics and water flow, generating surface runoff [8,14] or subsurface lateral flow [22].

One of the ways to evaluate erosion effects on soil water dynamics is to combine field and laboratory methods with numerical modeling. To achieve accurate modeling results of (un)saturated soil water dynamics, suitable soil parameters, representative of an experimental area, are needed. The soil hydraulic properties (SHP) that are required for the simulation of water flow can be obtained in situ or by using laboratory methods such as a combination of HYPROP [23] and WP4C [24] instruments. They can also be indirectly estimated using statistical models called soil pedotransfer functions (PTFs). The program Rosetta [25] uses basic soil properties, including soil texture, bulk density, and organic matter content, to estimate soil water retention and hydraulic conductivity curves (SWRC & SHCC) and their probability distributions [26]. Experimentally, a combination of HYPROP and WP4C was reported to deliver high-resolution data of soil water retention and hydraulic conductivity over a wide range of soil tensions [27]. Such a data set allows selection of an appropriate hydraulic model for SWRC and SHCC fitting.

The SWRC and SHCC are the most essential functions for the application of physically-based numerical models for simulating water flow and solute transport in soils. The models use either single-, dual-, and even multiple-continuum approaches [28,29,30]. For parameterizing the soil hydraulic properties, models for unimodal, bimodal, or multimodal pore size distributions have been proposed [31]. Two-domain concepts exist for mobile-immobile solute transport and for the non-equilibrium type of preferential flow [32,33]. The latter assumes that the porous soil system can be divided into two superimposed domains, one representing the porous soil matrix (micropore domain) and the other, the fracture or macropore domain, resulting in two sets of coupled flow and transport equations [34]. To compare soil hydraulic parameter models, Haghverdi et al. [35] studied the performance of five unimodal water retention models and eleven combinations of alternative (e.g., Peters–Durner–Iden, PDI) and bi-modal expressions. Overall, Haghverdi et al. [35] showed that alternative expressions provided a better fit than unimodal expressions.

More complex multimodal modeling approaches have a large implication in arable structured soils for capturing heterogeneous pore structures. The most commonly used models in the eroded landscape include models such as WEPP [36], SWAT [37], or RUSLE [38]. However, the HYDRUS model allows close-up estimation of the effects of soil hydraulic parameters. The HYDRUS numerical code has widely been used for modeling water and solute dynamics in the (un)saturated soil zone in one, two, or three dimensions [39]. The model has been successfully used for estimating water flow in hillslope landscapes [22] or in structured arable soils [30] by applying uni- and bimodal soil hydraulic functions and single- and dual-porosity/permeability flow models [28,30]. For example, Dettman et al. [40] used inverse modeling procedure in HYDRUS-1D to obtain soil hydraulic properties in peatland based on the commonly applied van Genuchten–Mualem model and the bimodal model by Durner [31]. They showed that fitting only the macropore fraction of the bimodal model as immediately drainable additional pore space could be a practical approach to account for the macropore effect. A single-porosity model in HYDRUS-1D was sufficient to describe flow and advective diffusive non-reactive transport of Brilliant Blue in undisturbed columns from structured vineyard soil during the evaporation and during outflow experiments [30]. However, dye staining and X-ray imaging revealed a complex pore-architecture network with large vertical and horizontal biopores, which indicated a dual-porosity system with contrasting hydraulic properties. In the hillslope landscape, Rieckh et al. [19] studied water and dissolved C fluxes in four differently eroded soils affected by erosion-induced pedological and soil structural alterations that were previously analyzed on soil cores in the lab [41]. The single porosity Richards’ equation in HYDRUS-1D was used to estimate vertical water movement, while solute fluxes (DOC and DIC) were obtained by combining calculated water fluxes with observed solute concentrations. The spatially distributed leaching results supported the hypothesis that the effects of soil erosion influence fluxes through modified hydraulic and transport properties and terrain-dependent boundary conditions. The above cited HYDRUS code references proved the applicability of numerical modeling, ranging from simple uniform water flow to more complex transport simulations that include soil structural effects on soil water dynamics.

In this study, it was hypothesized that erosion-affected soil structural differences exist in vineyard soils between the top and bottom hillslope positions as observed by Jakšík et al. [42] and Nikodem et al. [43] for arable soils. These structural differences affect the local soil water balance [20], leading to differences in crop yield [44] and can be quantified in the form of bi- or multimodal soil hydraulic functions [45], which were tested through comparing those two functions.

The objective of the present study was an adequate identification of hydropedological factors that influence soil water dynamics [46] in erosion-affected and intensively managed agroecosystems (vineyards). The specific aims were to compare SHP using uni- and bimodal SWRC models to represent the soil pore system at the hilltop and bottom position of three eroded vineyards with similar land management and soil types. The effects of applying uni- versus bimodal soil hydraulic functions on the soil water balance will be studied together with the soil water dynamics at the top versus bottom slope positions for a one-year time period.

2. Materials and Methods

2.1. Experimental Site and Soil Properties

The experimental study was performed in central Croatia (in Zagreb and the nearby southwestern surrounding area). Three selected locations are vineyards located in the same climatic zone (humid continental climate), with the same management practices, i.e., each vineyard is intercropped with grass (2 m wide) parallel with the slope ranging from 10 to 20% (Figure 1). The selected vineyards were Jastrebarsko (45°41′22″ N 15°38′22″ E, noted location I), and two locations in Jazbina (45°51′27.0″ N 16°00′14.4″ E; 45°51′24″ N 16°00′22″ E, noted locations II and III, respectively). The study was performed during 2019 and 2020, focusing on the interrow area, where erosion is the most intense. The samples were collected at the top and the bottom hillslope positions from both the topsoil and the subsoil horizons. Besides the disturbed samples, a triplicate of undisturbed samples was taken as well using soil cores (250 cm³). On each location, soil pits were dug at both hillslope positions, and the soil profiles were described based on the FAO guidelines [47].

According to the WRB classification system [48], all soil profiles were classified into the reference soil group of Stagnosols. These represent the most widespread soils of continental Croatia [49,50]. Although all of the studied soils are classified as Stagnosols, they are not entirely uniform in their properties. Accordingly, principal and/or supplementary WRB qualifiers somewhat vary among the vineyards and between the slope positions (

Table 1. Selected soil properties, i.e., pH (H₂O), organic carbon content (%), particle size distribution (%), bulk density (g cm⁻³), and porosity (%), presented with horizon designation.

Sample Name	Location	Soil Profile	Depth (cm)	pH (H2O)	Organic C (%)	Particle Size Distribution (%)			Bulk Density (g cm−3)	Porosity (%)	Horizon Designation (FAO, 2006)	Soil Classification (IUSS Working Group WRB, 2014)
						Sand (2–0.063 mm)	Silt (0.063–0.002 mm)	Clay (<0.002 mm)
Apg_TOP_I	(I) Jastrebarsko	TOP	0–55	5.5	1.2	4	56	40	1.31	47	Apg	Eutric Protovertic Stagnosol (Aric, Humic, Inclinic, Loamic)
Bg_TOP_I			55–110	5.6	0.4	3	62	35	1.20	50	Bg
Apg_BOT_I		BOT	0–88	5.1	1.1	12	68	20	1.59	45	Apg	Dystric Stagnosol (Aric, Colluvic, Humic, Inclinic, Siltic)
Bg_ BOT_I			88–110	5.2	0.6	10	69	21	1.54	41	Bg
Apg_TOP_I	(II) Jazbina	TOP	0–50	6.2	1.2	11	67	22	1.41	47	Apg	Eutric Luvic Stagnosol (Aric, Humic, Endoloamic, Episiltic)
Btg_TOP_I			50–110	5.7	0.5	15	48	37	1.55	41	Btg
Apg_BOT_I		BOT	0–75	6.9	0.9	8	68	24	1.57	41	Apg	Eutric Stagnosol (Siltic, Aric, Colluvic, Inclinic)
Bg _BOT_I			75–110	5.1	0.8	6	69	25	1.70	36	Bg
Apg_TOP_I	(III) Jazbina	TOP	0–50	5.3	1.2	6	71	23	1.53	45	Apg	Dystric Luvic Stagnosol (Aric, Humic, Endoloamic, Episiltic)
Btg_TOP_I			50–110	5.2	0.2	7	61	32	1.79	36	Btg
Apg_BOT_I		BOT	0–65	5.4	1.0	7	73	20	1.61	43	Apg	Dystric Stagnosol (Aric, Colluvic, Humic, Inclinic, Siltic)
Btg_BOT_I			65–110	5.0	0.5	6	68	26	1.50	42	Btg

). The main differences among the analyzed profiles are as follows: (a) profile I-TOP had some shrink-swell cracks and notably more clay (averagely 37.5%) than other profiles (Table 1); (b) due to leaching, profiles II_TOP and III_TOP showed pronounced vertical increase in clay content (Table 1); (c) due to colluviation, at each location the Apg horizon is deeper at the bottom profile than at the top profile (Table 1). All disturbed soil samples were air-dried, crushed, and sieved through a 2 mm sieve. Soil particle size distribution was determined using the pipette method, with wet sieving and sedimentation after dispersion with sodium-pyrophosphate [51]. The soil pH was measured using the Mettler Toledo MPC 227 conductivity/pH-meter in water (pH H₂O) (HRN ISO 10390:2005). The organic matter content was determined by sulfochromic oxidation (HRN ISO 14235:1998). The selected soil chemical and physical parameters are presented in Table 1.

2.2. Soil Hydraulic Properties Estimation

SHP were determined on undisturbed soil cores (250 cm³; inner diameter 80 mm and height 50 mm) taken from Apg and Bg/Btg horizons at three vineyards (top and bottom hillslope positions). The soil sampling was performed in three replicates by carefully using a hammer and a sample ring insertion tool, after which soil cores were transported to the laboratory and stored at 4 °C till analysis. SHP were estimated by applying the simplified evaporation method [52] and using the HYPROP (METER Group, Pullman, WA, USA) automatized system [53]. The 18 undisturbed soil sample (250 cm³) were saturated from the bottom of the sample. After the saturation phase was completed, two tensiometers were placed inside the undisturbed soil sample on the HYPROP device at two depths and measured the soil water tension [54]. The evaporation method considers the change in the sample weight and matrix potential (i.e., soil water tension) in the soil sample during the evaporation drying process. It allows accurate characterization of the water-retention properties in the porous media [23]. The measurement is finished when the soil sample evaporates and the tensiometers reach their measurement peak (e.g., air entry point). The results allow the representation of SWRC and SHCC [55].

After the determination of data points for water content, hydraulic conductivity was calculated using the water flow velocity ($q_i$[cm d⁻¹]) between time points $t_{i-1}$ and $t_i$ through a horizontal plane in the middle of the two tensiometers as explained in Singh et al. [55]:

$$q_i=\frac{1}{2}\frac{\left(∆V_i/∆t_i\right)}{A}$$

(1)

where $∆V_i$ is the change in water volume in the sample (cm³), $∆t_i$ is the time interval between two consecutive measurement points, and$A$ is the cross-sectional area (cm²) of the column.

The data points for the hydraulic conductivity function were calculated by inverting Darcy’s equation [55]:

$$K_i\left(h_i\right)=\frac{-q_i}{\left\{\left(\frac{∆h_i}{z}\right)-1\right\}}$$

(2)

where $h_i$ (cm) is the time- and space-averaged tension, $∆h_i$ is the difference in tensions between the two tensiometers tips, and z (cm) is the distance between the tensiometer. The calculations assume that moisture tension and water content are distributed linearly through the column and, therefore, the arithmetic mean of the tensions at two points was used. This simplified assumption was shown to provide accurate results because linearity errors in fluxes and tensions cancel each other out [56].

After the HYPROP analysis, the undisturbed soil sample was removed from the base sensor and the tensiometers were carefully removed from the holes inside the sample. The next step was the preparation of subsamples for the WP4C device (METER Group, Pullman, WA, USA). The instrument can be combined with HYPROP to measure water potential in the dry range of the soil sample. WP4C is capable of measuring both matric and osmotic potential (applicable for salinity-affected soils) that completes the SWCC [27] and it is used for the water potential measurements in the dry range of the soil sample.

2.3. Fitting of the SWRC and SHCC Using Uni- and Bimodal Hydraulic Models

The hydraulic functions based on the HYPROP and WP4C measurement data were fitted using HYPROP-FIT software (Version 4.2.2.0). Two most commonly applied soil hydraulic models were used and compared: uni- and bimodal porosity models.

In the first approach, fitting was performed using van Genuchten–Mualem (VGM) model [57,58]:

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

(3)

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

(4)

$$K\left(\theta \right)=K_s{S}_e^l{\left[1-{\left(1-S_e^{\frac{1}{m}} \right)}^m\right]}^2$$

(5)

where h (cm) is pressure head, $\theta$, ${\theta }_r$, and ${\theta }_s$ (cm³ cm⁻³) are the present, residual, and saturated water contents, respectively. Parameter $S_e$ represents effective saturation (cm³ cm⁻³). Parameters $\alpha$ (cm⁻¹), n (unitless) and m (unitless) are the empirical parameters that describe shapes of the SWRC and SHCC, where m = 1 − 1/n. K is the unsaturated hydraulic conductivity function (cm day⁻¹).

In the second approach, a bimodal model [31] was used as well, which can describe a structured soil with developed macropores. Following Durner [31], the porous medium can be divided into i overlapping VGM functions weighted by the factor ${\omega }_i$.

$$S_e={\displaystyle \sum }_{i=1}^k{\omega }_i{\left(\frac{1}{(1+{\left({\alpha }_{i}h\right)}^{n_i}}\right)}^{m_i}$$

(6)

with the sum of ${\omega }_i$ and ${\omega }_k$ being equal to 1. By combining the bimodal retention functions with Mualem’s [58] pore-size distribution model, the bimodal unsaturated hydraulic conductivity can be described with the following equation [59]:

$$K\left(S_e\right)=K_s{\left({\displaystyle \sum }_{i=1}^k{\omega }_i{S}_e{}_i\right)}^l{\left(\frac{{{\displaystyle \sum }}_{i=1}^k{\omega }_i{\alpha }_i\left[1-{\left(1-S_{e_i}^{\frac{1}{m_i}}\right)}^{m_i}\right]}{{{\displaystyle \sum }}_{i=1}^k{\omega }_i{\alpha }_i}\right)}^2$$

(7)

Fitting quality for soil hydraulic parameters estimation is given in terms of the root mean square error (RMSE), which indicates the mean distance between the data point and the fitted function:

$$RMSE=\sqrt{\frac{1}{r}{{\displaystyle \sum }}_{i=1}^r{\left[y_i-{∆}_i\right]}^2}$$

(8)

where y_i and ŷ_i are measured and model-predicted quantities, respectively, i.e., water contents, θ or hydraulic conductivities, log₁₀(K). The model error for water retention (RMSE_θ) was calculated separately from the model error for hydraulic conductivity (RMSE_K).

2.4. Water Flow Modeling

Numerical modeling was performed using the HYDRUS-1D program (version 4.0) [39] that uses the Galerkin-type linear finite element scheme to solve the governing flow equations for uni- and bimodal porosity models. Here, a one-dimensional model was applied for hillslope as the aim was to directly compare top and bottom soil profiles between the three sites (I, II, and III) in terms of uni- and bimodal models and also the soil hydraulic properties responses.

Water flow simulations were based on a numerical solution of the Richards equation for variably-saturated media assuming root water uptake, i.e., grassed interrow strip with the Feddes approach [60] based on the potential uptake parameters according to Wesseling et al. [61]. Richards’ equation is assumed as follows:

$$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}\left[K\left(\frac{\partial h}{\partial z}+1\right)\right]-S$$

(9)

where t is time (T), z is the spatial coordinate (L), and S is the sink term representing root water uptake (L⁻¹), as stated already above.

Root depth was assumed as 15 cm with maximum root density as the grass cover is well developed and maintained for several growing seasons at all three sites. The simulations were performed for six profiles (three locations at top and bottom positions). The profiles were set down to 110 cm depth with the layering according to Table 1 (Apg/Bg/Btg horizon distribution). An atmospheric boundary condition (precipitation and evapotranspiration) with surface runoff assuming evapotranspiration and rainfall was selected at the soil surface, while free-drainage conditions were imposed at the bottom, at the 110 cm depth, as there was no groundwater table effect at selected plots. The initial condition for water flow was set as a pressure head of −100 cm for all profiles in order to have the same conditions for the comparison of the effects of SHP and layering on water dynamics. The daily evapotranspiration was calculated with the Penman–Monteith equation [62] implemented in the HYDRUS-1D based on the meteorological observations collected at Maksimir, Zagreb meteorological station (https://meteo.hr/ (accessed on 21 July 2021)). The simulations were carried out for a one-year period, i.e., the year 2020, using the one-dimensional soil profile with 101 nodes. Numerical simulations assuming uni- and bimodal SHP modeling approach with HYDRUS code [39] for unsaturated soil zone were applied to quantify and visualize differences in water dynamics among three sites for top and bottom positions in our study.

3. Results and Discussion

3.1. Soil Hydraulic Properties

The HYPROP-fitted SWRC for the three locations (I, II, and III) using the uni- and bi-modal water retention models for the first Apg (Figure 2), and for the second Bg/Btg horizons (Figure 3), match the measured data points for most samples, depths, and hillslope locations, except on location II at bottom slope position (Apg_BOT_II), where a larger range between data from three replicates was observed for the first horizon. Relatively large differences in θ_s between Apg and Bg/Btg horizon were estimated at the bottom of location II, although the particle size distribution and organic carbon content were similar. As mentioned above (Figure 2), SWRC for Apg horizon showed a larger data range between three replicates, which indicates the possibility of unreliable hydraulic parameter estimation at this particular location. As the cylinders cover a relatively small area (Ø 80 mm), replicates are necessary, especially in agricultural soils, where it could indicate soil heterogeneity. These differences in replications could be a result of soil erosion which can lead to the soil heterogeneity at the bottom; however, this needs to be further investigated. At location III, θ_s was largest at the top Apg horizon, which also had the highest porosity/lowest bulk density values, especially compared to the top Btg horizon. Although the top Btg horizon has the largest amount of clay particles, θ_s was the lowest, which can be caused by the highest bulk density values i.e., the lowest porosity [63], and being connected to the difference in organic carbon content [64]. Even though the data points (Apg_BOT_II) near saturation are more dispersed compared to the other soil horizons, the data at the dry range (around pF 4.2) are more uniform, thus giving validity to the measurements. This indicates the importance of covering the full range of soil water retention curve measurements. The evaporation-based HYPROP method provided a simultaneous fit of the soil water retention and the hydraulic conductivity functions (Figure 4 and Figure 5). As previously reported, the evaporation method provides high-resolution water retention and hydraulic conductivity data [56], highly applicable for the majority of soils [65,66,67].

The van Genuchten model and its bimodal form curves followed the data agreeably; although, we can clearly notice that bimodal performed better in particular cases where nonlinearity in soil hydraulic data was observed (e.g., Figure 2 and Figure 3; Apg_TOP_II; Bg_BOT_II; Btg_BOT_III). In some cases, fitted unimodal SWRC did not follow the dry range data derived from WP4C measurements (e.g., Figure 2 and Figure 3; Apg_BOT_I; Bg_BOT_I; Bg_BOT_II). The same behavior is reproduced at the SHCC curves, which in some cases showed bimodal character (e.g., Figure 4, Apg_TOP_I, Apg_TOP_II, and Btg_BOT_III). The soil water retention dry range data present an important part of the SWRC fitting as it can largely influence its shape and soil water content at pF 4.2, the wilting point, with similar already being documented [27]. Here, we had extended our wet range data using the WP4C to obtain the retention points near the wilting point, and reproduce more reliable results. For even more reliable results of the SHCC, some authors recommended that the evaporation data should be additionally combined with separate conductivity measurements near saturation by using, for example, tension disc infiltrometers [68].

The final SHP and its statistical evaluation (RMSE) derived from the HYPROP-FIT program using uni- and bimodal soil hydraulic models are presented in

Table 2. Unimodal SHP obtained with HYPROP-FIT software at three locations (I, II, III) at the hillslope (TOP/BOT) with the indication of fitting goodness (RMSE).

Sample Name	θr (cm3 cm−3)	θs (cm3 cm−3)	α (cm−1)	n (-)	Ks (cm day−1)	l (-)	RMSE_θ (cm3 cm−3)	RMSE_K (cm day−1)
Apg_TOP_I	0	0.476	0.00226	1.228	0.0431	−1.94	0.01	0.3181
Bg_TOP_I	0	0.518	0.00212	1.263	0.0879	−1.311	0.0102	0.3184
Apg_BOT_I	0	0.400	0.00845	1.175	1.95	0.5	0.0139	0.1735
Bg_BOT_I	0	0.404	0.00156	1.301	0.0807	2.31	0.0121	0.2812
Apg_TOP_II	0.035	0.378	0.0162	1.184	0.867	−6	0.0168	0.308
Btg_TOP_II	0	0.439	0.0157	1.147	2.67	−1.901	0.0084	0.2537
Apg_BOT_II	0	0.452	0.00641	1.173	0.648	−2.663	0.0248	0.3383
Bg_BOT_II	0	0.356	0.0086	1.131	0.676	−3.243	0.0091	0.3422
Apg_TOP_III	0	0.481	0.00583	1.261	1.07	−2.461	0.0112	0.1505
Btg_TOP_III	0	0.42	0.00367	1.163	0.814	−4.298	0.0097	0.2968
Apg_TOP_III	0	0.478	0.00211	1.3	0.173	−3.409	0.0148	0.2101
Btg_TOP_III	0	0.441	0.00433	1.221	0.291	−3.974	0.014	0.2561

and

Table 3. Bimodal SHP obtained with HYPROP-FIT software at three locations (I, II, III) at the hillslope (TOP/BOT) with the indication of fitting goodness (RMSE).

Sample Name	θr (cm3 cm−3)	θs (cm3 cm−3)	α1 (cm−1)	n1 (-)	Ks (cm day−1)	l (-)	${\mathit{\omega }}_2$ (-)	α2 (cm−1)	n2 (-)	RMSE_θ (cm3 cm−3)	RMSE_K (cm day−1)
Apg_TOP_I	0.044	0.481	0.00124	1.342	0.0493	−4.416	0.051	0.0436	2.975	0.0086	0.2186
Bg_TOP_I	0	0.529	0.00126	1.328	0.14	−3.327	0.057	0.0383	2.96	0.0073	0.2039
Apg_BOT_I	0.033	0.414	0.0699	2.04	5.15	−0.229	0.8	0.00081	1.527	0.0037	0.078
Bg_BOT_I	0	0.41	0.0122	1.122	1.52	2.792	0.391	0.0005	6.999	0.0073	0.2114
Apg_TOP_II	0.085	0.385	0.00314	1.353	0.397	−6	0.183	0.0804	2.406	0.0145	0.2679
Btg_TOP_II	0.125	0.441	0.00019	1.946	1.07	−3.773	0.575	0.0256	1.346	0.0075	0.2569
Apg_BOT_II	0	0.461	0.0512	1.935	0.612	−4.8	0.887	0.00177	1.23	0.0225	0.3632
Bg_BOT_II	0	0.364	0.00016	1.509	3.46	−5.34	0.46	0.0453	1.164	0.0047	0.2603
Apg_TOP_III	0.012	0.489	0.00228	1.343	1.28	−1.852	0.128	0.0319	2.259	0.0063	0.1576
Btg_TOP_III	0.074	0.422	0.00952	1.472	0.375	−5.049	0.691	0.00031	1.355	0.008	0.3035
Apg_TOP_III	0	0.489	0.00123	1.347	0.37	−5.174	0.088	0.0331	1.85	0.0121	0.1985
Btg_TOP_III	0	0.448	0.00152	1.294	0.308	−5.299	0.098	0.0418	2.642	0.0099	0.1317

. Both uni- and bimodal fits of the soil water retention models resulted in relatively similar values of θ_s (0.356 to 0.518 cm³ cm⁻³ for unimodal and 0.364 to 0.529 cm³ cm⁻³ for bimodal, respectively). Similarly, fitted K_s values ranged from 0.0431 to 2.67 cm day⁻¹ for uni- and from 0.0493 to 5.15 cm day⁻¹ for bimodal models, respectively. These soil properties indicate relatively large water retention capacity and low permeability, thus having the potential for surface runoff generation as they are located on the sloped landscape. As far as for model fitting evaluation, soil water content’s RMSE ranged from 0.0084 to 0.0248 cm³ cm⁻³ for the uni-model and from 0.0037 to 0.0225 cm³ cm⁻³ for the bimodal model, respectively. The model fitting for the hydraulic conductivity ranged from 0.1505 to 0.3422 cm day⁻¹ for uni- model and from 0.078 to 0.3632 cm day⁻¹ for the bimodal model, respectively. The RMSE evaluation clearly indicates that both models performed very well in the goodness of fit to observed data. However, particular care should be taken when evaluating SHP and SWRC/SHCC data. In each case, more appropriate soil hydraulic and water flow models should be selected, as selecting less applicable models can further have a significant effect when applying numerical modeling and quantifying hydrological responses [69]. The selection of a less appropriate model, even if the fittings (data vs. model) are similar, can lead to differences in the water balance when simulating longer periods (this is illustrated in Section 3.2). For the undisturbed soil samples collected from multiple irrigated turfgrass plots in southern and central California, Haghverdi et al. [35] reported that the van Genuchten model with five free parameters had RMSE of 0.007 cm³ cm⁻³, while Van Genuchten bimodal with five free parameters had RMSE of 0.004 cm³ cm⁻³. It was concluded that the alternative variants (i.e., PDI and bimodal expressions) outperformed the original unimodal expressions.

With the exception of a few horizons (see above), the water retention did not show the effects of multimodal pore size distributions. Since the HYPROP method focuses on the drying SWRC, the SWRC range close to full water saturation was not fully captured, such that the variations in the range of the structural macro-porosity are hardly described here. The bimodality observed with the present method is related to the textural pore size ranges, i.e., non-uniform pore size distributions caused by textural variations [31,70].

For the SWRC as a capacity function, the structural macroporosity is on the one hand not fully in the range of measurement by the HYPROP, and on the other hand, the pore volume fraction occupied by macropores is often relatively small compared to the bulk pore volume [32] and is difficult to detect considering the measurement limitations. However, for the hydraulic conductivity function as an intensity function, the bimodality is more pronounced because the hydraulic conductivity is highly dependent on pore size and pore connectivity and much less on the pore volume [71,72]. Therefore, the SHCC is multimodal for almost all samples, but notably for Apg_TOP_I, Apg_TOP_II, and Btg_BOT_III (Figure 4 and Figure 5).

3.2. Water Flow Modeling

With the fitted uni- and bimodal SHP, the soil water balance and the soil water dynamics during a one-year period (2020) were simulated for each profile (TOP vs. BOT) at three locations (Figure 6 and Figure 7) using HYDRUS-1D assuming local equilibrium. The simulations were used to visualize the effect of uni- and bimodal soil hydraulic functions in terms of soil water dynamics in the erosion-affected landscape. The modeling results show an appreciable representation of cumulative flux distribution in profiles, i.e., surface runoff, infiltration, and bottom flux. In general, both models provided similar fluxes distribution at each position between infiltration, surface runoff, and drainage (bottom flux), with a few exceptions mostly connected to the positions where nonlinearity in SHP was noticed (e.g., I_BOT and II_TOP positions). The surface runoff amount is directly connected to the predicted value of the saturated hydraulic conductivity, K_s, which was very low (below 3 cm day⁻¹) in most cases, with the exception of K_s of 5.15 cm day⁻¹ for the bimodal model at the bottom profile at location I. This set of SHP (location I bottom profile) with similar infiltration generated considerably larger bottom flux with the bimodal function (Figure 6). This example will be further discussed below. When comparing the top and bottom positions at each site, for the first two locations (I and II), larger surface runoff is generated at the top position while it is the other way around at the third (III) location). This indicates that local scale SHP was of greater importance and had more relevance as compared to fluxes at the hillslope scale. Details of water balance in Figure 7 further illustrate the similarity of flux components and differences at particular locations (e.g., I_BOT).

Figure 8 shows the daily temporal resolution of a bottom flux at the bottom profile at location I, which indicates the largest nonlinearity and bimodality in water retention data (Figure 2 and Table 3). The simulation using the bimodal SHP predicted a larger and highly dynamic bottom flux (Figure 8) as compared to the unimodal representation. The presented difference supports the discussion related to the importance of appropriate hydraulic model selection. Large bottom flux is also connected to the profile horizon depth where Apg is 88 cm deep; thus, its K_s value of 5.15 cm day⁻¹ is a cause of a large bottom flux in the bimodal approach. This large difference is caused by the shape of the two functions (SWRC, and SHCC), which control the relation between smaller values of the soil water capacity and higher hydraulic conductivities. In the case of the bimodal function during drainage, this leads to increased bottom flux. This example demonstrates the relevance of adequate parameterization of soil hydraulic functions in order to capture the hydrological response of erosion-affected soil systems.

The SHCC estimations via the van Genuchten–Mualem model can result in relatively low performance in the near saturation range due to the inability to account for macropore flow [73,74,75]. Furthermore, K_s is a highly variable soil hydraulic property dependent upon the pore geometry at the scale of interest [74] and seasonal variability [75]. Significant variabilities in K_s estimations might occur when using different estimation procedures [76] and measurement techniques [77], ultimately reflecting on the SHCC estimations. Further studies are needed to determine how adding extra near saturation hydraulic conductivity measurements to HYPROP+WP4C affects the numerical simulation results.

Note that with the present method, the SHP in still unsaturated soil close to full pore water saturation reflecting the effect of the macropores could not be identified with the HYPROP method and relatively small sample size. As reported by Yeh et al. [78], under the same saturated hydraulic conductivity function, the wetting front of the bimodal model moves down faster. This results in changes in the pressure head, water content, and internal stress of the soil. The results show that the water content and suction stress changes in the bimodal model are higher than those of the unimodal model due to the difference in water retention capacity. The present parameter determination and flow modeling still assume local equilibrium conditions in the pressure head of the pore water in the different pore systems. In the case of non-equilibrium-type preferential flow, when infiltrating water bypasses most of the unsaturated soil matrix [32], the kinetics of the capillary-driven macropore–matrix exchange flux needs to be considered, in addition to the bi- or multimodal porosity functions. In this case, the flow modeling requires using a dual-permeability approach with two interacting flow domains, in which the two-pore systems are treated separately, each having separate hydraulic functions. The validity of the increasingly complex modeling approaches, however, needs to be further tested by comparing with field or experimental tracer data, e.g., on soil columns.

4. Conclusions

The study illustrated that erosion-affected soil structural properties governing hillslope hydrology in the arable landscape, in this case, vineyards, were evident and had a significant impact on SHP and, consequently, soil water dynamics. Both unimodal and bimodal soil hydraulic models fitted the data agreeably; although, it can be clearly noticed that the bimodal model performed better in particular cases where data showed non-uniform pore size distributions. HYDRUS-1D simulations showed, in general, that both models provided a similar distribution of flux components between infiltration, surface runoff, and drainage (bottom flux) in most cases. Overall, the differences generated when using the bimodal hydraulic functions can lead to a large discrepancy in water flow quantification. Location I bottom profile simulation revealed the influence of the shape of two functions (SWRC, and SHCC), which control the relation between soil water capacity and hydraulic conductivities in the case of bimodal function during drainage, leading to increased bottom flux. It is evident that the SHP and water dynamics in highly erosion-affected heterogeneous soils with developed structure and pore space (e.g., compacted soil with cracks and biopores) cannot be adequately explained using the unimodal porosity functions or by applying single porosity models. However, the validity of more complex approaches should be further tested and parametrization should be performed with extra care, as using the non-appropriate model can lead to errors in the water balance. The study indicates that local scale SHP were of more considerable importance and had more relevance as compared to fluxes at the hillslope scale, and that an accurate parameterization of SHP is required to capture the hydropedological response of erosion-affected soil systems.