Computational Framework for Combined Simulation of Moisture Transport and Root Growing Under Drip Irrigation

Abstract

This paper investigates the simultaneous modelling of moisture transport in soil and the growth of the root system of plants (specifically, corn) under drip irrigation. To describe moisture transport, a modified Richards equation is used, which, in particular, includes fractional-order derivatives to describe the fractal properties of the soil. We use a stochastic model of root system growth that takes into account the influence of moisture potential field on the growth process. Optimization methods, such as particle swarm optimisation (PSO) and genetic algorithms, are used to calibrate the model. The results of the analysis of the dataset that contains soil moisture content measurements acquired during cultivation of corn show that the combination of integeraccuracy in predicting the dynamics of the moisture potential field and, accordingly, its availability to plants on the test dataset. The best accuracy in the case of the training dataset was provided by the fractional-order model. The results of the study demonstrate the potential of the proposed approach for its application in decision support systems in agriculture. Moisture transport, root growth model, Richards equation, fractional-order differential models, particle swarm optimization, genetic algorithms 1ICST-2025: Information Control Systems & Technologies, September 24-26, 2025, Odesa, Ukraine * Corresponding author. These authors contributed equally.

1. Introduction

Drip irrigation can lead to significant moisture distribution heterogeneity in soil profile. In this regard, within the framework of the irrigation management methodology, which ensures the specified levels of soil moisture availability to plants in root-containing zones, it is important to simultaneously predict both the moistened zones and the root-containing zones with the distribution of the root system in them. Therefore, their scenario modelling based on the mathematical models calibrated for a specific soil and crop conditions has significant advantages over exclusively experimental studies, providing the opportunity to evaluate soil moisture distributions under different weather conditions, irrigation rates and organization of irrigation system.

Such modelling is the basis of decision support systems (DSS) (see, e.g., the review [1]), in Such DSSs can be considered as optimization superstructures over moisture transport models [3], which are based mainly on the Richards differential equation [4] in a two-dimensional statement, or are formed from experimentally determined wetting contours [5]. To further increase the accuracy of modelling, more complex models can also be used. They include, in particular, the models that take into account fractal properties of soils [6]. An alternative to such models, which describe physical processes in soils and are sensitive to the accuracy of parameters measurement, is the use of machine learning algorithms (see, e.g., [7]), the disadvantage of which, in turn, is the requirement for large amounts of input data that are not always available in production conditions.

One of the directions for increasing the accuracy of modelling soil moisture content and moisture consumption by plants during irrigation is the simultaneous modelling of moisture transfer and the growth of root systems, the peculiarities of development of which is influenced by moisture availability [8]. The classical approach here is to represent the root system as a tree with subsequent transformation of such a discrete model into a continuous function of root system density, which is used in moisture transport models that consider the soil as a continuous porous medium. The study of various approaches to such a transformation, their impact on the accuracy and speed of modelling, is studied in [9]. More complex models of this class can contain a description of the processes of moisture movement in the plant tissue and the interaction of the plant with the atmosphere in the form of differential equations [10]. They can also contain a description of the influence of chemical compounds concentrations on the processes of root system growth and its interaction with the soil [11].

It should be noted that the study of the practical application of the above-mentioned models and their testing using real measured data are limited. Mathematically, this requires solving inverse problems for parameter identification, the set of which is often limited exclusively to the parameters of the moisture transport equation [12]. Given the complexity of such problems, soft computing methods, such as genetic algorithms, are used to solve them [13].

Considering the case of the experiment in production conditions with limited capabilities for monitoring the state of the soil (the availability of automatic measurements only of the suction pressure of soil moisture), in this study we build an algorithm for calibrating the moisture transport model, combined with a stochastic model of root system growth, in order to most accurately determine the moistened zones, the size and distribution of the root system for specific observed conditions. In our case, the root system growth model takes into account the influence of soil moisture potential field obtained according to the moisture transport model. The latter model, evapotranspiration, specifically moisture consumption by plants. A combined model calibrated based on the measured data for one growing season can be used in future to refine irrigation regimes and the parameters of irrigation pipeline placement and the distant between emitters on them when growing the same crop in subsequent years.

