Predicting Water Distribution and Optimizing Irrigation Scheduling in Turfgrass Management Using HYDRUS-2D

1. Introduction

Turfgrass areas provide a multitude of ecosystem services, social benefits, and economic advantages, particularly in urban areas [1,2]. These include carbon sequestration, which is crucial in mitigating climate change. Studies have shown that turfgrass systems can sequester significant amounts of carbon with specific management practices that enhance this capability [3,4,5]. For instance, managing irrigation and fertilization regimes can influence carbon sequestration rates in turfgrass, highlighting the importance of sustainable practices to maximize these benefits [4,5]. In addition to carbon sequestration, turfgrass contributes to various ecosystem services such as heat dissipation and water infiltration [6,7,8]. These services are particularly vital in urban settings, where green spaces help mitigate the urban heat island effect and improve overall environmental quality. The presence of turfgrass in urban landscapes enhances aesthetic value and supports recreational activities, thereby promoting social well-being [9,10].

Nevertheless, turfgrass requires the provision of conditions that support year-round sports activities, promote plant growth, and maintain a visually appealing appearance. Also, a good quality turfgrass cover needs more water compared to most other uses, necessitating a management system that minimizes water consumption.

The task at hand is to optimize the functionality of these systems while simultaneously guaranteeing their sustainability, a challenge that can be difficult to attain since improving one aspect might negatively impact the other [11]. As urban areas and their managed turfgrass spaces continue to expand, the environmental concerns surrounding them are expected to intensify. Emphasis on sustainability and ecological maintenance management has become increasingly important due to the growing pressure to address climatic change, particularly in preparing for impending water shortages [12,13,14,15]. It is crucial to integrate diverse aspects of playability, esthetic requirements, efficiency, and sustainability into a comprehensive and compatible compromise [16,17]. The need for irrigation in regions with insufficient rainfall to maintain healthy and aesthetically pleasing turfgrass has been a topic of significant discussion [18]. The efficiency of turfgrass irrigation is influenced by various climatic, technical, vegetation-related, and soil-physical factors. The rootzone construction method and irrigation system used are vital for achieving sustainable turfgrass areas [19,20].

To improve the efficiency of turfgrass irrigation, it is crucial to comprehensively assess the key factors and optimize the utilization of water in turfgrass areas to maximize irrigation efficiency. An investigation into the distribution of irrigation water under various construction methods and associated irrigation systems (sprinkler or subsurface drip irrigation) is crucial for predicting the success of the system’s design. This is because it significantly affects both irrigation scheduling and water efficiency [21,22].

The selection of appropriate rootzone construction methods with the central function of water retention (making it available to plants) and water transport (draining off excess water to maintain durability and playability) is a key factor for ensuring efficient water usage and the long-term sustainability of turfgrass areas. Sand-dominated rootzones are widely used for high-quality sports facilities because of their good drainage properties, greater air-filled pore space after compaction, and more consistent playing characteristics [23]. To increase the water-holding capacity, rootzones are amended with peat, and soil surfactants enhance moisture retention in turfgrass rootzones [24]. Common rootzone construction methods for turfgrass areas are multilayered systems, and many elementary and complex soil physical relationships can be identified in the soil matrix. Water retention and movement issues primarily influence decisions regarding sports turf construction methods [25]. Irrigation systems with high distribution accuracy and efficient water application can provide the basis for vital turf growth and sustainable use of water resources [26]. An efficient irrigation system should minimize the losses from wind drift, surface runoff, percolation, and evaporation [27]. Sprinkler irrigation systems are standard for turfgrass irrigation. These systems, using rotary, multi-jet, and spray sprinklers, are the standard for turf sports fields. Still, climatic factors such as wind drift and technically induced distribution inaccuracies may cause low irrigation efficiency, intensive percolation, and nutrient leaching in the soil matrix, which is weak in sorption and water retention [28,29].

In contrast, subsurface drip irrigation (SDI) systems are characterized by direct water delivery into the plant root zone, which can increase the irrigation system´s efficiency [30]. Subsurface drip irrigation systems typically use either point or linear sources for water supply. Linear-source systems (“porous pipes”) feature uniform water flow along the entire length of the line, whereas point-source systems, commonly referred to as drip lines, are equipped with non-pressure-compensating or pressure-compensating drippers. These SDI systems employ polyethylene pipes placed at various depths and distances and are customized to the specific plants being irrigated [31].

The rootzone construction method, along with the choice of irrigation system (sprinkler, SDI, or hybrid approach), constitutes a highly intricate design in terms of irrigation water distribution [32]. However, numerical simulation models such as HYDRUS-2D can be used to predict soil water distribution for different construction methods, irrigation systems, and irrigation strategies quickly and efficiently [33,34] and have proven to be particularly reliable and versatile [35]. Both in research into subsurface irrigation in agriculture [36,37,38] and in other ecosystems (e.g., urban greens) [39], reliable agreement between HYDRUS simulations and the actual water distribution in field trials in various soils has been demonstrated [37,40,41]. The basis of modeling water behavior in unsaturated porous media lies in describing the media’s pore structure, as evidenced by the water retention curve and the unsaturated hydraulic conductivity function. These functions are indispensable for simulating variations in water content and flux [42].

A well-known parametric model that connects volumetric water content to matric potential was introduced by van Genuchten [43], which incorporates parameters such as ϴr, ϴs, α, n, and m. This model has been extensively used in various studies [44,45] and is integrated into the HYDRUS-2D model as follows:

