Determination of Soil Hydraulic Parameters and Evaluation of Water Dynamics and Nitrate Leaching in the Unsaturated Layered Zone: A Modeling Case Study in Central Croatia

Abstract

Nitrate leaching through soil layers to groundwater may cause significant degradation of natural resources. The aims of this study were: (i) to estimate soil hydraulic properties (SHPs) of the similar soil type with same management on various locations; (ii) to determine annual water dynamics; and (iii) to estimate the impact of subsoil horizon properties on nitrate leaching. The final goal was to compare the influence of different SHPs and layering on water dynamics and nitrate leaching. The study was conducted in central Croatia (Zagreb), at four locations on Calcaric Phaeozem, Calcaric Regosol, and Calcaric Fluvic Phaeozem soil types. Soil hydraulic parameters were estimated using the HYPROP system and HYPROP-FIT software. Water dynamics and nitrate leaching were evaluated using HYDRUS 2D/3D during a period of 365 days. The amount of water in the soil under saturated conditions varied from 0.422 to 0.535 cm³ cm⁻³ while the hydraulic conductivity varied from 3 cm day⁻¹ to 990.9 cm day⁻¹. Even though all locations have the same land use and climatic conditions with similar physical properties, hydraulic parameters varied substantially. The amount and velocity of transported nitrate (HYDRUS 2D/3D) were affected by reduced hydraulic conductivity of the subsoil as nitrates are primarily transported via advective flux. Despite the large differences in SHPs of the topsoil layers, the deeper soil layers, having similar SHPs, imposed a buffering effect preventing faster nitrate downward transport. This contributed to a very similar distribution of nitrates through the soil profile at the end of simulation period. This case study indicated the importance of carefully selecting relevant parameters in multilayered soil systems when evaluating groundwater pollution risk.

1. Introduction

Determination of soil hydraulic properties (SHPs) is very important for various disciplines such as hydrology, ecology, environmental sciences, soil science [1], agricultural or groundwater management [2], and for modeling water and solute transport in unsaturated soils [3,4,5,6]. Soil hydraulic properties describe macroscopic relations between the phase concentration, chemical potential, and transmission behavior of water, solute, and gas movements in the soil [1]. These relations depend on many factors including temperature, matrix surface characteristics, chemical composition of soil solution, pore geometry in soil, etc. SHPs are the basis for understanding the flow and transport processes, and have an impact on water retention, rate of water flow, the movement of the nutrients, chemicals, and pollutants in soil, but also determine the accessibility of plant for water uptake, crop growth, and environmental quality [7,8]. Understanding of these processes is essential for effective soil and water management [7]. For example, one study [9] focused on the effects of soil conditioners on hydraulic properties, leaching processes, and denitrification in a silty-clay soil. It was performed in a laboratory, using a plexiglass soil columns, with the aim to determine the amount of leached metals and nutrients after simulated stormwater events. Results showed decreased nitrate leaching into deeper soil layers by increasing denitrification. Another study [10] was conducted to track the effects of biochar on the hydraulic and physical properties of compacted soils. The authors showed that decreased infiltration rate and K_s after biochar application can have a positive impact on storing water in soil pores and reduce the possible leaching, and, thus, the environmental pollution. Hydraulic properties of the topsoil layer are also affected by land use and management, which indicates that land use should be taken into consideration when studying nutrient leaching in soil, especially on arable areas [11]. Since the determination of SHPs could be time-consuming and expensive, with undisturbed soil samples needed, pedotransfer functions (PTFs) [12] and remote sensing, e.g., [13,14,15] are often used as an alternative to estimations of soil hydraulic properties. However, the disadvantage of PTFs compared to SHPs is that soil properties should not be generalized and all factors affecting SHPs should be considered, e.g., [16,17].

The solute transport in the soil vadose zone is one of the most complex and demanding issues in numerical modeling, especially important for the estimations of groundwater pollution. Groundwater pollution has a rising trend worldwide and nitrate (NO₃⁻) is considered the most widespread pollutant [18]. Nitrogen fertilizers are commonly applied in agriculture to increase crop production [19,20], and its movement in the soil depends on soil type [21], soil organic matter and humus content in soils [22,23], temperature and precipitation [24], irrigation [25], and tillage [26,27]. The main issue with nitrates is their mobility in soil [28], due to low adsorption to the soil particles [29], but also the fact that they can persist in shallow groundwater for years [28]. Nitrogen is leaching when transformed in nitrate form due to its anionic nature [30], causing the potential risk for ecosystems (e.g., eutrophication) [31,32]. As a result of nitrate leaching in ground and drinking water, several health issues may occur, such as methemoglobinemia in infants or even stomach cancer in adults [33].