2. Moisture transport simulation framework

Governing equation for moisture transport. For the purpose of modelling, we use the Richards equation [4] stated in terms of water head in a two-dimensional approximation, similar to the presented in [14]:    H H H  C(h) +  s Ss  t = Dx(,L)x (kx (H ) x ) + Dz(,L)z (kz (H ) z ) − S, 0  x  Lx , 0  z  Lz , t  0 ( 1 ) P(x, z,t)

 g where h(x, z,t) = is the water head, m, H (x, z,t) = + z is the full moisture potential, m, P(x, z,t) is the suction pressure, Pa,  is the density of water, kg/m3, g is the  acceleration of gravity, m/s2, C(h) = is the differential soil moisture content, %/m,  (x, z,t) is h the volumetric soil moisture content, %,  s is the moisture content in water-saturated soil, %, Ss is the specific storage, 1/m, kx (H ), kz (H ) are the hydraulic conductivities in the corresponding directions considering the soil heterogeneous, m/s, S (x, z,t) is the source function, %/s, which models the extraction of moisture by plant roots and its supply by subsurface drip irrigation. P(x, z,t)  g Here Dx(,L) is the two-sided fractional derivative of Caputo type of an order  w.r.t. the variable x (a derivative w.r.t. the variable z is defined similarly) in the following form, which was used in the derivation of the fractional-order mass conservation law in [15]:

Dx(,L) H = (22− ) (x −1Dx,l H + ( L − x) −1 Dx,r,L H ), Dx,l H = (11− ) 0x Hx (x − )− d , Dx,r,L H = (11− ) Lx Hx ( − x)− d .

ion. Water retention curves  (h) of the soil are represented in ( 1 ) according to the van Genuchten model [16] in the form  (h) = r + ( 2 ) 1 + (100 h )n 1−1/n .

 

The soil is modelled as a layered structure with the values of the coefficients r; s;  , 1 / cm ; n that change from layer to layer. The dependency of the hydraulic conductivity k (h) on water head is represented according to the Mualem's model [17] in the form  s − r k(h) = k f  (h) 1 − (1 − n/(n−1) (h))1−1/n 2 , (h) =  (h) − r    s − r where k f , m/s, is the hydraulic conductivity of saturated soil (filtration coefficient),  is a fixed exponent.

The values of the coefficients in the models ( 2 ) and ( 3 ), their initial guess at least, are obtained using the least-squares fitting to the data of experimental studies conducted using a technique described in detail in [18]. Such a technique allow for the determination of k(h) only in a horizontal dimension. Thus, we consider a possible heterogeneous case assuming that kx (H ), kz (H ) are represented by Eq. ( 3 ) with  that is equal in both directions and k f that is different.

The forms of boundary conditions are given similarly to [19]. They include only gravitational flow condition h = 0 on the bottom of the domain; symmetric flow conditions H = 0 on its left z x and right side; and the condition of flux-controlled interaction with the atmosphere on the upper side: H z z=0 kz (h, z, x) = Qe (t, x) − Qp (t, x) ( 3 ) ( 4 ) where Qe (t, x) , Qp (t, x) are the fluxes, m/s, of evaporation and precipitation.

Several assumptions are made about evaporation and precipitation fluxes:

When h(x, 0,t)  Hmax , Qp (t, x) is set to 0: we assume that there is a surface water run-off is water pressure on the surface exceeds a given threshold Hmax ; 2. When soil surface is in a non-saturated state, we assume that the evaporation cannot exceed its maximal possible value Emax (x,t) = −kz (h(0,t), 0, x)h(x, 0,t) ; 3. the spatial extent of the canopy equals to the spatial extent of the root system.

The function S models the extraction of moisture by the root systems of plants the way it is described in [19]. The distribution of transpiration along the depth z is described according to [20] in the form

T (t)L(z) Sz (z,t) = zr  L(z)dz 0 where zr is the depth of the root-containing layer, T (t) is the transpiration rate, m/s.

