Simulation of irrigation in southern Ukraine incorporating soil moisture state in evapotranspiration assessments

Abstract


Introduction
The accuracy of evapotranspiration estimates is one of the determining factors for performing forecasts of soil moisture state in the irrigation management process (Wanniarachchi and Sarukkalige, 2022). In the areas of agricultural production, evapotranspiration can be considered as consisting of three components: evaporation of intercepted moisture, transpiration, and evaporation from the soil surface (Savenije, 2004). Changes in the corresponding fluxes are mutually influenced by fluctuations in meteorological parameters, vegetation dynamics, and soil moisture (Rodriguez-Iturbe, 2000;Shao et al., 2017). Thus, to quantify the intensity of evapotranspiration for more accurate modeling of moisture transport in the "soil-plantatmosphere" system, integrated models that consider soil and atmospheric physics along with plant physiology are needed (Overgaard et al., 2006). A widely used evapotranspiration model based on the Penman-Monteith equation (Monteith, 1965) successfully estimates it in the case of closed vegetation for various weather and soil conditions (see, e.g., Shao et al., 2022). In it, the processes of transpiration and evaporation from the bare soil, which are different intrinsically, are not considered separately. One of the first models in which the description of these processes was separated was the Shuttleworth-Wallace model (Shuttleworth and Wallace, 1985). Among the disadvantages of this model, the need to determine the values of numerous parameters that limits its application can be singled out (Gharsallah et al., 2013). The compromise between the model's complexity and prediction accuracy is provided, in particular, by the Priestley-Taylor energy balance model (Priestley and Taylor, 1972). The modified Priestley-Taylor model can also effectively evaluate evaporation and transpiration separately using the data on the downward energy flux (Qiu et al., 2019), which can be considered as the main factor determining the intensity of crops evapotranspiration (Gong et al., 2021).
In the studies of irrigated crops' growing processes, the best strategy according to Faybishenko (2007) is to choose the methods that take into account the highest number of input parameters. At the same time, the problem of determining the accuracy of evapotranspiration estimates remains urgent in each specific situation and such approaches as the usage of different methods' linear combination with fittable coefficients (Romashchenko et al., 2020) or machine-learning-related approaches (Elbeltagi et al., 2023) are used. The need for scenario modeling here stems from the fact that under irrigation the maintained ranges of moisture content in soil's root layer significantly influence the availability of moisture for plants. Hence, the processes that have a decisive influence on evapotranspiration and the parameters that quantify their intensity may change.
In this paper, evapotranspiration estimates are used as input to the models of moisture transport based on the Richards differential equation. The study of their accuracy is carried out by assessing the compliance of the simulated dynamics with sensor readings. The model parameters were identified using the readings in the initial part of the growing season. Then, to assess the accuracy, we perform extrapolation modeling over longer time ranges including the entire growing season.

Material and Methods
We investigate two models for evapotranspiration assessment -the Priestley-Taylor and Penman-Monteith methods -along with their modifications that consider the current state of soil moisture.