Water flow and solute transport models are used to describe and predict specific processes in unsaturated and saturated soil zone. Different models can be used to test and implement experiments with different soil types or crops, including the application of fertilizers and pesticides. The use of field and laboratory experiments combined with numerical models could prove to be an appropriate approach for the quantification of agrochemical transport through the soil profile in order to prevent their leaching into the groundwater. Nowadays, models can be used to predict various processes in the ecosystem. For example, the HYDRUS model is widely used for estimations of vadose zone processes [5]. It can be used for salt leaching evaluation [34], water and potassium movement simulations [35], for tracking nitrogen dynamics in unsaturated and saturated soil [36], and for modeling water flow and fertilization scenarios [37]. HYDRUS-2D was also used for modeling water and nitrate transport in a potato field [38]. The assumption was that a significant amount of nitrate can be leached below the root zone because the texture of the soil is coarse, and the potato plant has a shallow root system. According to this study, even a small increase in irrigation emitter discharge showed a significant increase in water percolation and nitrate leaching. Furthermore, another study [39] was performed to assess water flow and nitrate fluxes using HYDRUS-2D, focusing on the efficiency of zero tension plate lysimeters, which were implanted in silty-clay soils affected by a high groundwater table. The authors showed that when the groundwater table was low, water and solute diverged from the plate towards the drier surrounding soil. Using HYDRUS-2D, soil nitrate dynamics were evaluated in an intercropping dripped ecosystem [40]. The study was conducted in the corn and tomato intercropping ecosystem using different N-fertilizer applications. Their results showed that nitrate (NO₃-N) leaching was higher in the corn region than in a tomato region (i.e., more N-fertilizer was applied in the corn region), there was intensive NO₃-N exchange in the horizontal direction between different crop regions because NO₃-N moved from the tomato to corn region, and, correspondingly, NO₃-N leaching in tomato region was low. The soil NO₃-N storage in the root zone of corn was lower than in the tomato region, which indicated a large effect of agricultural management on N distribution. In a conducted study, HYDRUS-2D was also used to estimate the effect of plastic mulch cover on water and nitrate dynamics [41]. The results of numerical modeling showed large influence of plastic mulch cover on water and nutrient distribution in soil. The nitrate transport under mulch was slower compared to system without soil cover. The combination of using field experiments complemented with numerical modeling on soil water and nitrate distribution have proven to give interesting and relevant new insights on nitrogen dynamics in agroecosystems.

Here, we expanded field- and laboratory-derived soil data with the numerical modeling using HYDRUS 2D/3D software package. The main goal was to assess the behavior of water dynamics and nitrate leaching in a layered soil profiles within the similar soil type and agroecosystem management. The aims of this study were (i) to estimate soil hydraulic properties of the similar soil types at various locations; (ii) to determine annual water dynamics at different locations; and (iii) to estimate the impact of subsoil horizon properties on nitrate leaching.

2. Materials and Methods

2.1. Study Locations and Soil Properties

The study was conducted in the near proximity of the City of Zagreb, Croatia, at four locations where soil profiles (P−1 to P−4) were exposed (Figure 1; P−1–45°48′21.9″ N 16°02′37.4″ E; P−2–45°45′59.6″ N 16°02′59.9″ E; P−3–45°45′09.9″ N 15°58′57.3″ E; P−4–45°44′03.7″ N 15°58′08.3″ E). The study locations were chosen because of the same land use (currently unused plot left as meadow, but bordering agricultural soil under crop) and similar climatic conditions.

Disturbed and undisturbed soil samples were taken from two soil depths at each location (surface or topsoil layer and subsoil). Disturbed soil samples were used for texture, carbonate content (CaCO₃), porosity, bulk density, and organic carbon (C_org) determination. These samples were taken at P−1 location from 0–25 cm, 25–55 cm, and 55–100 cm depth; at P−2 location from 0–30, 30–85 cm, and 85–110 cm; at P−3 location from 0–30, 30–78 cm, and 78–110 cm; and at P−4 location from 0–42, 42–82 cm, and 82–110 cm. Undisturbed soil samples (250 cm³) were taken for soil hydraulic parameters (SHPs) estimation in two repetitions within the first and third sampling depth of disturbed samples, i.e., at P−1 location 10–15 and 75–80 cm; at P−2 15–20 and 95–100 cm; at P−3 15–20 and 90–95 cm; and at P−4 location from 25–30 and 90–95 cm. Long-term (1949–2019) average annual temperature for this area is 11.0 °C with an average rainfall of 861.5 mm, respectively [42]. The groundwater level depends on meteorological conditions and inflows from the highland region of Vukomeričke gorice. At low water levels, the depth of groundwater ranges from 7 to 11 m, and at high water levels from 4 to 9 m [43]. Soil properties of the investigated locations are shown in