As the default case, we set L(z) in the form [21] default case, we set the width of the root system equal to rp and assume that its density decreases linearly subject to the horizontal coordinate x [19]:

 rp − | x − xr | , rp − | x − xp | 0, Sx (x) =  rp2

 0, rp − | x − xp | 0.

The total moisture extraction function in the default case has the form [19]

ST (x, z,t) = Sz (z,t)Sx (x) .

To model subsurface drip irrigation we add to ST (x, z,t) the density of irrigation water flow Qss (x, z,t) = Qss0 (t) (xss ) (zss ) where Qss0 (t) is flow density from one emitter, 1/s, xss , zss are the  () is the Dirac delta ( 5 ) ( 7 ) ( 8 ) ( 9 ) ( 10 ) function. Finally, we obtain [19]

S (x, z,t) = ST (x, z,t) + Qss (x, z,t) .

To subdivide evapotranspiration ET into the evaporation flow Qe and transpiration T we use the value of an empirical parameter  dependent on the Leaf Area Index [22]:

T = (1 − e− )  ET , Qe = ET − T .

ET is calculated within the considered framework using the variety of the Penman-Monteith method described in [23].

The two-dimensional model based on ( 1 )-( 10 ) assumes that the distance between the emitters is sufficient for the formation of uniform wetting in the plane along the pipeline.

Numerical technique for the direct problem. The discretization of the model ( 1 )-( 10 ) is performed using a finite-difference scheme on a grid, which is uniform w.r.t. the space variables and non-uniform in time. The solution technique is described in detail in [14].

Parameters identification problem. In order to adapt the model to the actual plant growing and soil conditions, fitted multipliers are introduced to some of the inputs: kET for ET, k prec for precipitation, and kirr for irrigation flow [24].

To take into account errors due to the difference in the behaviour of collected soil monolith samples and the actual soil in the field, additional parameters that can be fitted include two parameters, whose determination accuracy is the lowest: the filtration coefficient k f , or two coefficients k f ,x , k f ,z in two dimensions, and the residual moisture content  r . Additionally, the parameters  , n that influence the shape of  (h) described by ( 2 ) can also be fitted along with an empirical parameter  used to split ET into evaporation and transpiration according to ( 10 ).

The fitting procedure assumes that the values Hi , i = 1,...., N of water heads are known in the moments of time Ti in the points (xi ,zi ),i = 1,...., N Qp (t) , and Qss0 (t) . The goal function to be minimized has the form

1 N  H (xi ,zi ,Ti , v) − Hi , Hi  0, F2 (v) =   Hi (12)

N1 i=1 

 0, Hi = 0 where N1 is the number of non-zero measured water head values.

Choosing one of the mentioned goal functions, the fitting is performed using the Particle Swarm Optimization (PSO) algorithm given in detail in [24].

3. Root system growth simulation framework

Root growth model. In this paper we use a root growth model derived from the one described in [25]. In it, each root is represented by a tree of branches. Each branch is represented by a list of segments, where each segment corresponds to branch growth during one time step.

For the sake of simplification, considering only the case of one modelled plant, on each time step, while solving the problem ( 1 )-( 10 ), two root growth procedures are performed for each branch in the tree that describes a root: the production and the branching.

The production procedure consists in adding a segment to a branch if the length of the branch after the growth is less that the given threshold Lmax . Denoting as ( x0 , y0 ) the final point of the last segment of the branch, the newly created segment spans from ( x0 , y0 ) to a point ( x1, y1 ) = ( x0 , y0 ) + pst ((dx , d y ) + ( gx , g y )) m/s, t is the time step, ( gx , g y ) is the gradient of moisture potential field calculated according to the model ( 1 )-( 10 ),  is the coefficient of root hydrotropism intensity.

The branching procedure is applied at each time step to all root branches, for which the following conditions are fulfilled: where ( dx , d y ) ps is the growth speed, the level of the branch in the tree is less than a given number  max ; the number of sub-branches connected to the branch is less than a given number bmax ; the length of the branch is bigger that la + lb + (nbr + 1)ln where la , lb are the lengths of the apical and basal zones where branching does not occur, ln is the average spacing between sub-branches.