The Penman-Monteith method
Having a weather station equipped with sensors of temperature and relative humidity of air along with solar radiation and wind speed, estimation of potential evapotranspiration by the Penman-Monteith method can be performed the following way (Cannata, 2006 Soil heat flux G is set to be linearly dependent on the flux of solar radiation n R (mJ m -2 s -1 ).
The atmospheric component of evapotranspiration is calculated as (Cannata, 2006 Soil surface resistance rs can be estimated based on the so-called "relative soil moisture index" Rsm2 (Sellers et al., 1992;Kustas et al, 1998), which is defined as the ratio between the volumetric moisture content on the soil surface ( where Rns and Rnc are the energy (W m -2 ) received by soil and vegetation surfaces; αs and αc are the coefficients for λEs and λTr. Rns, Rnc, αs, and αc can be defined as (Gong et al., 2021) ,, kLAI ns n nc n ns where τ is the fraction of radiation that reaches the surface of soil; LAI is the Leaf Area Index; αs0 and αc0 are the coefficients for soil and vegetation defined as (Gong et al., 2021) where α0 is the reference coefficient (1.26), τc is the critical value of τ when the closure of vegetation is maximal.
Leaf aging coefficient fs can be determined as (Gong et al., 2021)  where Topt is the optimal temperature for crop growth ( • С).
The soil moisture stress index fsw is used to combine the data on soil and atmospheric conditions and can be determined according to the model presented in Deardorff (1978) as (Gong et al., 2021): is the effective moisture saturation in the upper soil layer with the depth of 0.1 m; θ, θs, and θw are the actual volumetric moisture content, saturated moisture content, and wilting point, correspondingly.
To finally obtain the expression for αe, G is described as a function of Rns (Choudhury et al., 1987) where fG is the fraction of G in Rns. Combining all the above mentioned, for αe we have (Gong et al., 2021) The method described by Venturini and co-authors in Venturini et al. (2008) modifies the Priestley-Taylor SM is the volumetric moisture content in the soil, and SMc is the field capacity.

Methods for modeling moisture transport
The main aim of our study is to experimentally test the effectiveness of the combined use of the abovementioned evapotranspiration models and differential moisture transport models while simulating changes in moisture content in irrigated soil. For this purpose, considering the Penman-Monteith method, we determine the values of s r according to (1) with a=8.2, b=5.9 according to Kustas et al. (1998). The value of Δ is multiplied by F calculated according to (2). As an alternative, we use the Priestley-Taylor model with similar modifications. To model moisture content dynamics, we consider the classical integer-order one-dimensional head-based Richards equation according to the method described in Romashchenko et al. (2020). The corresponding equation has the form H is the hydraulic conductivity (m s -1 ). In (3)  At the upper boundary 0 z  , in the case when the soil is saturated, the Dirichlet boundary condition is set (van Dam and Feddes, 2020). In other cases the Neumann boundary condition is set in the form The function S that models water uptake by plant roots has the form (Molz and Remson, 1970) where v is the depth of the root layer, ( ) L z is the function of root length distribution density, T is the transpiration.
With the known value of actual evapotranspiration E T , it is subdivided on the components e Q and T according to Gigante et al. (2009) the following way: The finite difference method (Samarskii, 2001) is used to numerically solve the initial-boundary value problem for Equation 3 as described in Romashchenko et al. (2020). Additionally, we consider the one-dimensional space-time-fractional equation of moisture transport in the form (Romashchenko et al., 2021)  are the Caputo derivatives of fractional order subject to time t and depth z (Podlubny, 1999).
The numerical technique for the fractional-differential model is described in Romashchenko et al. (2021).
In computational experiments, soil's water retention curves are determined by selecting the parameters of the van Genuchten model (van Genuchten, 1980)  is the fixed power. To compensate for the errors in the measurement of irrigation water flow, precipitation flow, and evapotranspiration estimates, the corresponding flows were multiplied by coefficients, the values of which are, similarly to the described in Romashchenko et al. (2020), fitted by the particle swarm optimization (PSO) algorithm (Zhang et al., 2015). These coefficients are fitted in a way to minimize the total sum of squares deviations of the simulated water head dynamics from the measured one. To compensate for the errors in laboratory determination of the saturated hydraulic conductivity, as well as errors arising from soil heterogeneity, the value of the saturated hydraulic conductivity was also determined by the PSO algorithm. An approach for assessing the impact of evapotranspiration models' accuracy on the effectiveness of irrigation management The use of the considered technique in irrigation management involves estimating the time and rate of subsequent watering by predictive modeling with evapotranspiration estimated using forecast weather data. To estimate the seasonal irrigation rate and the volume of actual evapotranspiration the irrigation simulation is performed for the entire growing season with watering assigned to maintain in the given range the average water head in the root layer. The seasonal effectiveness of irrigation can be assessed using the so-called relative yield function (Kovalchuk and Matiash, 2006), which simulates the decrease of yield due to not-optimal irrigation. In particular, according to Kovalchuk and Matiash (2006), for maize grown in southern Ukraine, such a function has the form 2 ( , , ) 0.444 2.02 0.556 where u is the actual seasonal irrigation rate (mm), w is the biologically optimal irrigation rate (mm) defined as an irrigation rate for a regime, in which evapotranspiration is maintained at the level close to the potential one, p is the precipitation (mm). Given the close relationship between w and evapotranspiration, we transform (4) where W is the actual seasonal water supply (mm), E T is the total potential evapotranspiration (mm).
After conducting a scenario modeling of irrigation assignments throughout the growing season and obtaining simulated values of W and E T , the effectiveness of the method of evapotranspiration assessment for irrigation management can be estimated by the value of ( , ) f W ET . The height of plants for the whole period was assumed to be equal to 1 m, and the depth of the root system was assumed to be 0.5 m.

Input data
In the scenario modeling of the entire growing season using the data collected in 2021, irrigation was applied when the average relative volumetric moisture content in the root layer of the soil decreased below 70% of field capacity that corresponds to the moisture content level of 24.8%. Irrigation was simulated until the corresponding level rises above 100% of field capacity (35.4%). The simulation was performed using the data acquired during the period of intensive irrigation from 25.07.2021 to 01.09.2021. In the period from 19.05.2021 to 25.07.2021 water supply to plants was provided mainly by precipitation. To test the sensitivity of computational procedures to inaccuracies in forecast meteorological data, modeling was also performed using the data on temperature, humidity, and wind speed, from the weather station located (46°51'N, 34°24'E) at a certain distance from the field. In the absence of predicted data on solar radiation, it was assumed that accurate measurements are performed once per 5 days, and then a constant value is used in the simulation. The soil at the experimental sites corresponds to the southern low-humus heavy loam chernozem on loess. A single-layer soil model with the following values of the van Genuchten model's parameters was used for the data collected in 2021: θ0=0.094, θ1=0.5059, α=0.00919, n=1.475, m=0.3223. Parameters' values were obtained using Rosetta software based on soil particle size distribution data. The actual irrigation regime, according to the data obtained in 2021, allowed the decrease of moisture content down to 68% of field capacity compared to the maximum reduction to 78% of field capacity according to the data collected in 2018. The minimum recorded value of suction pressure was -132 kPa compared to -40 kPa according to the data collected in 2018.
To identify the parameters of the models in the case of the data collected in 2018, a time interval of 11 days from the initial moment was used. For the data collected in 2021, we used a 7 days interval. For both cases, 273 one watering was carried out within these intervals. An interval of 50 days was used to test the influence of different evapotranspiration assessment methods' usage on modeling accuracy. The population size of the PSO algorithm was 20 particles, 20 iterations were performed with the following values of the parameters: ω=φr=φp=0.8. The smallest errors for the data set collected in 2018 were achieved when modeling a domain with the depth of 3 m with the Dirichlet condition H=-3.8 m at its lower boundary on the finite-difference grid with 50 nodes.
In the case of data collected in 2021, a 1 m deep domain was used with the Neumann condition at the lower boundary on the finite-difference grid with 20 nodes.

Results and Discussion
The root-mean-squared error (RMSE) and average relative error ( are the measured and the modeled water head values, N is the number of measurements) for the intervals of 11 and 50 days are given in Tables 1 and 2. The results show that for 2018 the differences between the modeling errors using different algorithms for estimating evapotranspiration are insignificant. Given the maintained high level of soil moisture content, its consideration has no significant impact on the obtained dynamics of evapotranspiration. Regarding the use of the model that contains the derivatives of fractional order, the results confirm the conclusions given in Romashchenko et al. (2021). Thus, the accuracy of the parameters' identification for the fractional-order model is ~10% higher, but when modeling for longer intervals it decreases faster than in the case of the classical model. Performing simulation based on the data collected in 2021 we observed up to ~10% reduction in RMSE using the modifications of the Penman-Monteith and the Priestley-Taylor methods that take the data on soil moisture content into account. This decrease can be explained by a larger range of water head changes than in the case of the data collected in 2018. Accuracy when using the Priestley-Taylor model was higher here and an additional increase in accuracy when using the fractional-differential model was observed for both considered intervals. Thus, in regard of RMSE and average relative errors when modeling one or several irrigation cycles in two cases of different average moisture content and corresponding pressures we confirmed that, even when input data are collected in production condition, the incorporation of soil moisture assessments into the considered evapotranspiration models allows obtaining the increase in simulation accuracy in accordance with experimental results (see, e.g. Ding et al., 2013;Gong et al., 2021) on the importance of soil moisture factor. Some of the obtained results on the dynamics of water heads are shown in Figure 1, 3 and 4. The dynamics of average volumetric moisture content in the 0.5 m layer of the soil is given in Figure 2 and 5. The results show that taking into account soil moisture in the Penman-Monteith method leads to an overall increase in the simulated moisture content level compared to the basic version of the method (Figure 2 and 5). When using the Priestley-Taylor model, the opposite trend was observed. The reasons for the latter could be the error accumulation in the performed long-range simulations that, in turn, could be caused by lower accuracy of the Priestley-Taylor model in water stress condition as reported in Shao et al. (2022). In the case of the data collected in 2021, there was an underestimation of water intake in computational experiments while performing predictive modeling. This is the main reason for the greater efficiency in this case of the modified Priestley-Taylor model. As can be seen from Figure 4, the reason for the high modeling errors is the lack of response of the sensor located at the depth of 0.4 m to the second and third irrigation. A similar trend was observed for the sensor located at the depth of 0.6 m. The simulation results were consistent with the data of the sensor located at the depth of 0.2 m ( Figure 3) and with subsequent dynamics of water head changes according to the other two sensors. Thus, the used modeling procedure has stable response to inaccuracies in input data subject to the change of evapotranspiration assessment method similarly to the reported about other factors in Bohaienko et al. (2022).

Influence of evapotranspiration models' accuracy on the efficiency of irrigation management
The values of total water supply, total actual evapotranspiration, and the relative crop yield function (5) in the scenario modeling of the 2021 growing season are given in Table 3. At the maximum maize yield equal to 14.5 t/ha (UIPVE, 2019); total evapotranspiration, calculated from the actual meteorological data according to the Penman-Monteith equation; and the actual measured water supply, the expected yield according to Equation (5) differs from the actual one by 1.5%, which confirms its applicability in the considered case. The simulation technique's sensitivity to the meteorological data when using the Penman-Monteith method was low (not more than 4% deviation in the seasonal parameters between the cases of actual and forecast meteorological data). When using the Priestley-Taylor method, it increases (deviation <9%). The use of the considered modifications of evapotranspiration assessment methods leads to seasonally simulated scenarios with lower levels of both water supply and evapotranspiration compared to the usage of original methods. This decrease is more significant when using forecast meteorological data, to which the modified formulas are more sensitive (average deviations of parameter values equal to ~14%). The values of the relative crop yield function when using modified methods are 24-61% higher.

Conclusion
The results of water head dynamics modeling under sprinkler irrigation according to the two data sets collected growing different crops in different meteorological conditions demonstrate the possibility to increase modeling accuracy by ~20% using the Priestley-Taylor method modified to take into account the data on soil moisture content when compared with the classical Penman-Monteith method. The efficiency of this scheme increases with the increase of the range in which moisture content changes. For the data set collected under the irrigation regime aimed at maintaining moisture content in the root layer in the range of 80%-100% of field capacity, modified methods did not improve the accuracy of modeling. These results confirm the well-known conclusion that using such irrigation regimes, moisture content level is optimal for the grown crops, and evapotranspiration is maintained at the level close to the potential one. The results of scenario modeling for the entire growing season and the estimation of yield under the formed water regime showed that, according to the used yield model, the application of evapotranspiration estimates that incorporate moisture content data generates irrigation regimes with lower water supply.