Table 1. Soil properties at investigated locations (P−1 to P−4).

Locations	Depth [cm]	Texture				Munsell Color Chart	CaCO3 Content [%]	Porosity [% vol.]	Bulk Density [g cm−3]	Corg [g kg−1]
		Coarse Sand	Fine Sand	Silt	Clay
P−1	0–25	1.9	6.3	74.3	17.5	2.5Y5/3	-	52.2	1.22	24.01
	25–55	0.8	7.7	73.9	17.6	2.5Y6/3	7.3	49.3	1.33	11.43
	55–100	1.2	9.2	75.9	13.7	2.5Y6/6	5.7	47.6	1.38	8.12
P−2	0–30	0.6	2.2	77.8	19.4	2.5Y6/2	22.7	49.8	1.31	22.97
	30–85	1.6	1.1	82.8	14.5	2.5Y6/3	26.2	43.9	1.43	8.87
	85–110	1.3	6.1	81.7	10.9	2.5Y7/4	27.7	46.3	1.45	7.08
P−3	0–30	8.0	16.4	62.0	13.6	2.5Y5/3	19.9	49.4	1.33	16.36
	30–78	0.8	3.9	75.9	19.4	2.5Y7/3	22.0	44.3	1.47	10.67
	78–110	0.6	34.8	57	7.6	2.5Y6/3	25.0	49.1	1.38	4.52
P−4	0–42	1.2	15.6	67.0	16.2	2.5Y5/3	12.4	46.3	1.43	14.09
	42–82	0.3	9.5	75.8	14.4	2.5Y7/4	20.4	44.1	1.50	5.57
	82–110	0.9	3.7	81.1	14.3	2.5Y7/3	28.3	43.8	1.50	4.52

. The presence of skeletal fraction (11.6%) was determined only in the topsoil (0–30 cm) at the P−3 location. According to IUSS Working Group WRB soil classification [44], soil profile at P−1 location is Calcaric Phaeozem (Siltic) with soil horizons A-AC-C, P−2 location is Calcaric Regosol (Aric, Fluvic, Raptic, Siltic) with soil horizons Ap-C-2C, P−3 is Calcaric Fluvic Phaeozom (Aric, Raptic, Siltic) with soil horizons Ap-2C-3C, and P−4 is Calcaric Fluvic Phaeozem (Aric, Raptic, Siltic) with soil horizons Ap-2C-3C.

Bulk density was determined according to standard method [45], soil particle size distribution method by sieving and sedimentation [46], and the organic matter content by sulfochromic oxidation [47]. Soil profiles were characterized according to FAO Guidelines for soil profile description [48]. Total porosity was determined indirectly, by calculations using soil density data [49].

2.2. Soil Hydraulic Parameters Estimation

Undisturbed soil cores (250 cm³) were taken in two repetitions from two selected depths. In P−3 profile (the topsoil layer), the measurement was performed only once due to the large pores inside the cylinder (earthworm hole—the sample was excluded). Soil samples were taken from two selected horizons—top layer (0–25 cm, 0–30 cm, 0–30 cm, and 0–42 cm at P1 to P4, respectively), and subsoil (55–100 cm, 85–110 cm, 78–110 cm, and 82–110 cm at P1 to P4, respectively).

Soil hydraulic parameters (SHPs) were estimated on the undisturbed soil cores (250 cm³) using the simplified evaporation method [50] and the HYPROP automatized system which are applicable to most soil types [51]. Two tensiometers were placed inside undisturbed soil sample at two depths and measured the soil water tension [52]. The evaporation method considers the change of the sample weight and matrix potential in the soil sample during the evaporation drying process [52] and allows an accurate characterization of the water-retention properties in the porous media [53]. Fitting quality for soil hydraulic parameters estimation is given in the terms of the root mean square error (RMSE), which indicates the mean distance between data point and the fitted function [54]. R² is a coefficient of determination used as a fitting parameter for soil hydraulic parameters as well. The same method was used on a large number of undisturbed soil samples of different textures and indicated its reliability [50].

SHPs were described using the van Genuchten–Mualem model (VGM) [55]:

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

(1)

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

(2)

$$S_e=\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}$$

(3)

$$m=1-\frac{1}{n} n > 1$$

(4)

where ${\theta }_r$ and ${\theta }_s$ denote residual and saturated volumetric water contents (L³ L⁻³), respectively, h is the pressure head (L), $S_e$ is the effective saturation (–), α (L⁻¹) and n (–) are shape parameters, and $l$(–) is a pore connectivity parameter. Pore connectivity parameter, l, was fixed to a value of 0.5, as recommended for most soils [56].