$${\theta }_{\left(\psi \right)}={\theta }_r+{\displaystyle \frac{{\theta }_{\mathrm{s}}-{\theta }_{\mathrm{r}}}{{\left(1+{\left|\alpha ·\psi \right|}^n\right)}^m}}$$

where ${\theta }_{\psi }$ is the water content at matric potential $\psi$; ϴ_r is the residual water content (cm³ cm⁻³); ϴ_s is the saturated soil water content (cm³ cm⁻³); α, n, and m describe the shape of the function without physical meaning; and m is usually fixed as 1 − 1/n [36].

The previously discussed formulation can be integrated with Mualem’s equation [46] to elucidate the unsaturated hydraulic conductivity function, and it is also incorporated in the HYDRUS-2D model:

$$K_{\psi }=K_{\mathrm{s}}·S_e^{ L}·{\left[1-{\left(1-S_e^{ 1/m}\right)}^m\right]}^2$$

with

$$Se={\displaystyle \frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}}$$

where ${\mathrm{K}}_{\psi }$ is the hydraulic conductivity at matric potential $\psi$; K_s is the saturated hydraulic conductivity; S_e is the effective water content; L is a parameter describing the pore structure of the soil, usually set to 0.5; and m is fixed as m = 1 − 1/n [33,47,48].

Rootzone materials, characterized by a high content of coarse pores, often exhibit poor re-wettability after drying due to air entrapment or issues with the wetting of organic materials. This phenomenon, known as hysteresis of the water retention curve [47,49], has been demonstrated to be particularly significant under SDI in several studies [36,48,50]. Consequently, when rewetting occurs at a specific matric potential, the material's water content will not achieve the same level as measured during the drying process. Thus, it is imperative to incorporate the effects of hysteresis into the retention function of the HYDRUS model to attain more precise model results, which requires that both the main drying and main wetting curves are known [51,52,53].

The HYDRUS code incorporates the description of the main drying curve through the soil hydraulic parameters ϴ_r, α, n, and ϴ_s. Differences between measured and calibrated values are often due to not considering hysteresis in the wetting curve [45]. Based on the study of Kool and Parker [54] for hysteretic description in the wetting curve, the soil hydraulic parameter α_w was employed. Their research indicated that α was the sole parameter to change for considering hysteresis, and they determined that the range of α_w/α is sufficiently narrow that assuming α_w/α ≈ 2 should provide a useful approximation in cases where data are lacking (where the subscript w indicates wetting). Additionally, the soil hydraulic parameter ϴ_s^w was implemented based on the study of Simunek et al. [33], who posited that α_w and ϴ_s^w are the only independent parameters necessary to describe hysteresis in the retention curve. Huang et al. [55] used hydraulic parameters α_w and n_w as variables to describe hysteresis and found in their comparison with the Kool and Parker model that their formulation yielded higher accuracy.

Physically based models, such as HYDRUS-2D, often require only minimal calibration when all the required input parameters (e.g., hydraulic parameters for water flow in soil) have already been independently determined [56]. Calibrating a model involves adjusting the input parameters of the odel to fall within reasonable ranges until the simulated results closely match the measured data [57]. In cases where hysteresis is a relevant factor, the parameters describing hysteresis may also need calibration.

Understanding the sensitivity of the soil hydraulic parameters helps to decide which parameters may be used for calibration (i.e. the most sensitive ones). In a study [58,59], which tested the sensitivity of soil hydraulic parameters in the HYDRUS-2D model, it was revealed that parameter n was the most sensitive parameter for the simulation output (soil water content), followed by the saturated soil water content ϴ_s while K_s was the least sensitive parameter.

According to Ghazouani et al. [60], HYDRUS tended to overestimate the modeled water contents when compared to the measured ones, which could be attributed to an imperfect parameterization of the soil hydraulic functions. Nevertheless, when accurately described and parameterized, HYDRUS demonstrated its ability to provide reliable values [61].

Widely used measures for evaluating the model quality are based on statistical indicators, namely the root mean square error (RMSE), mean absolute error (MAE), correlation coefficient (R²), and Nash-Sutcliffe efficiency (NSE), which indicate acceptable quality at values > 0.5 [36,61].

Although HYDRUS-2D is a widely used numerical model for predicting water, heat, and solute movement and distribution, no work has been carried out to verify the accuracy of this model for simulating irrigation water distribution under sprinkler and SDI systems across various multilayered rootzone construction methods for turfgrass areas.

This study aimed to compare experimentally measured greenhouse data with predicted data on irrigation water distribution. The objectives were i) to verify and enhance the accuracy of the HYDRUS-2D model, ultimately providing an affordable and rapid tool for evaluating efficient irrigation scheduling in turfgrass areas, ii) to determine the sensitivity and the necessity of a precise parameterization of the soil hydraulic parameters to achieve a high model accuracy, iii) to use the calibrated model to determine if dividing the total irrigation quantity into multiple smaller applications will enhance irrigation efficiency, and iv) to decide based on the calibrated model if a hybrid-irrigation approach (sprinkler-irrigation plus SDI) inhances irrigation efficiency compared to sprinkler or SDI alone.

2. Materials and Methods

This study used three different rootzone construction methods (2A, 2B, and 3) and two different irrigation systems (sprinkler and subsurface drip irrigation) to verify and enhance the model quality and evaluate efficient irrigation scheduling in turfgrass areas.