The branching procedure consists in uniformly random selection of a branching point pbr located between l = la and l = lbr − lb where l is the length of the branch up to the selected point and lbr is the total length of the branch; uniformly random selection of a branching angle b between the newly created branch and its predecessor in the range − / 2, / 2 where  is the given maximum branching angle; • • • • • creation of a new branch, which is connected with its predecessor within the tree and contains one segment that starts at the point pbr , is directed according to the selected angle b , and has the length pst .

After modifying the root system model in the form of a tree, its representation in the model ( 1 )-( 10 ) in the form of a root density function is recalculated assuming that the length pst of newly created segments is less than the size of a cell in a finite-difference grid: • • the depth and the width of the root system are calculated from the extent of all points of the moisture extraction function is considered equal to ( x, z )  ci, j where ci, j is the (i, j ) cell of the grid, Tp is the total number of points in root ST (x, z,t) = T (t)N p,i, j / Tp for s, N p,i, j is the number of points that are located within ci, j .

Simulation procedure and inverse problem. The question posed in this paper is whether it is possible, using the above-described procedure for modelling the growth of the root system, to obtain such a form of the root system that would improve the accuracy of modelling the dynamics of water head according to the goal function ( 11 ) in comparison with the use of generalized root system density functions, in particular ( 6 )? Will the simulated form of the root system correspond to the actual one?

To answer the first question, the following modelling procedure is proposed.

In the first stage, the model is calibrated using the root system density function ( 6 ) based on suction pressure measurements collected during several irrigation cycles mid-season

In the next step, assuming as fixed the pre-calibrated values of the parameters of the model ( 1 )-( 10 ), we search for a density function that minimizes the goal function ( 11 ) by generating root system models.

The calculation of the goal function values is as follows. Having fixed values of the parameters of the root system growth model, its development is simulated starting from the start of vegetation until the start of suction pressure measurement. In this case, it is assumed that the pressures in the modelled root-containing zone are maintained by irrigation in a given range. It is assumed that at the start of measuring the pressures for the purpose of modelling, the root system has reached its maximum size and does not change further. After that, using the modelled root system, a simulation of moisture potential dynamic with known data on water inflow is carried out for the period, in which the pressures were measured, and based on this simulation, the value of the goal function ( 11 ) is calculated.

The initial data for this procedure are an estimate of evapotranspiration during the growing season; measurements of suction pressure during one or more irrigation cycles in the period of the maximum root system development and the highest water consumption; data on precipitation and irrigation rates during this period.

A problem for the minimization of the goal function ( 11 ) by fitting the parameters of the root system growth model is approximately solved by a genetic algorithm, in which potential solutions are coded by floating-point numbers in the form ( ps ,la ,lb ,ln ,bmax , , Lmax , ) . Maximal tree depth  max is here considered as given.

It is worth noting that, since the root system development model is stochastic, the goal function is random. Accordingly, the obtained solution to the minimization problem in the form of a set of parameters of the root system growth model means only that for such parameter values, the model during the period of the highest water consumption (if the obtained accuracy is higher than in the base case). This improvement is not guaranteed for a specific set of root growth model parameters, but is achieved for a certain root system shape obtained as a result of the modelling.

4. Simulation setting

The initial data for modelling were the measurements obtained using the iMetos micrometeostation, on which, in particular, Watermark suction pressure sensors were installed, from 31.05.2024 to 23.09.2024 during an agronomic experiment on growing corn for seed on an experimental field within the Makariv district of the Kyiv region of Ukraine.

Irrigation pipelines were placed at the depth of 30 cm directly in the rows of plants, the distance between which was 70 cm.

The layout of the Watermark sensors location within the soil massif is shown in Fig. 1.

The values of the parameters of the van Genuchten and Mualem models obtained after processing the data of laboratory studies of soil monoliths are given in Table 1.

For the purpose of calibrating the model ( 1 )suction pressures in the peri per 1 hour. Two irrigations were carried out during this period: on 04.07.2024 from 15:00 to 19:30, a total of 83 m3/ha, and on 06.07.2024 from 12:00 to 19:00, a total of 95 m3/ha.