2.3. Governing Flow and Transport Equations

Numerical modeling was performed using the HYDRUS 2D/3D which uses Richards’ equation to solve water flow and the advection-dispersion equation to solve the nitrate leaching [5]. The simulation of water flow in a two-dimensional plane is described using Richards’ equation for Darcy’s water flow in (un)saturated porous medium:

$$\frac{\partial \theta \left(h\right)}{\partial t}=\frac{\partial }{\partial x_i}\left[K\left(h\right)\left(K_{ij}^A\frac{\partial h}{\partial x_j}+K_{iz}^A\right)\right]-S\left(h\right)$$

(5)

where $\theta$ is volumetric water content (L³ L⁻³), h is soil water potential (L), x_i are the spatial coordinates (i = 1,2) (L), Z is vertical coordinate (L), K is unsaturated hydraulic conductivity (L T⁻¹), $K_{ij}^A$ is anisotropy tensor (−), t is time (T), and S is the water which the plant absorbs (T⁻¹). Two-dimensional modeling was used in order to visually present water and nitrate distribution with depth at selected dates which is only available in the 2D/3D version and profiles can be directly compared.

In a general form, the equation for solving water flow of (un)saturated medium can be expressed as:

$$\frac{\partial \theta }{\partial t}=\nabla \left(K\nabla H\right)-S_w$$

(6)

where $\theta$ is volumetric water content (L³ L⁻³), K is unsaturated hydraulic conductivity (L T⁻¹), H is hydraulic head (L), $S_w$ is the water which the plant absorbs (T⁻¹), t is time (T), and $\nabla$ is spatial gradient operator.

Potential evapotranspiration is calculated in HYDRUS using Penman-Monteith combination Equation (7) for reference evapotranspiration [57]. Reference evapotranspiration (ET₀) is defined as the rate of evapotranspiration from a hypothetic crop height of 12 cm, an albedo of 0.23, and a fixed canopy resistance of 70 sm⁻¹, closely resembling evapotranspiration from an extensive surface of green grass of uniform height, completely shading the ground, not short of water and actively growing.

$$E{T}_0=\frac{0.408\mathsf{\Delta }\left(R_n-G\right)+\gamma \frac{900}{T+273}{U}_2\left(e_a-e_d\right)}{\mathsf{\Delta }+\gamma \left(1+0.34U_2\right)}$$

(7)

where $E{T}_0$ is the reference crop evapotranspiration (mm day⁻¹), $R_n$ is net radiation at the crop surface (MJ m⁻² day⁻¹), $G$ is the soil heat flux (MJ m⁻² day⁻¹), $T$ is the average air temperature (°C), $U_2$ is the windspeed measured at 2 m height (m s⁻¹), $\left(e_a-e_d\right)$ is the vapor pressure at dew point (kPa), $e_a$ is the saturation vapor pressure at temperature T (kPa), $e_d$ is the vapor pressure at dew point (kPa), and 900 is a conversion factor.

The advection dispersion equation was used to model nitrate transport under the assumption that this is a non-reactive transport of substances:

$$\frac{\partial \theta c}{\partial t}=\nabla \left(\theta D\nabla c\right)-\nabla \left(qc\right)$$

(8)

where c is concentration of the substance in the liquid phase (M L⁻³), D is the dispersion coefficient tensor (L² T⁻¹), and q is the volumetric flux density (L T⁻¹).

The HYDRUS 2D/3D model uses soil hydraulic functions [55,56] to describe soil water retention curves—$\theta \left(h\right)$ and hydraulic permeability of unsaturated soil—$K\left(\theta \right)$.

2.4. Numerical Modeling Using HYDRUS 2D/3D: Initial and Boundary Conditions

HYDRUS is applicable for the prediction of multiple soil-related processes which can be used for various purposes such as optimization of irrigation systems [58,59], estimations of the influence of plants on the amount of water in soil, modeling of groundwater oscillations [60], simulations of the transport of various substances, e.g., salts, nitrates, metals, pesticides [61], as well as to determine the impact of land use and environmental change on soil processes [62].

HYPROP device was used to determine hydraulic parameters (residual and saturated water content and hydraulic conductivity) in undisturbed soil samples measuring soil water potential to create a moisture release curve. Based on the obtained data, which were then used as an input for the HYDRUS 2D/3D, water flow and solute transport modeling was performed.