2.1. Experimental Setup and Measurements

Experimentally measured data were collected under greenhouse conditions using bare soil profiles (plots with no grass cover).

The research area (11.94 m × 4.71 m including non-considered edge areas of the actual test plots) was designed as a completely randomized two-factorial split plot with three replications for each treatment. The total area comprised six main plots (each 4.11 m × 1.70 m; three for SPR and three for SDI irrigation). Each of the six main plots was divided into three plots (1.70 m × 1.37 m) for the three construction types (two 2-layered and one 3-layered design). Two different irrigation cycles with discharge amounts of cycle 1 = 10 mm (the usual irrigation amount for turfgrass) and 2 = 20 mm (data cycle 2 not shown) were applied, both with an intensity of 10 mm h⁻¹. After irrigation cycle 1, all plots underwent a drying phase (approximately four weeks with an average evaporation rate of 3.22 mm day⁻¹) to achieve a uniform initial soil volumetric water content. The initial volumetric water content was between 9 and 10 vol.% without significant differences between the plots.

A total of five materials—STD, UHFC (rootzone layer), URTS F (intermediate layer), URTS C, and DS (gravel drainage layers)—were used for the three construction types. Table 1 lists the physical properties of each material.

STD is generally used for sports turf areas with a common requirement for functionality and sorption capability, consisting of sand, peat, and humus-rich topsoil. The UHFC (ultra-high-function) material, also consisting of sand and peat but without topsoil, has high functionality at the expense of sorption capability with a lower silt content than STD. It also exhibits higher hydraulic conductivity while maintaining a water-holding capacity (field capacity) similar to that of STD.

URTS F is an intermediate layer material consisting of medium-sized sand (grain size up to 0.5 mm) with very low silt content, with a much higher hydraulic conductivity and lower water-holding capacity (field capacity) than STD and UHFC.

DS is a gravel drainage layer material with a 0–8 mm particle size distribution, higher hydraulic conductivity, and lower silt content than the rootzone materials (STD and UHFC). URTS C is a much finer drainage material consisting of very coarse sand with a grain size of up to 2 mm, much higher hydraulic conductivity, and lower bulk density than DS. The construction layers had a total thickness of 24 cm, consisting of rootzone layers (0–12 cm), plus gravel drainage layer (12–24 cm) for the 2-layered constructions, or plus intermediate layer (12–18 cm) and drainage layer (18–24 cm) for the 3-layered constructions, and with an underlying 26 cm sand carpet base with a grain size up to 2 mm for surface drainage (Figure 1). Three rootzone construction methods were applied. Two construction types (2A and 2B) according to the German Standards for Sports Grounds [66] consisted of a 2-layered design, including a 12 cm peat-amended sandy rootzone (STD or UHFC) over a 12 cm gravel drainage layer (DS) with a total thickness of 24 cm. One construction method (3) comprised a 3-layer design based on the guidelines of the United States Golf Association (USGA, 2018) that included a 12 cm rootzone (UHFC) over a 6 cm intermediate sand layer (URTS F) and an underlying 6 cm drainage layer (URTS C).

Figure 1 provides an overview of the construction types, associated materials, and irrigation systems.

In each case, all materials were installed with a bulk density of 95% of the standard Proctor density [66] (STD = 1.55, UHFC = 1.46, URTS C = 1.60, URTS F = 1.41, DS = 1.80 g cm⁻³) and two systems used for irrigation. Due to the plot size and to maximize the water distribution uniformity, sprinkler-irrigated plots were hand-watered with a discharge rate of 10 mm h⁻¹. The subsurface drip irrigation system (SDI) used was a line-source drip system consisting of porous pipes (radius r = 0.9 cm) with a discharge rate of 3 L h⁻¹ m⁻¹ and 9.09 L h⁻¹ m⁻², respectively, resulting from a line spacing of 33 cm. The SDI installation depth was 16.5 cm, and the operating pressure was 0.2 bar.

The collected dataset for each evaluated sprinkler or SDI irrigated variant (2A, 2B, and 3) comprised 30 measurements (n=30), which were derived from a) three distinct observation depths (3, 6, and 12 cm) and b) ten specific observation times (0.00, 0.17, 0.33, 0.50, 2.00, 4.00, 8.00, 12.00, 24.00, 48.00 hours after irrigation initiation). Data were collected in two separate operations during both cycles and to ensure an identical initial volumetric water content (VWC), irrigation cycle 2 commenced following a drying phase of three weeks. Soil samples were taken with a soil sampler through lateral openings at −3, −6, and −12 cm using a 50 mm diameter metal tube, which was driven into the soil and withdrawn. The soil gravimetric water content was measured by oven-drying the soil samples at 105 °C, and the VWC was calculated by multiplying it with bulk density. For each replication and observation time, a total of 15 soil samples were collected from five measurement points (center of the plot above the SDI line and at 8.25 and 16.5 cm right and left of the center; at similar positions for the SPR variants) and at depths of 3, 6, and 12 cm. The five water content values for one depth for each sampling were averaged. After sampling, the voids were refilled with the same soil material.

2.2. Simulation Model, Initial and Boundary Conditions

HYDRUS-2D (H2D) finite element model version 5.04 [56] was used to simulate the distribution of irrigation water in bare soil profiles. The simulation was initiated by creating a model setup with the correct material layers, initial conditions, and boundary conditions within H2D that corresponded to the experimental setup.