Data irrigation on 13.07.2024 from 16:30 to 19:30, a total of 59 m3/ha, were used to check the scalability of the moisture transport model and to solve the optimization problem for selecting the parameters of the root system growth model. When solving the latter problem, the shape of the root system on 1.07.2024 was determined by modelling its growth starting from 31.05.2024, assuming that the average moisture content in the root zone is maintained in the range from 19.2% (corresponding to 20 kPa) to 25.3% (field capacity). modelling of water head dynamics according to the integer-order model ( =  = 1 ) in the period between 1.07.2024 9:00 and 17.07.2024 14:00 using root density function ( 6 ) without calibration; calibration of the integer-order model fitting the values of the empirical coefficients kET , k prec , kirr ,  and the test of model calibration of the integer-order model fitting the values of the empirical coefficients kET , k prec , kirr ,  , and filtration coefficients k f calibration of the fractional-order model fitting the values of the empirical coefficients kET , k prec , kirr ,  , filtration coefficients k f , and the orders  ,  of the fractional derivatives determination of the optimized values of parameters for the root system growth model along with the specific form of the root system using the results of calibration for the integer-order model (experiment 3); determination of the optimized values of parameters for the root system growth model along with the specific form of the root system using the results of calibration for the fractional-order model (experiment 4).

In all cases, the goal function had the form ( 11 ), and the grid size was 50x75 cells. Additional computational experiments with the goal function (12) showed that its use does not lead to a significant change in the accuracy of the solutions. Conducting experiment 3 in the heterogeneous case with the selection of different values of the filtration coefficients in different directions also failed to obtain solutions that were more accurate than in the homogeneous case k f ,x = k f ,z = k f . Also, ineffective were the experiments where, in addition to the filtration coefficients, we fitted the values of the residual soil moisture r and the parameters  , n of the van Genuchten model.

When calibrating the model using the PSO algorithm, the swarm contained 60 particles for the integer-order model and 100 particles for the fractional-order model. 50 iterations were performed. The values of the inertial weight  , cognitive coefficient  p and social coefficient  g were equal to 0.6. The regularization parameter in the goal function ( 11 ) was equal to 0.05. The search was carried out in the following ranges of parameter values:  , [0.8,1] , kET , k prec [0.01,10] , kirr [0,10] ,  [0, 2] , k f [10−9 ,5 10−6 ] (m / s) .

When determining the op shape by the genetic algorithm, the population size was equal to 20. 10 iterations were performed with a crossover probability of 0.75 and a mutation probability of 0.05. The search was carried out in the following ranges of parameter values: ps  5 10−8 ,5 10−6  (m / s) , la [0.05, 0.5] (m) , lb ,ln [0.01, 0.1] (m) , bmax [10,100] ,  [0,120] , Lmax [0.2,1] (m) ,  [0.1,10] . The maximum depth  max of the tree was equal to 3.

5. Simulation results

Simulation according to the above-described setting was performed using the software and input files accessible through https://github.com/sevaboh/root_growth . Computations were performed on the SCIT5 cluster of the Institute of Cybernetics of NAS of Ukraine.