Numerical simulations were conducted for the period from 1 December 2017 (date of the field sampling) to 30 November 2018, in two-dimensional square domain 50 cm wide and 200 cm deep. The last horizon properties were used for the zone below 1 m depth as no measurements were done in the field. The goal was to assess water and nitrate dynamics within and below the root zone. The domain represents the default volume within which the flow of the observed fluid takes place. Transport parameters of the soil were adjusted to 2 m depth as commonly used in similar modeling, e.g., [63]. The boundary water flow conditions were set as atmospheric boundary conditions at the top and free drainage conditions at the bottom of the profile. The initial conditions were set at the level of field capacity (FC) of soils studied as observed during the profile excavation. The initial conditions for nitrate concentration were set to 0 mmol cm⁻³. For nitrate leaching, third-type solute transport boundary condition was applied with imposed 1 mmol cm⁻³ of solute pulse at the simulation initiation. This type is used when the solute flux along a boundary is specified and its mass conservative. This was done in order to set the same solute initial conditions at all locations, as the aim was to compare the influence of different SHPs and layering on water dynamics and nitrate leaching. For space discretization, Upstreaming Weighting FE type method with two-dimensional network was used. Time derivations were solved by Crank–Nichols scheme. Millington and Quirk tortuosity model was used to model solute transport.

3. Results

3.1. Variability of Soil Hydraulic Parameters

Water retention and hydraulic conductivity curves of samples were determined (Figure 2 and Figure 3) using HYPROP device. Soil hydraulic properties (SHPs) for the investigated locations were determined using the HYPROP-FIT program, as shown in

Table 2. van Genuchten–Mualem (VGM) hydraulic parameters measured by HYPROP-FIT for locations P−1 to P−4. The θ represents retention and K hydraulic conductivity. RMSE and R² indicate the VGM model fit with the measured data.

Profile	Sampling Depth [cm]	Depth [cm]	θr [cm3 cm−3]	θs [cm3 cm−3]	α [1 cm−1]	n	Ks [cm day−1]	RMSE_θ	RMSE_K	R2_θ	R2_K
P−1	10−15	0−25	0.000	0.535	0.10310	1.200	909.2	0.0100	0.3246	0.9301	0.9369
	10−15	0−25
	75−80	55−100	0.000	0.486	0.01150	1.278	17.2	0.0073	0.1320	0.8777	0.8665
	75−80	55−100
P−2	15−20	0−30	0.000	0.459	0.00959	1.200	5.8	0.0088	0.3562	0.8561	0.8883
	15−20	0−30
	95−100	85−110	0.000	0.471	0.00778	1.288	6.2	0.0042	0.0808	0.8663	0.8531
	95−100	85−110
P−3	15−20	0−30	0.000	0.499	0.08460	1.200	990.9	0.0033	0.1877	0.9220	0.9367
	90−95	78−110	0.123	0.452	0.01350	1.736	13.8	0.0034	0.1386	0.8469	0.8392
	90−95	78−110
P−4	25−30	0−42	0.000	0.428	0.47400	1.200	113.3	0.0066	0.3957	0.9098	0.9232
	25−30	0−42
	90−95	82−110	0.000	0.422	0.00218	1.210	3.0	0.0087	0.2609	0.8137	0.8245
	90−95	82−110

. In previous studies, using HYPROP resulted in an equal or even more accurate estimation of the soil retention curve compared to the scarce data obtained by the traditional approach (sandbox and pressure extractor) [2]. The amount of water in the soil under fully saturated conditions (θ_s) in the study area (Table 2) varied from 0.422 (profile P−4) to 0.535 cm³ cm⁻³ (profile P−1). The variability of θ_s depended on the physical and chemical properties of the soil, and their coupling resulted in different water retention in the soil at different pressures. In our study, the repetitions used for the estimations of soil water retention curves (SWRC) showed similar values, which indicated the reliability of applied method. Figure 2 shows retention curves for locations P−1, P−2, P−3, and P−4, and Figure 3 shows the ratio of the relative amount of water in the soil and the hydraulic conductivity which are measured in the undisturbed soil samples. Blue and red circles represent the observation point of two undisturbed samples. In Figure 3, subsoil at the P−4 location shows errors at the beginning of the measurement (red circles), which can be a result of low K_s values thus affecting pressure differences. In the second set of measurements (blue circles), these errors were smaller and the measurement was more precise.