The observed volumetric soil water content in the soil profile was taken as the initial water content for all simulated scenarios. The model´s soil surface was subjected to atmosphere boundary conditions for soil water evaporation (E_V) corresponding to greenhouse conditions with 0.013 cm h⁻¹. A free drainage boundary condition was imposed at the bottom of the soil profile. On the right- and left side of the soil profile, a no-flux boundary was used (Figure 2).

Sprinkler irrigation was scheduled in accordance with the experimental setup with an intensity of 1 cm h⁻¹. For SDI, a variable flux boundary was used around the SDI emitter (Figure 2). During irrigation, the drip pipe boundary had a constant water flux, which was obtained by dividing the emitter discharge flow rate of 3 L (h × m)⁻¹ by the surface of the drip pipe as:

$$q={\displaystyle \frac{\mathrm{Emitter}\mathrm{discharge}\mathrm{flow}\mathrm{rate}}{\mathrm{pipe}\mathrm{surface}\mathrm{area}}}=\left({\displaystyle \frac{3000{\mathrm{cm}}^3 {\left(\mathrm{h}\right)}^{-1}}{2·0.9 \mathrm{cm}·\mathsf{\pi }·100\mathrm{cm}}}\right) =5.306{\mathrm{cm}\mathrm{h}}^{-1}$$

During the no-irrigation period, the flux was kept at zero.

The simulation model domain was 50 cm deep. The high differencesbetween the soil hydraulic parameters of the different materials, as well as the significant temporal-spatial variability under SDI, necessitated a node spacing of 0.50 mm, with progressively closer spacing down to 0.25 mm around the SDI pipe and each layer interface.

2.3. Model Quality Evaluation Criteria

The evaluation of the HYDRUS-2D model encompassed the assessment of model quality, sensitivity analysis, and model calibration based on the predicted and measured volumetric soil water content dataset (n=30) of each variant within the upper 12 cm. This specific depth was chosen as it typically represents the rooting zone of cool-season turfgrass [67].

The correlation coefficient (R²), root mean square error (RMSE), mean absolute error (MAE), and Nash-Sutcliffe efficiency (NSE) were calculated using the following equation to assess the model quality parameters:

$$\mathrm{R}²={\left({\displaystyle \frac{S_{XY}}{S_X{S}_Y}}\right)}^2$$

$$RMSE=\sqrt{{\displaystyle \frac{1}{N}}{\displaystyle \sum }_{i=1}^N\left(X_i-Y_i\right)²}$$

$$MAE={\displaystyle \frac{1}{N}}{\displaystyle \sum }_{i=1}^N\left|\left(X_i-Y_i\right)\right|$$

$$NSE=1-{\displaystyle \frac{{{\displaystyle \sum }}_{i=1}^N\left(X_i-Y_i\right)²}{{{\displaystyle \sum }}_{i=1}^N\left(X_i-X_{av}\right)²}}$$

where X_i = the measured volumetric water content (cm³ cm⁻³), Y_i = the predicted volumetric water content (cm³ cm⁻³), X_av = the average of the measured volumetric water content (cm³ cm⁻³), and N = the number of observations.

The correlation coefficient should be close to 1. The root mean square error (RMSE) is a commonly utilized metric for evaluating the agreement between measured and simulated values and should be close to zero. A widely accepted standard that does not consider over- or under-forecasting is the mean absolute error (MAE). In an ideal scenario, the MAE should be nearly zero. Both the root mean square error (RMSE) and mean absolute error (MAE) share the same units as the measured and predicted values [68]. Lastly, the Nash-Sutcliffe efficiency NSE [69] is a normalized statistic commonly used to evaluate hydraulic models and compares residual and measured variance [70]. The parameter equals zero when the square of the differences between the measured and predicted values equals the variability in the measured data. If the NSE value is negative, the measured mean is a more accurate predictor than the model [71]. If the model gives perfect results, the NSE = 1. An acceptable-quality model should have an NSE > 0.5 [68].

2.4. Input Parameter, Parameter Sensitivity and Model Calibration

The water retention drying curves were determined using Eijkelkamp standard sandbox apparatus [72] to assess the water content at pF 1, 1.8, and 2.5. The water content at pF 4.2 was determined using a pressure plate apparatus, as per the DIN EN ISO 11274 (2019) [65]. The soil water retention curves were then parameterized according to the van-Genuchten equation [43] by adjusting the θ_r (residual water content), α, and n (shape factors) values using the EXCEL solver function to match the measured water content and water suction drying curve values. The saturation water content θs was fixed at the total porosity (TP). The material´s hysteretic behavior was evaluated based on the capillary rise of water in the materials employed in this study within experimental containers composed of 10 rigid plastic rings, each measuring 2 cm in height [36]. The rings were filled with the material of interest (95% bulk density of Proctor density D_PR), and a flooding depth of 1 cm was maintained for 48 h before the water content was determined gravimetrically. It was assumed that the water tension in the rings is defined by the distance to the water table, and the water content at equilibrium is analogous to the water retention curve [47]. The wetting water retention curve was initially parameterized by adjusting and establishing the parameter α to α_w [33]. The parameters used in HYDRUS-2D are presented in Table 2.