Some of the obtained modelling results are given in Table 2. The measured and simulated water head dynamics in the case of the best overall modelling accuracy (integer-order model with root system shape fitting) are shown in Fig. 2 (sensor 1, the case of the largest absolute error), Fig. 3 (sensor 4 with low absolute modelling error), Fig. 4 (sensor 3 located directly near the emitter). The simulated water head field 1 hour after the irrigation was finished and the optimized shape of the moisture extraction function S are illustrated in Fig. 5. The results demonstrate that: • • • • • • the use of experimentally determined filtration coefficient values leads to fundamentally inadequate modelling results (errors >10 kPa for all sensors); the use of the fractional-order differential model, and the corresponding increase in the number of fitted parameters, allowed obtaining better modelling accuracy for the training dataset (period 1). However, for the testing dataset (period 2) the relative error F2 (v) increased in comparison with the integer-order model. This may result from overfitting, as Watermark sensor measurements represent averaged pressures over a vicinity rather than point-wise values. An additional argument in support of this assumption is the fact that in the case of the fractional-order model it was not possible to obtain an increase in modelling accuracy by fitting root system ; the smallest average absolute error for the testing dataset was achieved in experiment 5, which included the selection of root system shape for the integer-order moisture transport model. This effect was achieved primarily due to the smallest errors among all experiments in modelling the pressures at the depth of 5 cm (sensors 1 and 12). The soil in this depth was the driest with a slight influence of irrigation on moisture content. Accordingly, the pressures and absolute errors of their modelling were the largest; similar simulated water head distributions were obtained with significant differences in the values of the parameters of the model ( 1 )-( 10 ), which were selected in different ways. The reason behind this observation may be the different-directional effect of different factors on processes in the same zones. E.g., when the multiplier for evapotranspiration is increased together with the multiplier for the irrigation water flow and evapotranspiration is redistributed between plant transpiration and evaporation from soil surface, a similar water head distribution in the root-containing zone could be achieved; significantly different values of the fractional derivative orders were obtained for different directions. This indicates the possible heterogeneity of the structure of the studied soil; when compared to the measurements, in the case of all experiments conducted with the selection of filtration coefficients, a more intensive spreading of moistened front from the emitter was simulated (Fig. 3) along with a slower drying of the soil after the third irrigation (Fig. 4). The latter may indicate a change in atmospheric factors, which led to additional inaccuracies in the assessment of evapotranspiration compared to previous periods.

The generated optimized root system shapes for the cases of integer and fractional-order models are shown in Fig. 6. The found values of the root system growth model parameters in the case of using the integer-order moisture transfer model were equal to

( ps ,la ,lb ,ln ,bmax , , Lmax , ) = (2.98 10−6 , 0.24, 0.075, 0.068, 68, 79, 0.79, 6.63) , and, in the case of the fractional order model, to (1.88 10−6 , 0.18, 0.05, 0.065, 37, 61, 0.804, 7.59) .

Similar values of all parameters, except for the maximum number bmax of branches of the next tree level. This leads to the generation of similar root system shapes (Fig. 6) characterized by the predominant distribution breadthwise in the 0-30 cm layer, which occurs at the initial stages of development when the moisture from a deeply placed irrigation pipeline does not rise sufficiently to the upper soil layers. There is also an increased concentration of root mass in layers below 40 cm, where irrigation moisture is concentrated under the influence of gravity.

The difference in the density of the systems shown in Fig. 5 (a) and Fig. 5 (b), due to the normalization procedure used for obtaining the root system density function based on these models, does not have a significant impact on the modelling of moisture transport.

During the agronomic experiment, the weight of roots in 10 cm thick layers was measured. These measurements were used to assess the adequacy of the modelled shape of the root system to the actual one. Since the model does not take into account the different sizes of the root system segments, it was assumed that their weight is proportional only to their length. The accumulated part of the root system volume (for the model - by length, for the experimental data - by weight) in the soil layer, which spans from the surface to a certain depth, is shown in Fig. 7. The volume of the root system located in 30 cm thick layers is given in Table 3.

From the obtained data it is clear that, despite significant deviations when the step of 10 cm was used (Fig. 7), the volume of the root system generated according to the model in 30 cm thick layers adequately describes the experimental data (Table 3). Less smooth changes in the modelled root stimation of hydrotropism intensity. The volumes of the root system, which was generated using the fractional-order moisture transfer model, fully correspond to the experimental ones. When using the integer-order moisture transfer model, a larger than experimentally determined volume of the simulated root system is located in the 60-90 cm layer, but a smaller one in the 30-60 cm layer.

6. Conclusions

This study presents a novel approach to predict the dynamics of soil moisture availability to plants, which proposes to explicitly model the process of root system growth and take into account the potential effects of structural fractality of the soil considered as a continuous porous medium.

When modelling a dataset collected during corn cultivation under subsurface drip irrigation, the most accurate predictions were obtained when using the integer-order model together with the generated optimized root system shape.