The hydraulic conductivity (K_s) varied from 3 cm day⁻¹ (profile P−4) up to 990.9 cm day⁻¹ (profile P−3; Table 2). The lowest K_s was measured in the subsoil of P−4 profile compared to other profiles (P−1 to P−3). Lower total porosity (43.8%) and bulk density of 1.5 g cm⁻³ (the highest measured) caused a very low hydraulic conductivity. The highest value of K_s was measured in the topsoil layer of P−3 profile. Compared to the other locations, this profile has the highest amount of sand particles (24.4%), lower bulk density (1.33 g cm⁻³), and a total porosity of 49.4%. Total porosity was similar at all locations and for all horizons, except the subsoil at the P−4 location, but hydraulic conductivity had a wide value range depending on the location. Bulk density increased (Table 1) with soil depth at all locations, while hydraulic conductivity decreased with soil depth with the exception of the P−2 location (Table 2).

In the topsoil layer, R² for measuring water retention and hydraulic conductivity was greater than 0.90 at the P−1, P−3, and P−4 locations (Table 2). In the subsoil of P−1, P−3, P−4, and whole soil profile at the P−2 location, R² was between 0.81 and 0.89, which also indicated high model efficiency. RMSE_θ < 0.01 and RMSE_K < 0.4 (Table 2) indicated the applicability of van Genuchten–Mualem model for all locations. It also indicated that the applied method was reliable for the assessment of soil hydraulic parameters [2].

3.2. Water Flow and Nitrate Transport Modeling

Using the HYDRUS 2D/3D program, the soil water content and nitrate transport were simulated at the investigated locations over a period of 365 days. Soil water content varied during the study period from 0.140 to 0.440 cm³ cm⁻³ (Figure 4), depending on the intensity, layout, and distribution of precipitation in the area (Zagreb). According to simulations presented in Figure 4 for the 60th and 365th day, the highest value of the soil water content was recorded at location P−2 and it was between 0.380 and 0.440 cm³ cm⁻³ for the two presented days. Soil water content was higher during other periods over the year because θ_s for this location was 0.459 in the topsoil layer and 0.471 cm³ cm⁻³ for the subsoil. The maximum water retention of 0.440 cm³ cm⁻³ inside the P−2 soil profile was recorded on day 365. In the topsoil layer at location P−4, soil water content had the lowest values, and water content in the subsoil was higher. The P−3 location had a moderate water content.

According to obtained results from HYPROP and HYPROP-FIT and especially the wide range of K_s values, we expected large differences in nitrate leaching at different locations and also the faster nitrate transport followed by the higher K_s values. After the soil annual water content was successfully simulated, the NO₃⁻ transport was modeled. Three days of the investigated year are presented, i.e., 60th, 121st, and 365th day (Figure 5). Nitrate concentration ranged from 0 to 5.68 × 10⁻³ mmol cm⁻³. For all three days for which the simulations are presented, the fastest nitrate transport was recorded at the P−3 location. As expected, the slowest nitrate transport for all simulated days was recorded at the P−2 location, which is influenced by the lowest K_s through the soil horizons compared to other investigated locations.

Figure 6 shows the transport of nitrate through the soil with respect to depth and climatic parameters (precipitation, evapotranspiration, and temperature) for the 60th and 365th day. At the end of simulation period, it is evident that the nitrates are transported to a soil depth of 200 cm at all locations, which indicates the possibility of groundwater pollution. It can also be perceived that the mass of nitrate on the last simulated day is lower in comparison to the start of the research. Figure 7 shows the influence of time and precipitation on water and solute fluxes for the entire simulated period. The water fluxes (Figure 7a) in the simulated soil column are in line with precipitation. Nitrate leaching throughout the soil column had a delay in comparison with precipitation (Figure 7b), showing the time needed for nitrates to reach soil profile bottom boundary. The differences in leaching rate can be seen and showed the earliest leaching at the P−3 location.

4. Discussion

4.1. Variability of Soil Hydraulic Parameters and Layering Impact