Sensitivity analyses and model calibration focused on the materials used in the top two layers, including STD, UHFC (rootzone layers), URTS F (intermediate layer), and DS (gravel drainage layer) within each construction method (2A, 2B, and 3) and irrigation system (SDI and SPR). The HYDRUS-2D model was used to analyze the sensitivity of the soil hydraulic parameters θ_r, α, and n (representing the soil water retention drying curve parameterization) while excluding the exact analytically determinable parameter θ_s (fixed at total porosity), as well as the parameters α_w and θ_s^w (representing the opportunities in the HYDRUS code for describing material´s hysteresis behavior). The study evaluated the impact of a 20% incremental increase in these parameters on simulation output (VWC) after 10 mm irrigation (cycle 1) regarding the deviation in model efficiency (NSE).

The HYDRUS-2D software package’s inverse solution function was employed for model calibration. This process utilized the four parameters identified as most sensitive through sensitivity analysis (scenario F1–F4). Additionally, based on previous research [33,45,54,55], a combined calibration approach was applied to the parameters α_w and n, as well as α_w and ϴ_s^w (scenario F5–F6). Six different scenarios (F1–F6) were used for the calibration procedures. The HYDRUS inverse solution was used to independently calculate the calibrated values of the materials used in the top two layers (STD, UHFC, URTS F, DS) for each combination of construction method (2A, 2B, 3) and irrigation system (SPR, SDI).

2.5. Irrigation Scheduling Development

The HYDRUS-2D code was also used in hypothetical instances to examine the effect of irrigation treatments, which are given in Table 3. The parameters for the boundary and flow domain remained unchanged from previous descriptions. Based on previous research [73] that demonstrated the favorable water retention characteristics of the three-layered construction method 3, this methodology was employed to evaluate the efficacy of various irrigation approaches regarding water usage.

The analysis focused on two key factors influencing irrigation efficiency, i.e. drainage flux and soil water storage within the simulation domain (upper 50 cm). Four different irrigation approaches were assessed. Each approach involved applying a total irrigation amount of 10 mm through either sprinkler (SPR 1–4), SDI (SDI 1–4), or HYBRID (HYBRID 1–4) in up to five irrigation events within a 12 h period.

3. Results

3.1. Model Quality Evaluation

The model's performance was evaluated using the model efficiency parameter NSE, and the results were graphically visualized by comparing the observed values against the simulation output differences, presented through spatial maps (ordinary Kriging method).

Figure 3 shows the observed values (left) against the difference values (right) of construction methods 2A, 2B, and 3 under sprinkler irrigation, each averaged across the entire observation time 0-48 hours following 10 mm irrigation. Figure 4 is analogous to Figure 3, but in each case under SDI.

The model’s effectiveness, measured by Nash- Sutcliffe efficiency (NSE), ranged from 0.58 (2B) to 0.68 (2A) and 0.72 (3) for sprinkler-irrigated variants (generally achieving NSE values above 0.5) and from 0.27 (2A) to 0.48 (3) and 0.52 (2B) for SDI variants. The highest model performance was observed in the 3_SPR variant (NSE = 0.72), whereas the poorest performance was observed in the 2A_SDI variant (NSE = 0.27).

Additionally, graphical analysis showed that, for sprinkler-irrigated variants 2A and 2B, the model tended to calculate higher predicted values than observed values. At an observation depth of 6 cm, the average overestimation was 0.85 vol. % for 2A and 1.22 vol. % for 2B. At a depth of 12 cm, the average overestimation was 0.94 vol. % for 2A and 1.55 vol. % for 2B. In the case of SDI, an average overestimation within variant 2A of 1.27 vol. % (6 cm depth) and 3.11 vol. % (12 cm depth) was observed. In contrast, variant 3_SDI exhibited a different pattern, with the model generating lower predicted values than the observed ones. The average underestimation was 1.54 vol. % at a depth of 6 cm and 2.27 vol. % at a depth of 12 cm. Similar to the model’s performance, variant 3 showed the most minor differences among sprinkler-irrigated variants, while 2B displayed the least variation among SDI variants.

3.2. Sensitivity Analysis

The sensitivity analyses of the 5 parameters ϴ_r, α, n (normal drying WRC), ϴ_s^w , and α_w (wetting WRC) presented in Figure 5 illustrate the respective deviations in model efficiency (NSE) under a soil hydraulic parameter perturbation of +20% within each variant's upper two layers (Layer 1 and Layer 2). The results demonstrated that the parameter perturbation elicited contrasting sensitivity responses across the evaluated variants.

Generally, it could be observed that the materials within Layer 1 reacted more sensitively under a soil hydraulic parameter alteration compared to Layer 2 with an average NSE increase (absolute value across all parameters and variants) of 0.14 compared to 0.05, respectively (Table S1 in the supplement).

Further, among the five parameters examined, α consistently exhibited the least sensitivity with a maximum NSE deviation of −0.04 (3_SDI, Layer 1) and 0.04 (2A_SDI, Layer 2), while parameter n demonstrated the highest sensitivity with a maximum NSE deviation of −0.62 (3_SPR, Layer 1). An analysis of parameter n variations revealed that four out of six variants in both layers showed negative NSE deviations. The exceptions were variants 2A_SPR and 2A_SDI in Layer 1 and 2B_SDI and 3_SDI in Layer 2.

Parameter ϴ_s^w had the second highest sensitivity, and the NSE deviation ranged from 0.14 (3_SDI) to −0.33 (2A_SDI) in Layer 1 and from 0.09 (2A_SDI) to −0.11 (2B_SDI) in Layer 2.