Although the fractional-differential model showed the best results for the training dataset, its application to the testing dataset led to an increase in errors, which may indicate overfitting.

Modelling and optimising the root system improved the accuracy of moisture content modelling, especially for the upper soil layers, where, in the case of using the generalized root system density function, the errors were the largest. The simulated optimized shape of the root system demonstrates its preferential distribution breadthwise in the upper layers and concentration in the moistened zone under the emitter.

The proposed approach can be used to refine irrigation regimes and the parameters of irrigation pipelines placement along with the distance between emitter on a pipeline, which is especially important for the conditions of limited monitoring of soil moisture status.

Declaration on Generative AI

The authors have not employed any Generative AI tools.

References

[1] H. Hammouch , M. A. El-Yacoubi , H. Qin and H. Berbia , A Systematic Review and MetaAnalysis of Intelligent Irrigation Systems , IEEE Access 12 ( 2024 ) 128285 - 128304 . doi: 10 .1109/ACCESS. 2024 .3421322

[2] M. Kandelous , T. Kamai , J. V drip irrigation design and management parameters for alfalfa , Agricultural Water Management 109 ( 2012 ) 81 - 93 . doi: 10 .1016/j.agwat. 2012 . 02 .009

[3] G. Arbat , J. Puig-Bargués , J. Barragán , J. Bonany , F. Ramírez de Cartagena, Monitoring soil water status for micro-irrigation management versus modelling approach , Biosystems Engineering 100 ( 2 ) ( 2008 ) 286 - 296 . doi: 10 .1016/j.biosystemseng. 2008 . 02 .008

[4] L. Richards , Capillary conduction of liquids through porous media , Physics 1 ( 5 ) ( 1931 ) 318 333. doi: 10 .1063/1.1745010

[5] E. Holzapfel , M. Marino , A. Valenzuela , Drip irrigation nonlinear optimization model , J. Irrig. Drain Eng . 116 ( 4 ) ( 1990 ) 479 - 496 . doi: 10 .1061/(ASCE) 0733 - 9437 ( 1990 ) 116 : 4 ( 479 )

[6] V. Bohaienko , V. Bulavatsky , Fractional-Fractal Modeling of Filtration-Consolidation Processes in Saline Saturated Soils , Fractal Fract. 4 ( 2020 ) 59 . doi: 10 .3390/fractalfract4040059

[7] L. Umutoni , V. Samadi , Application of machine learning approaches in supporting irrigation decision making: A review , Agricultural Water Management 294 ( 2024 ) 108710 . doi: 10 .1016/j.agwat. 2024 .108710

[8] V. Clausnitzer , J.W. Hopmans , Simultaneous modeling of transient three-dimensional root growth and soil water flow , Plant and Soil 164 ( 1994 ) 299 - 314 . doi: 10 .1007/BF00010082

[9] D. Leitner , A. Schnepf , J. Vanderborght , From hydraulic root architecture models to efficient macroscopic sink terms including perirhizal resistance: quantifying accuracy and computational speed , Hydrol. Earth Syst. Sci . 29 ( 2025 ) 1759 1782. doi: 10 .5194/hess-29- 1759 - 2025

[10] G. Lei , W. Zeng, Th. H. Nguyen , J. Zeng , H. Chen , A. K. Srivastava , Th. Gaiser, J. Wu , J. Huang , Relating soil-root hydraulic resistance variation to stomatal regulation in soil-plant water transport modeling , Journal of Hydrology 617A ( 2023 ) 128879 . doi: 10 .1016/j.jhydrol. 2022 .128879

[11] T.H. Mai , A. Schnepf , H. Vereecken , Continuum multiscale model of root water and nutrient uptake from soil with explicit consideration of the 3D root architecture and the rhizosphere gradients , Plant Soil 439 ( 2019 ) 273 292. doi: 10 .1007/s11104-018-3890- 4 convergence issues during inverse identification of soil hydraulic properties , Advances in Engineering Software 173 ( 2022 ) 103278 . doi: 10 .1016/j.advengsoft. 2022 .103278