According to the obtained results, it is evident that soil water capacity depends on organic matter and clay particles content, i.e., at the location that has higher value of clay particles in the horizon, the soil water capacity was also higher, and vice versa. Organic matter, e.g., [64], and clay particles, e.g., [65], have a positive impact on water retention. Although all four locations were covered with grass, differences between locations and depth were observed. Grass is mostly rooted in the first 15 to 20 cm of soil depth, thus affecting physical and chemical soil properties. Roots, besides affecting soil hydraulic conductivity and soil water storage, also control connectivity and preferential flow pathways and impact the infiltration [66]. In general, topsoil layers mostly have better soil structure, are more aerated, have a larger amount of organic matter, e.g., [67], lower bulk density, e.g., [68], and thus higher porosity, e.g., [69], while soil structure, organic matter, and porosity are mostly depleting with depth and bulk density is increasing. Results showed that total porosity was similar at all investigated locations. However, total porosity does not give data on the presence and architecture of micro-, meso-, and macropores which can cause differences in water and solute transport through the soil, and it is known that macropores have a big influence on increased water flow in soil [70]. It is also established that as bulk density increases, the K_s decreases [71,72]. Although these physical properties can cause very fast hydraulic conductivity, these values are still quite high. Biopores, roots, preferential flow, or carbonate content may have had an impact on K_s. Biopores in the soil are present in a wide size range and can influence preferential flow. In most cases, large-sized biopores are air-filled due to the inability to retain water or lack of capillary forces. After rainfall, the water replaces the air in the pores but the soil cannot retain it and it is rapidly transported to the deeper layers. The P−3 profile subsoil also had the highest amount of sand, and the lowest values of silt and clay, but still much lower K_s values. Some of the reasons can be higher bulk density and compaction with depth or lower amount of C_org, all of which have a negative effect on soil structure and thus lower K_s. At the P−2 location, K_s was higher in the subsoil than in the topsoil layer, which may be caused by soil degradation after rainfall events [73]. At this depth, the amount of sand particles is also very low and thus the porosity is lower than of other soils and layers. P−2 had the least permeable soil with low hydraulic conductivity throughout the soil profile which can be a result of a higher content of silt and clay particles since the saturated hydraulic conductivity is negatively affected by the clay content [8], i.e., K_s increases with particle sizes [74]. The spatial variation of K_s is greater for the soil surface compared to the subsurface which may be the cause of minor differences found in the subsoil [75]. Variations of bulk density and hydraulic conductivity with depth could also be explained by the impact of soil organic matter [65]. Soil physical and chemical properties such as the content of organic matter [76], soil texture and structure [77], compaction [78], and porosity [79] have an impact on soil hydraulic conductivity. Carbonate content could also have an impact on K_s, especially in the subsoil. CaCO₃ is transported and precipitated through the soil and its values increase with depth (except at the P−1 location), thus affecting pore size and consequently soil hydraulic conductivity [80]. Although all soil profiles had the same land use for the past few years, physical and chemical properties (and thus soil hydraulic properties) were slightly different. The greatest difference was between the hydraulic conductivity values, but also, percentage of sand, silt and clay particles, porosity, bulk density, and C_org, all influencing water and solute transport. Even at the same location, there was the difference between the top- and subsoil layer (particularly pronounced for P−1 and P−3 profiles). These results show that soil properties, leaching, and transport estimation should not be generalized solely based on soil type. According to these results, when mapping the soil moisture, estimation could be very challenging because it can vary highly in spatial and temporal distribution [81] and using soil texture for pedotransfer functions might not be sufficient for accurate mapping of heterogenous soil hydraulic properties [11].

4.2. Water Dynamics and Nitrate Leaching Potential

The maximum water retention simulated inside the P−2 soil profile recorded on the last day of simulations (365th) possibly was result of high clay and C_org content, which have high adsorption capacity. However, climatic conditions (increased precipitation and low temperatures resulting in lower evapotranspiration), high content of silt, and low K_s values could have also contributed to the soil water retention. These results are in line with another study whose results showed that higher water retention and water content are results of high silt and clay content [82]. On the other hand, the lowest values of soil water content (topsoil layer at P−4) are possibly affected by the highest values of bulk density because of a negative correlation between water content and bulk density [64]. Higher water content in the subsoil at the same location could be a result of lower porosity since the compaction reduces total porosity and increases the volumetric water content [83]. Moderate water content at the P−3 location may be affected by the negative correlation between the water content and coarse sand and total sand [64].

The fastest nitrate transport recorded at the P−3 location can be explained by the highest sand content, i.e., the lowest values of silt and clay, thus resulting in higher K_s values in the topsoil compared to the other investigated locations. As mentioned above, K_s increases with particle size [74], and thus increases the possibility of nitrate leaching. In another study [84], it is also shown that soil with a coarse texture, such as sands, allows faster nitrate leaching to groundwater. In our case, groundwater was not raised to 2 m depth, but we wanted to compare various locations. If the groundwater at these locations would go to 2 or less meter depth, the risk of nitrate leaching would be even higher and the transport would possibly be even faster. On the other hand, the slowest nitrate transport was recorded at the P−2 location, which also has the lowest K_s values through the whole soil profile. The movement of nitrogen is closely related to soil water flow [40]; thus, K_s influences its transport. All investigated locations (P−1, P−2, P−3, and P−4) had an approximately uniform nitrate transport through the soil. Similar volume and velocity of leached nitrate could be a result of reduced hydraulic conductivity in the subsoil. It is evident that lower K_s values in the subsoil had an impact on the topsoil layer and regulated nitrate transport. The coefficient of the permeability of the deepest layer is a controlling factor of permeability of the whole profile [85]. This could consequently affect the nitrate transport through the soil. The presence of distinguished layers in the soil also controls the infiltration of water and its redistribution throughout the soil [66]. According to this study, it is evident that soils within the same climate region can vary in water flow and nitrate leaching. Thus, at these locations, more soil samples should be taken for additional accuracy. It is important to take more samples per depth and also throughout the profile for soil hydraulic properties and transport estimation. By taking samples from more than one layer, it is possible to see differences within depth and it is easier to track how soil hydraulic properties change throughout the soil depth and thus causing changes in flow direction while increasing or decreasing solute transport. On the last day of the study (365th), the model suggests that the nitrate would be leached below 200 cm of soil depth at all locations.