Summarizing, based on the absolute averaged NSE deviation (across all parameters and variants), the four most sensitive parameters are as follows (Table S1 in the supplement):

n (0.35) > ϴ_s^w (0.14) > ϴ_r (0.11) > α_w (0.04) within Layer 1;

n (0.11) > ϴ_s^w (0.05) > α_w (0.04) > ϴ_r (0.02) within Layer 2.

3.3. Model Calibration

The results of model calibration calculated using the HYDRUS-2D inverse solution are presented for each material (STD, UHFC, DS, and URTSF) and scenario (F1–F6) in Table 4. The influence of model calibration scenarios F1–F6 on model efficiency (NSE) is illustrated in Figure 6.

Notably, substantial variations are observed depending on whether the calibration is incorporated into Layer 1 (L1) or Layer 2 (L2) or if it is implemented across both upper layers, Layer 1 and 2 (L1+L2). In general, the sprinkler-irrigated variants, which exhibited already Nash-Sutcliffe efficiency (NSE) values exceeding 0.50 under laboratory-measured default settings, demonstrated only a minor improvement across all model calibration scenarios (between +0.07 and +0.01).

Within the SDI variants, generally, model calibration scenario F5 pushed the NSE values to a range of > 0.50, whereby the greatest improvement could be achieved with the combined implementation of both layers (L1+L2), and the values increased in 2A_SDI from 0.27 to 0.79, in 2B_SDI from 0.52 to 0.81, and in 3_SDI from 0.48 to 0.76

The impact of the model calibration scenarios on the average model efficiency values across the SDI variants was ranked in descending order as follows (details in Table S2 in the supplement):

F5_L1+L2 (0.79) > F5_L1 (0.68) > F5_L2 = F4_L1 (0.63).

This indicates, that scenario F5 (i.e. calibration α_w and n) resulted in the highest gain in efficiency and should be used for further simulations.

Figure 7 illustrates the regression analysis of the measured values (n=30) in comparison with the simulated volumetric water content under the default settings (without calibration) and across the highest NSE improvement scenario F5 (sprinkler F5_L1 and SDI F5_L1+L2), including the correlation coefficient (R²), slope of the regression lines (SLOPE), root mean square error (RMSE), and mean absolute error (MAE).

Correlation coefficient (R²), RMSE, and MAE values across the sprinkler-irrigated variants were between 0.72 and 0.79, 1.45 and 1.75 vol. %, and 1.09 and 1.46 vol. %, respectively, under default settings and developed to 0.73–0.83, 1.42–1.49 vol. %, and 1.02–1.15 vol. %, respectively, under model calibration scenario F5_L1.

Within the SDI variants, the values were between 0.44 and 0.81, 1.58 and 3.03 vol. %, and 1.20 and 2.28 vol. %, respectively, under default settings and developed to 0.75–0.84, 1.11–2.18 vol. %, and 0.77–1.71 vol. %, respectively, under model calibration scenario F5_L1+L2.

Generally, the correlation coefficients (R²) were greater across the calibrated SDI variants (0.75–0.84) than across the calibrated sprinkler-irrigated variants (0.73–0.83). Nevertheless, the slopes of the regression lines were closer to 1 within the sprinkler-irrigated variants compared to SDI variants and ranged from 0.81 (3_SPR) to 1.04 (2B_SPR) and 0.65 (3_SDI) to 0.95 (2A_SDI), respectively.

3.4. Irrigation Scheduling Evaluation

The results of the HYDRUS-2D code of the three-layered construction method 3, which were utilized in hypothetical instances to examine four different irrigation approaches (IA1: 10 mm irrigation at time zero; IA2: 75% at time 0 and 25% after 6 hours; IA3: 50% at time 0, and 25% each after 6 and 12 hours; IA4: 50% at time 0, and 12.5% after 3, 6, 12, and 18 hours; see Table 4). The results are shown in Figure 8 (and Table S3 in the supplement). These observations were made at various times (4–48 h) after 10 mm irrigation (sprinkler, SDI, or hybrid). Different behavior was observed across the irrigation approaches concerning the soil storage of irrigation water and drainage flux.

Up to 12 hours (8 hours for IA1), the water storage was very similar within the irrigation approaches with a tendency of generally minimal higher water storage in the SDI variants. Generally, the water storage was higher in IA1 and IA2 compared to IA3 and IA4. Strong differences could be identified between the water storage in the different irrigation systems after 24 hours: the SPR variant always showed the lowest water storage followed by the SDI irrigation. The hybrid irrigation system always had the highest water storage. The hybrid-approach results at 48 h were 3.33 mm (IA4), 3.61 mm (IA3), 2.69 mm (IA2), and 2.39 mm (IA1) (Figure 4). Since the cumulative drainage is naturally determined by the water storage, these values were always inversely related to the storage (Figure 4).

In summary, focusing on the maximum residual soil water storage and minimum drainage flux of each irrigation technique (SPR, SDI, HYBRID) and approach (IA1—IA4) at an observation time of 48 h, the following order was observed for soil water storage:

HYBRID-IA3 (3.61 mm) > SDI-IA4 (2.53 mm) > SPR-IA3 (0.38 mm),

and for drainage flux:

HYBRID-IA3 (2.57 mm) < SDI-IA4 (3.59 mm) < SPR-IA3 (5.46 mm).

The difference of water storage plus drainage flux to the 10 mm irrigation is due to the actual evaporation during the 24 hours.

4. Discussion

Turfgrass areas, regardless of whether they are golf courses or sports fields, offer numerous ecosystem services, including carbon sequestration, oxygen production, water and air purification, heat dissipation, and water infiltration, which are vital for mitigating climate change and improving urban environmental quality [1,2,74]. Additionally, they enhance aesthetic value and support recreational activities, contributing to social well-being in urban landscapes. [9,10].

Nevertheless, adequate drainage is essential for all turfgrass areas, as it facilitates the removal of substantial volumes of water and prevents the cancelation of events due to waterlogging and unstable playing surfaces. Conversely, the rootzone construction method and associated irrigation systems must support efficient irrigation water usage, uniform soil moisture distribution, and adequate soil moisture retention to maintain acceptable turfgrass quality and create a complex system regarding soil water dynamics.

For this purpose, HYDRUS-2D is an affordable and widely utilized numerical predictive model and a valuable tool for evaluating efficient irrigation scheduling. To analyze and enhance the accuracy of this model, we compared experimentally measured greenhouse data with the predicted data of the model concerning irrigation water distribution.

The results of this study demonstrated the efficacy of the HYDRUS-2D model in simulating irrigation water distribution, with the need for model calibration among multilayered rootzone construction methods, particularly in the context of subsurface drip irrigation.

4.1. Theoretical Aspects

The model performance varied depending on the construction method and the irrigation system used. Initial model simulations of water content utilizing default (laboratory-measured) parameters demonstrated acceptable to high performance for sprinkler irrigation (NSE 0.58–0.72) but suboptimal performance for SDI (NSE 0.27–0.52). The model tended to overestimate the two-layered sprinkler irrigation and SDI variants, consistent with Ghazouani et al. [60], who reported overestimation of modeled water content values compared to measured values, with 2A_SDI displaying the lowest model accuracy. In contrast, the three-layered SDI variant notably underestimated soil water content, particularly in the lower parts of the rootzone (observation depths of 6 and 12 cm).

Sensitivity analysis of the parameters indicated that the shape factor n exhibited the highest sensitivity within both layers and across all variants, whereas shape factor α (drying WC) demonstrated the lowest sensitivity, which is consistent with the results of Inoue et al. [58]. Soil hydraulic parameter ϴ_s^w shows the second highest sensitivity, aligning with the findings of Abbasi et al. [59]. The sensitivity analysis further demonstrated the need for the precise parameterization of soil hydraulic parameters, as construction method 2B/Layer 2 exhibited a decrease in the Nash-Sutcliffe efficiency (NSE) value of −0.23 when the shape factor n was increased by 20% under sprinkler irrigation. In contrast, this modification increased the NSE value of +0.19 under subsurface drip irrigation (Figure 5).

The calibrated α_w (wetting WRC) values resulted in α_w/α ratios averaging 1.90 across all variants (calibration scenario 1), indicating substantially greater α_w values than the corresponding drying curve values. This finding aligns with the study of McCoy and McCoy [45] and approximates the value of 2 (α_w/α ratio) suggested by Kool and Parker [54] for estimating wetting curve parameter α_w. This supports using this approximation when measured wetting curve data are unavailable. The shape factor n shows changeable behavior within calibration scenarios F2 and F5 with a partially low to strongly pronounced increase and decrease depending on the construction method and associated irrigation system; these findings do not agree with the results of McCoy and McCoy [45], which showed lower n values when considering the typical characteristics of the hysteresis wetting curve response. Nevertheless, our findings suggest that, particularly under subsurface drip irrigation conditions, the shape parameter n_w(wetting curve) should be additionally considered within the HYDRUS code, along with α_w, to describe the hysteresis behavior of the material more accurately.

Following calibration, the model strongly agreed with the measured data, achieving a Nash-Sutcliffe efficiency (NSE) value of up to 0.81 for subsurface drip irrigation (SDI) and 0.75 for sprinkler irrigation.

Calibration demonstrated that under sprinkler irrigation, the variants exhibited NSE values > 0.5 already under laboratory-measured default settings and revealed only slight improvement after calibration. Calibration is therefore not necessary, at least under the conditions investigated. If calibration is to be carried out anyway, consideration of layer 1 is sufficient.

In contrast, under SDI, a significant improvement in model quality could be observed by considering both layers (Layers 1 and 2). A plausible explanation is that during SDI, the model primarily utilized the wetting retention curve in Layer 2. However, during drainage and redistribution, both layers and their corresponding wetting and drying curves and the connections between them necessitated consideration within the model calibration procedure.

The soil hydraulic parameters for each material were independently calibrated using the HYDRUS inverse solution method. However, the results demonstrated that the model structure significantly influenced the calibrated material's soil hydraulic properties, as evidenced by the high variability in the calculated soil hydraulic parameters. For instance, across construction methods 2A_SPR and 2A_SDI, the shape factor α_w exhibited values ranging from 0.152 to 0.388 (Table 4). Furthermore, the models revealed pronounced gradients at the interface between neighboring layers (Layers 1 and 2) due to markedly different soil physical properties of the materials employed. For instance, in construction method 3, Layer 1 (UHFC) demonstrated a saturated hydraulic conductivity (K_sat) of 649 mm h⁻¹, while Layer 2 (URTSF) showed a value of 1465 mm h⁻¹. Similarly, substantial water content and water suction gradients are observed at the interface of Layer 1 and 2 in construction method 2A, where Layer 1 exhibited a K_sat of 220 mm h⁻¹ and Layer 2 had a K_sat of 916 mm h⁻¹. These findings indicate that complex model structures involving multilayered rootzone construction methods comprising adjacent layers with significantly different soil physical properties necessitate model calibration, particularly under subsurface drip irrigation conditions.