In the simulations of nitrate transport in relation to soil depth, it is evident that the transport rate is defined by the hydraulic characteristics of the soil, primarily by the hydraulic conductivity. During the investigated year, the amount of precipitation was higher than the multi-year average, which possibly increased nitrate leaching into deeper soil layers. The issue suggested by these results could be especially significant during the winter months when there is a lot of precipitation and prolonged rains often occur, although high intensity rainfall events common during the summer months could also contribute to nitrate leaching. After these rainfall events, soil pores can rapidly transport water and solute, and thus increased nitrate leaching is expected. The obtained results in a study conducted to estimate the effect of temperature and precipitation on nitrate leaching showed a significant impact of distribution and amount of precipitation on nitrate leaching [24]. Further, nitrate leaching losses were increased during the autumn–winter periods when evapotranspiration was lower [19]. The mass of nitrate is lower at the end of the conducted simulations (Figure 6), which is consistent with results obtained in a similar study [86]. Their results showed that the content of total nitrogen and nitrate nitrogen (NO₃-N) is decreased with the increase of soil depth.

The amount and velocity of transported nitrate (HYDRUS 2D/3D) were affected by reduced hydraulic conductivity of the subsoil, e.g., [75], as nitrates are primarily transported via advective flux. Despite the large difference in SHPs of the topsoil layer, the deeper soil layers with similar SHPs imposed a buffering effect and prevented faster nitrate downward transport. This contributed to very uniform distribution of nitrates through the soil profile at the end of simulation period. As shown in Figure 8, in conceptual modeling example, subsoil hydraulic parameters directly influence nitrate percolation to the deeper soil layer, e.g., [85], and possibly groundwater quality. In the modeling example, NO₃⁻ was simulated for 30 days for soil profiles with (i) increased topsoil layer porosity (K_s = 168 cm day⁻¹, α = 0.1 and n = 1.8 VGM properties values), and (ii) uniform (low) porosity profile (K_s = 1.68 cm day⁻¹, α = 0.01 and n = 1.2 VGM properties values). The simulated example indicated the importance of correctly estimating SHPs throughout the entire soil depth when nitrate leaching is studied. The importance of subsoil parameter proper estimation is also evident, as often the deeper soil layer properties are not evaluated, or data is missing. Neglecting the subsoil parameters and relying only on topsoil layer hydrology could lead to misrepresentation of nitrate transport to groundwater, which was also alerted in several studies [16,17].

5. Conclusions

The methods used in this study provide a basis for assessing the risk of agrochemicals water pollution, and for the possible application of predictive numerical models in research of agrochemicals transport in agroecosystems. Soil hydraulic parameters were estimated using the HYPROP system and evaluated in the HYPROP-FIT program. Hydraulic parameters of the soil at investigated locations under same land use (meadow, i.e., currently unused plot bordering agricultural soil under crop) showed that the least permeable soil was at the P−2 location, and that the highest value of K_s was measured in topsoil at the P−3 location, both resulting from the site-specific soil properties. The soil water content varied during the study period and ranged from 0.140 to 0.440 cm³ cm⁻³ in the presented simulation snapshots, depending on the intensity, layout, and distribution of precipitation. As one-year simulations showed, all nitrates present in soil will be leached below 200 cm of soil depth, which can cause groundwater pollution. Although climatic conditions and land use were the same at all locations, their physical and chemical properties, and thus the soil water content, as well as the possibility of nitrate leaching differed. A wide range of hydraulic conductivity values was noticed in the topsoil layers, but these differences were reduced in the subsoil which acted as buffering layer. Simulations resulted in a similar amount and velocity of transported nitrates at all four locations. It is evident from the results how many differences could be present on soil hydrology even in the same soil profile, which also affects different groundwater pollution risks. Therefore, more soil samples throughout soil depths and more repetitions should be taken so that the complete scenario can be modeled. This case study indicated the importance of carefully selecting the relevant information in multilayered soil systems when evaluating groundwater pollution risk, especially in highly diverse agroecosystems.