Notwithstanding the observed local variations, the NSE values consistently remained positive, indicating substantial agreement between the measured soil water content and the model’s simulated values [69] despite the complexity of the conditions to which the model was subjected, that is, multi-layered rootzone construction methods and irrigation systems with high spatial variability.

4.2. Practical Significance

Despite discrepancies between the calculated and measured values of soil water content, these variations typically remained below 4 vol. % and are consistent with standard field measurement uncertainties. These uncertainties are significantly influenced by the specific measurement technology employed (e.g., measurement with a frequency domain reflectometry (FDR) sensor), the measurement point's location, and the measurement's timing. Regarding practical significance, it is essential to note that the model demonstrates efficacy in multilayer construction methods with highly complex soil physical relationships and under irrigation systems with considerable spatial variability (i.e., SDI), significantly influencing irrigation water dynamics [75]. Consequently, this tool facilitates the adaptation of irrigation management to the specific onsite structural conditions of sports turf areas and enables the implementation of efficient irrigation practices.

In the context of irrigation scheduling evaluation, the model data of the three-layered variant (construction method 3) indicated that, under sprinkler irrigation, dividing the total irrigation water quantity into multiple smaller applications (approaches 3 and 4) generally resulted in higher soil water storage and lower drainage flux compared to single large applications. Furthermore, the model data under the hybrid irrigation approach indicated the most efficient use of resource water, characterized by the highest soil water storage and lowest drainage flux compared to the other irrigation approaches (Figure 8). The practical applicability or the requirement profile of an irrigation technique should also be emphasized, as maintenance practices typically conducted on turfgrass areas necessitate surface irrigation (i.e., sprinkler) for establishing seeding or flushing in granulated fertilizers. In the context of turfgrass irrigation, the primary focus remains on the targeted application of irrigation water into the root zone, which should ideally be applied with minimized evaporation loss, surface runoff, and moisture penetration of the turf to reduce disease pressure. The results indicate that under a hybrid irrigation approach, these requirements are attainable, further demonstrating the potential of HYDRUS-2D as a cost-effective and valuable tool for evaluating and optimizing turfgrass irrigation strategies.

4.3. Limitations of the Study and Further Research

While the calibrated models exhibited strong overall performance, some discrepancies between measured and predicted values remained, particularly for the three-layered construction under SDI. This suggests that further refinement of the HYDRUS code, such as implementing n_w as an additional shape factor, may be necessary to account for the material's hysteresis behavior more precisely. Moreover, to evaluate the fundamental model performance across various rootzone construction methods and irrigation systems, the study was conducted on bare soil profiles, and the incorporation of turfgrass and its associated root water uptake patterns would undoubtedly influence water distribution and irrigation efficiency.

Future research should incorporate root water uptake models to provide a more comprehensive understanding of soil-plant-water dynamics and relationships in multi-layered construction methods for turfgrass areas, wherein an increasing complexity for model calibration could be anticipated.

5. Conclusions

This study demonstrated the efficacy of the HYDRUS-2D model in simulating irrigation water distribution across various turfgrass rootzone construction methods. The key findings include the following: (i) The initial model performance varied depending on the construction method and irrigation system, with acceptable to high performance for sprinkler irrigation (NSE 0.58–0.72) but suboptimal performance for subsurface drip irrigation (SDI) (NSE 0.27–0.52). (ii) Sensitivity analysis revealed that the shape factor n exhibited the highest sensitivity, whereas the shape factor α showed the lowest sensitivity across all variants. Model calibration significantly improved the performance, achieving Nash-Sutcliffe efficiency (NSE) values up to 0.81 for SDI and 0.75 for sprinkler irrigation, respectively. This calibration appears necessary to enhance model quality among multilayered rootzone construction methods, particularly in the context of subsurface drip irrigation. The calibrated α_w/α ratios averaged 1.90 across all variants, aligning with previous research and supporting this approximation when the measured wetting curve data are unavailable. (iii) The evaluation of irrigation scheduling approaches revealed that dividing the total irrigation water quantity into multiple smaller applications generally resulted in higher soil water storage and lower drainage flux than single large applications. (iv) A hybrid irrigation approach combining sprinklers and SDI systems showed the most efficient use of water resources.

While the calibrated models demonstrated good overall performance, some discrepancies between the measured and predicted values persisted, particularly for the three-layered construction under the SDI. This finding suggests that further refinement of the model is necessary. These findings indicate that HYDRUS-2D has the potential to be an efficient tool for evaluating irrigation strategies in turfgrass areas. However, future research should incorporate root water uptake models to provide a more comprehensive understanding of soil-plant-water dynamics in multilayered construction methods for turfgrass areas.

Overall, this modeling approach has the potential to optimize irrigation scheduling, thereby enhancing water-use efficiency in turfgrass management. Consequently, this optimization improves the sustainability, capacity for carbon sequestration, and ecosystem services of turfgrass, which are particularly crucial in urban settings.