Optimizing Water Consumption in Richards’ Equation Framework with Step-Wise Root Water Uptake: A Simplified Model

This work arises from the need of exploring new features for modeling and optimizing water consumption in irrigation processes. In particular, we focus on water flow model in unsaturated soils, accounting also for a root water uptake term, which is assumed to be discontinuos in the state variable. We investigate the possibility of accomplishing such optimization by computing the steady solutions of a θ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document}-based Richards equation revised as equilibrium points of the ODEs system resulting from a numerical semi-dicretization in the space; after such semi-discretization, these equilibrium points are computed exactly as the solutions of a linear system of algebraic equations: the case in which the equilibrium lies on the threshold for the uptake term is of particular interest, since the system considerably simplifies. In this framework, the problem of minimizing the water waste below the root level is investigated. Numerical simulations are provided for representing the obtained results. Root water uptake is modelled in a Richards’ equation framework with a discontinuous sink term. After a proper semidiscretization in space, equilibrium points of the resulting nonlinear ODE system are computed exactly. The proposed approach simplifies a control problem for optimizing water consumption. Root water uptake is modelled in a Richards’ equation framework with a discontinuous sink term. After a proper semidiscretization in space, equilibrium points of the resulting nonlinear ODE system are computed exactly. The proposed approach simplifies a control problem for optimizing water consumption.


Introduction to the Physical Problem
Soil water content changes across the unsaturated zone because of inflows of water via infiltration from the surface, due, e.g., to rainfall and irrigation, and outflows caused by evaporation and root water uptake as well as percolation beneath the root zone. Richards' equation represents the most commonly used model for soil moisture dynamics into the unsaturated zone, approximated as a porous medium. Such a model combines the Darcy's law extended to unsaturated soils with mass conservation; it relies on the assumption that the gas phase does not need to be modeled, because its pressure is assumed to be constant (see Bear and Cheng 2010). The solution of this equation, both analytical and numerical, is complex because of the strong nonlinear relationships that link soil moisture and soil hydraulic conductivity to matric-potential.
In this paper, analogously to Roose and Fowler (2004), where the considered state variable is the relative water saturation, we are going to deal with the -based form of Richards' equation, as follows: being the state variable the volumetric water content , K and D state-dependent hydraulic functions, named hydraulic conductivity and soil-water diffusivity, respectively. In particular, represents the relationship between the soil water potential and the water content .
The advection diffusion equation (1) needs to be endowed with proper boundary conditions, for example of Dirichlet type: In this framework, the problem of describing the water uptake by roots arises, and it has been faced by different empirical and mathematical models. In analogy with Ohm's law for electrical current (Gardner 1991), the uptake by roots can be considered as a resistance to the fluid flow into the unsaturated zone; the same author still observes that "the uptake in the lower portion of actual root systems was less than one would expect from the number of roots found, even after allowing the gravitational factor"; therefore, "water uptake from the lower portion of the root system does not occur for some time and then only after significant quantities of water have been removed from the soil layers above". Still in Gardner (1991), the moving sink model is proposed: in practice, "uptake may be considered to take place in a relatively narrow zone, or sink, which moves down through the soil at a rate initially determined by the available water and the transpiration rate".
The starting point for modeling the root water uptake is Richards' equation, to whom a sink term needs to be added; in general, such sink term could depend explicitly on time (as in Mathur and Rao 1999) and on space (Mai et al. 2018). As in Broadbridge (1999), Difonzo et al. (2021), here it is assumed that plant root density does not depend explicitly on time; conversely, as in Jarvis (1989) we assume that the sink term depends not only on the state variable, but also on the depth (implying an implicit dependence on root length density), i.e., Still following Jarvis approach, we will assume that the sink term-acting, as in Jarvis (1989), by a stress index function-has different weights according to depth layers.
In order to simplify Eq.
(3), we assume that the Gardner's relation holds, as in Broadbridge et al. (2017), i.e., there exists a constant such that as a consequence Eq. (3) is replaced by The Gardner condition (4) is satisfied, for instance, for D and K chosen as in Table 1 of Broadbridge et al. (2017). Without loss of generality, we consider for instance the following conductivity function: which is an increasing function, being m a suitable positive constant, and K s the saturated hydraulic conductivity.
Let us now focus on the uptake term R. As summarized in Broadbridge et al. (2017), R( ) is usually very low and close to zero at low moisture contents near the wilting point where water uptake closes down. As soil moisture increases, R( ) strongly increases too, while it becomes almost constant when the water content approaches its maximum value. Accordingly, R( ) may be represented by a convex function ( R �� ( ) < 0 ) and generally by some sigmoid function with R �� = 0 near the wilting point.
Similarly to D' Abbicco et al. (2016), the sink term can be modeled as a Hill function, with the following form: being a certain threshold depending on the crop under consideration. As in D' Abbicco et al. (2016), if the order p of the Hill function is high enough (i.e., if the curve is steep enough at the threshold), the dynamics arising from (3) can be treated by a step function, whose numerical implementation is much more natural by proper event driven methods (see for instance Berardi 2014;Del Buono and Lopez 2015).
We point out that this approach is particularly effective for crops in which the root water uptake is function of the water content alone, such as horticultural crops (see Broadbridge et al. 2017, for instance), and does not take into account variations depending on the depth. In the following, we will merge these observations into the model (3).
For general and complex hydraulic functions, a numerical solution to Richards' equation is often the only way to obtain a solution.
The literature on numerical methods of Richards' equation is extremely rich, both in fundamental and in recent papers, and different issues arise in the numerical treatment of (4) K � = D; 1 3 such problem; it is worth citing, for instance, studies on time-stepping techniques (mainly low-order implicit methods, as implicit Euler, see Pop 2002;Kavetski et al. 2001;Keita et al. 2021), or different spatial discretization techniques such as mixed finite element or finite volume methods-endowed also with rigorous error estimates- (Arbogast et al. 1993(Arbogast et al. , 1996Eymard et al. 1999Eymard et al. , 2006Manzini and Ferraris 2004;Radu et al. 2004Radu et al. , 2008, discontinuous Galerkin (e.g., Clément et al. 2021;Li et al. 2007), gradient discretization schemes (Eymard et al. 2014), and virtual element methods (da Veiga et al. 2021). Equally challenging is the problem of solving the nonlinear system arising from the time integration: to this purpose, many schemes have been proposed, as L-scheme (List and Radu 2016;Mitra and Pop 2018), Newton and nested-Newton (Bergamaschi and Putti 1999;Casulli and Zanolli 2010), modified Picard (Celia et al. 1990), and a combination of Newton and Picard methods (Lehmann and Ackerer 1998). Nevertheless, analytical solutions in Richards' equation have been eagerly investigated (see, for instance: Leij et al. 2021;Broadbridge et al. 2017;Meunier et al. 2017;Basha 1999), both because they are an exact solution to a certain model, and because they provide a perfect benchmark to numerical methods.
In particular, the research of exact steady-state solutions for Richards equation has a long history and has nowadays a significant relevance. This study has been introduced in Warrick (1974), where the steady-state assumptions are justified as an approximation of high-frequency irrigation: a semi-infinite column or a saturated bottom boundary conditions are therein considered, together with different extraction functions. Such interest is confirmed in the milestone paper Philip (1989), motivated since, even for multidimensional unsaturated flows, the steady solutions are approached rapidly if close to the water source term, concluding that the application interest for such investigations is more practical than it might seem at first glance; on the other hand, the advantage of handling just linear systems is therein highlighted. More recently, in Basha (1999), the steady solution is derived using perturbation expansions methods, with a possibly nonzero water extraction term: as in other papers (see for instance Ursino 2000), the solution is represented as the superposition of the terms: the zero-order part, which is the linear solution, and the first-order correction term, which accounts for the nonlinearity. This concept is also reported and generalized in Severino and Tartakovsky (2014), claiming that "flow in both root systems and ambient soils can be represented by a sequence of steady states": also the aforementioned paper makes use of matched asymptotic expansions for steady flows, relaxing or eliminating classical previous assumptions (infinite soil column and constant hydraulic head prescribed on its surface). On the other hand, the study of steady-state solutions of the Richards equation is ancillary to the analysis of their stability; for instance, in Van Duijn et al. (2004), the authors study the stability of steady, vertically upward and downward flow of water in a homogeneous layer of soil, for several classes of soils.
The starting point for this paper, which considers step-wise root water uptake functions, is the semi-discretization of second-order partial derivative in (3), obtaining a system of differential inclusions; the study of steady-state solutions is therefore accomplished analytically simply by solving a linear system, under suitable linear constraints. Differently from previous papers on steady-state solutions, this study is aimed to an optimal choice of boundary conditions for minimizing irrigation and percolation beneath the root zone, while keeping at the same time the desired absorption at each depth under consideration: in practice, as detailed in the following, this work can be considered as a simplified optimal control approach for managing irrigation and water consumption. On the other hand, the presence of discontinuous thresholds allows to represent interesting properties of the eco-hydrological system, allowing to mimic a sort of self-regulation of the root system for controlling the uptake.
The paper is organized as follows. In Sect. 2, the model is presented: by applying the method of lines we move to study a system of nonlinear ODEs; in Sect. 3.1, the equilibrium points are described as the solution of an algebraic linear system and we discuss which conditions on the boundary values allow to reach a realistic equilibrium; in particular, in Sect. 3.2, we study the case in which water supply is due only to irrigation: here, the problem of minimizing the water waste below the root depth is also investigated. In Sect. 4, we specify the aforementioned discussion to the case of a 5-point mesh; finally, in Sect. 4.1, some simulations are presented.

The Proposed Model
In Jarvis (1989), an empirical model is used, which assumes that the soil depth is split into different layers, and water uptake is given as a function both of the potential transpiration rate and of a weighted stress index which accounts for the effects of the vertical distributions of roots and soil water content, for each layer. In Gardner (1991), the uptake process is represented by a model based on nonlinear behavior of the root membranes, and described by a distributed sink moving downward through the soil profile. Differently from , , Seus et al. (2018), Suk and Park (2019), here a homogeneous soil is considered for simplicity, albeit the root zone is generally made of overlapping soil layers with different hydraulic properties. On the other hand, as specified in Kuhlmann et al. (2012), heterogeneity of media could significantly affect the reliability of classical macroscopic models for root water uptake.
In the light of Gardner (1991), Jarvis (1989), here we propose a model in which the sink term R in (3) takes different form according to the depth: we suppose to consider a mesh defined by N + 1 points z i = iΔz for i = 0, … , N , and in each layer [z i−1 , z i ] we study Eq. (5) assuming that the sigmoid function R (defined as in (7), for instance) is replaced by a step-wise function defined on [ r , s ] , that is, for any z ∈ (z i−1 , z i ]: where ̄i ∈ [ r , s ] and i is a positive constant that is the greater the larger is the density of roots at layer [z i−1 , z i ] ; here, s denotes the saturated water content and r the residual water content, and both are hydraulic properties of the soil under consideration. On the other hand, we assume that there exists z ∈ [z N−2 , z N−1 ] such that for z >z there are no roots, and then the sink term becomes zero, that is, R ≡ 0 in [z N−1 , z N ] , and for any z ∈ (z N−2 , z N−1 ] , it holds where S N−1 is defined by (8) as done for the upper layers; here, for any i ∈ {1, … , N − 2} , by R(̄i, z) ∈ [0, i ] we mean that the sink term in ̄i is represented by a multivalued function which takes values in the interval [0, i ] ; in particular, for any time t and any depth z ∈ (z i−1 , z i ] with i ∈ {1, … , N − 1} such that (t, z) =̄i equation (5) reads as On the other hand, the discontinuity in (8) suggests to study equation (5) as a discontinuous differential equation using Filippov theory, in the frame of Agosti et al. (2016), Agosti et al. (2015), , where novel numerical techniques have been proposed; nevertheless, this treatise goes beyond the goals of our work. A sink term defined by means of a differential inclusion, similarly to (8), was already considered, for instance, in Kumar et al. (2014) where a dissolution-precipitation model is studied: therein, the sink term does not depend on the depth and a regularization technique is used in order to deal with the discontinuity. In our work, following the idea in Jarvis (1989), we suppose that the sink term depends explicitly on the depth, and then also the discontinuity point is not constant in the whole domain, but it is assumed to be constant and equal to ̄i in each interval (z i−1 , z i ] ; this can be justified, for example, by the variation of the root density distribution at different depths; however, the same approach works if one assume that the point of discontinuity remains the same at each depth, or even in the case in which it depends continuously on the depth. In principle, our approach could be applied to a 2D spatial domain, but explicit computations would become cumbersome, and therefore we do not consider it. Now, let us resort to the method of lines (MoL), the classical tool for numerically integrating PDEs including Richards' equation (Berardi and Vurro 2016;Tocci et al. 1997). Consider the standard semi-discretization for the second and the first derivative, i.e., Therefore, the right-hand side of (5) can be approximated as follows: In particular, by using the notations a − = (2 − Δz)∕4 and a + = (2 + Δz)∕4 , equation (11) reads By (10), we conclude that if i =̄i for some i ∈ {1, … , N − 2} , then (12) reads as (10) 0 ≤ This work is aimed at studying a simple approach for optimizing water consumption in irrigation. In particular, starting from the semi-discredized Richards' equation with a sink term (11), here the goal is to compensate vertical infiltration with root water uptake, in such a way that percolation does not occur beneath root zone, and ensuring at the same time sufficient uptake by roots.
Loosely speaking, the computation of equilibrium points for the system (11) (see Definition 1) is based on the idea that water moves and infiltrates through the layers yet keeping a balance between uptake and percolation. The solution behavior at the thresholds i is a key factor in such balance. To this purpose, we are interested in looking for the equilibrium points of such model. Indeed, the study of equilibrium points of a system paves the way to the analysis of turnpike behaviors in optimal control problems. In practice, the turnpike property describes the fact that trajectories of optimally controlled systems "most of the time" stay close to an equilibrium point. A deep review about the turnpike problem can be found in Grüne and Guglielmi (2018) where the authors investigate such properties for optimal control problems with linear dynamics and a cost function consisting of quadratic linear and constant terms.
The optimization of water consumption can be interpreted as the problem of minimizing the increase in water content N in the last point of the mesh. We will accomplish this goal by minimizing the flux in z N−1 .
In the last level [z N−1 , z N ] of the mesh, there is no root water uptake, and then by Eq.
(3), the Darcian flux in z N−1 can be written as: Due to the Gardner condition (4) and neglecting the argument for the sake of readability, the flux can be expressed as Then, after a second-order discretization, we can approximate By using this representation, we will investigate how to manage the irrigation in order to have the desired root water uptake at each level with the minimum value of q that corresponds to the minimum water consumption.

Research of the Equilibrium Points as Solution of an Algebraic Linear System
In this section, we assume constant boundary conditions 0 and Z for Eq. (3), and we study the corresponding steady solution = (t) ; by steady solution, we mean that ∕ t = 0 for t ∈ (t 0 , T) for suitable T > t 0 ≥ 0 . By the spatial discretization introduced in Sect. 2, such steady solutions can be studied as equilibrium points of the linear ODEs system defined by (11) according to the following definition: Definition 1 Let 0 , N ∈ [ r , s ] ; we say that the vector ( 1 , … , N−1 ) describes an equilibrium point for system (12) if the (N + 1)−tuple ( 0 , 1 , … , N−1 , N ) satisfies the algebraic equation: for any i = 1, … , N − 1 such that i ≠̄i and the following inequality holds for any i = 1, … , N − 1 such that i =̄i.
In Definition 1 the requirement (16) is an immediate consequence of (13). Due to our definition of the sink term R at each level [z i−1 , z i ] , it is clear by (11) that one should investigate the existence of equilibrium points in every possible configuration, each of which is distinguished from the other by the relative position of each i with respect to the corresponding threshold ̄i . More precisely, there are 3 N−2 possible configurations.
In the following, we consider in each possible configuration the set of indexes , m} , and we show that there exist a nonsingular diagonal block matrix A ∈ ℝ (N−1)×(N−1) and a vector b ∈ ℝ N−1 such that is an equilibrium point for system (11) if and only if the condition is satisfied for any i = 1, … , m ; in particular, in formula (17), the matrix A depends only on (i) parameter , (ii) mesh width Δz , whereas the vector b ∈ ℝ N−1 , depends on (i) parameter , (ii) mesh width Δz , (iii) boundary conditions 0 , N , iv) thresholds ̄k i for i = 1, … , m. From now on, we denote For technical reasons, we need to distinguish different cases.
1. Firstly, let us discuss the case m = 0 , that is, no water level i coincides with the corresponding threshold ̄i if the equilibrium is achieved. In such case, by (12), we deduce that the vector (u 1 , … , u N−1 ) has to satisfy the following equation: where A and b are, respectively, an (N − 1) × (N − 1) matrix and a vector of ℝ N−1 ; they have the following representation: and We remark that in this case, for any where A k j has the same representation as in (20) and Moreover, one obtains immediately Thus, we can introduce the block diagonal matrix that is and the vector b defined as 3. If k 1 = 1 and k m ≠ N − 2 , then we reduce to m systems; similarly to the previous case, for each j ∈ {2, … , m} , we have the following system of dimension where A k j and b k j have the same form as in (21); on the other hand, the vector Analogously to the previous case, we define the block diagonal matrix and the vector b m ∈ ℝ N−1 given by 4. If k 1 ≠ 1 and k m = N − 2 , then we reduce to m systems; in particular, the vector u 1 = (u 1 , … , u k 1 −1 ) satisfies the following equation (20) and On the other hand, for each j ∈ {2, … , m} , the vector u j = (u k j−1 +1 , … , u k j −1 ) satisfies the same Eq. (21). Finally, the identity (22) holds. Also in this case, we define the block diagonal matrix and the vector b m ∈ ℝ N−1 given by 5. If k 1 ≠ 1 and k m ≠ N − 2 , then we reduce to m + 1 systems; in particular, the vector u 1 = (u 1 , … , u k 1 −1 ) satisfies the same equation as in (26), and the vector u m+1 = (u k m +1 , … , u N−1 ) satisfies the same equation as in (24). Finally, for each j ∈ {2, … , m} , the vector u j = (u k j−1 +1 , … , u k j −1 ) satisfies Eq. (21).
Here, we define the block diagonal matrix and the vector b m ∈ ℝ N−1 given by In all the cases 1-5, our definition of A ∈ ℝ (N−1)×(N−1) and b ∈ ℝ N−1 allows to conclude that the vector u = (u 1 , … , u N−1 ) satisfies the equation Regardless of the configuration, the matrix A which we constructed is always the direct sum of Toeplitz tridiagonal matrices; thus, by Brugnano and Trigiante (1992), we know that a sufficient condition for the nonsingularity of A is that a + a − ≤ 1∕4 that is always verified for any choice of and Δz . As a consequence, Eq. (27) has always a unique solution that is Thus, according to Definition 1, the vector ( 1 , … , N−1 ) defined by is an equilibrium point at time t if and only if condition (16) is satisfied for any i ∈ {k 1 , … k m }.

Remark 1
We observe that the matrix (−A) −1 is a block diagonal matrix too, obtained as suitable direct sum of the matrices (−A k j ) −1 where the matrix A k j is defined in each case 2-5 in an appropriate way. In particular, for any j ∈ {1, … , m} , the inverse matrix (−A k j ) −1 can be computed explicitly by using a recurrence formula (see Usmani 1994). The same formula can be applied to evaluate (−A) −1 in case 1.

Remark 2
For a fixed choice of bottom boundary condition N and thresholds ̄i , the discontinuity of R in ̄1 allows to look for 0min and 0max in [ r , s ] such that for any 0 ∈ [ 0min , 0max ] at the equilibrium, it holds 1 =̄1 and condition (16) is satisfied for i = 1 . Then, it is possible to study the steady solutions for (3), for t ∈ (t 0 , T) even if the boundary condition 0 = 0 (t) is time-dependent with 0min ≤ 0 (t) ≤ 0max for any t ∈ (t 0 , T) . Indeed, since N remains constant, such a restriction on 0 (t) guarantees that for any t ∈ (t 0 , T) , it holds 1 (t) =̄1 , and then, even the water content at each level i does not depend on time. In fact, by formula (28), for any t ∈ (t 0 , T) , the equilibrium point ̃( t) = ( 1 (t), … , N−1 (t)) can be written as: since for any t ∈ (t 0 , T) , it holds 1 (t) =̄1 , then, the vector b(t) is not dependent on 0 (see formula (23) or (25) according to the value assumed by N−2 ) and thus, both the matrix A and the vector b(t) are constant with respect to t in the interval (t 0 , T) ; therefore, we can conclude that for any i ∈ {1, … , N − 1} , also i remains constant. In Sect. 4.1, an estimate of 0min and 0max is given for a certain choice of the parameters (see Fig. 2).

A Realistic Equilibrium
Once we have evaluated the equilibrium points i for our system (11), it is necessary to discuss for which choice of boundary conditions ( 0 , N ) and thresholds ̄i , such solution can really describe a realistic equilibrium in the chosen configuration. Indeed, for any fixed configuration, we have to ensure that the obtained equilibrium points are consistent with that configuration and with the soil physical properties. In practice, if in the selected configuration we are assuming S i ( i ) = 1 for some i ∈ {1, … , N − 1} , then we have to verify that our solution i satisfies the inequality i >̄i ; similarly, if S i ( i ) = 0 , we have to check the inequality i <̄i for our solution; finally, if i =̄i , we have to ensure that condition (16) is satisfied; finally, we have to require that each i is not smaller than the residual water content r and not greater than the saturated water content s . All the described conditions lead to impose suitable constraints on the thresholds ̄i ; once thresholds that meet these conditions have been selected, we identify a feasible region for 0 and N which guarantees the existence of an equilibrium point in the desired configuration; such feasible region is always described by a system of linear inequalities in the variables K( 0 ) and K( N ) . In the following, we discuss all these conditions in each of the cases 1-5 which we introduced in the previous section.
1. In case 1, none of the equilibrium points j coincides with the threshold ̄j . Since the conductivity function K is increasing, in order to reach a realistic equilibrium, we need to impose u j <ū j if j <̄j and, similarly u j >ū j if j >̄j , i.e., the following inequality has to be satisfied for any j = 1, … , N − 2 On the other hand, we have to guarantee that when the equilibrium is reached, the water content i is not greater than the saturated water content and not smaller than the residual water content that is equivalent to require for any j = 1, … , N − 1 , where u s = K( s ) and u r = K( r ) . By using (17), we can express each u j in terms of u 0 , u N ,ū 1 , …ū N−2 ; thus, (29) and (30) can be written, respectively, as (30) u r ≤ u j ≤ u s , and Once we have fixed all the thresholds ̄j , such inequalities can be read as constraints on 0 and N which allows to get the equilibrium in the desired configuration. 2. In case 2, we consider the set of indexes {k 1 , … , k m } such that k i =̄k i . Specifically, for each i ∈ {2, … , m − 1} , inequality (18) has to be satisfied that is the following inequalities hold: where Additionally, since k 1 = 1 and k m = N − 2 condition (18) implies and where and On the other hand, for any index j different from each k i , we have to keep in consideration the relative position of j with respect to ̄j . Loosely speaking, the following inequality has to be satisfied for any j ∉ {k 1 , … , k m } , j ≠ N − 1 ; by (17), this is equivalent to require for any i = 2, … , m , for any index j such that k i−1 < j < k i , where Finally, we have to ensure that u r ≤ u j ≤ u s , for any j = 1, … , N − 1 . By (17), this is equivalent to require and for any i ∈ {2, … , m} and k i−1 < j < k i , where c j was defined in (35), and finally 3. In case 3, the conditions which allow to get a realistic equilibrium can be obtained following the same ideas as in case 2. In particular, we need to fulfill the linear inequalities (31) and (32) related to the thresholds k 1 , … , k m , and the conditions (34), (36) for each k i−1 < j < k i to varying of i ∈ {1, … , m} . Additionally, condition (18) in the last threshold point k m can be rewritten as: where Finally, for each k m + 1 < j < N − 2 , the conditions and (30) correspond to the following conditions on the bottom water content N : where and being c j is the same defined in (40). 4. In case 4, the same threshold conditions (31) and (33) have to be satisfied, together with the linear inequalities (34) which describe the relation between the water content ̄j and the corresponding threshold j , and conditions (36), (37) which ensure that the water content remains not greater than the saturated water content and not smaller than the residual water content at each level [z j−1 , z j ] with j ≥ 2 . It remains to discuss suitable conditions for the water contents j for any 1 ≤ j ≤ k 1 . In particular, the condition (18) on the first threshold k 1 implies the following restriction on 0 : where On the other hand, the following inequalities have to be satisfied by u 0 : for any index j such that 1 ≤ j < k 1 , where and where c j was defined in (44). 5. In case 5, the equilibrium points describe a realistic equilibrium if conditions (31), (34), boundary values 0 and N satisfy the linear inequalities (39) and (41), related to the indexes k m ≤ j ≤ N − 1 , and the conditions (42), (43) and (45) related to the indexes 1 ≤ j ≤ k 1 .
As described so far, in each of the possible configurations, there are some necessary inequalities that the thresholds ̄1 , … ,̄N −2 have to satisfy; then, we can investigate which couples ( 0 , N ) fulfill the remaining system of inequalities. When not empty, the corresponding set of couples (u 0 , u N ) is a convex subset of [u r , u s ] × [u r , u s ] ; indeed, it represents the feasible region of a system of linear inequalities. On the other hand, in each configuration, it is possible to investigate if there exists a couple of boundary conditions ( 0 , N ) such that the flux q in z N−1 is non-positive; in fact, this means that the vertical infiltration is completely compensated by the root water uptake. By (14), we obtain that the condition q ≤ 0 is satisfied if Then, by the previous description, we note that if at least one of the thresholds ̄i is reached, then u N−1 and u N−2 can be expressed in terms of ū k m and u N , where k m is the index of the last reached threshold; then, the non-positivity of the flux is uniquely determined by the value of u N and it is independent on the value of 0 ; instead, if no threshold is reached, the condition q ≤ 0 introduces a new necessary relation between N and 0 . A more detailed analysis of the flux in z N−1 is presented in Sect. 4 (see Fig. 4).
Remark 3 In Sect. 3.1, we discussed which bounds on boundary conditions 0 , N and on thresholds ̄i are necessary in order to get a realistic equilibrium in each configuration. If we consider a mesh with a very large number N of points it would be more meaningful to write such estimates in terms of the densities, u * 0 ∶= u 0 ∕Δz , u * N ∶= u N ∕Δz and ū * i ∶=ū i ∕Δz , ( i = 1 … , N − 2 ); this allows in some cases to obtain necessary conditions on u * 0 , u * N and u * i which are independent on Δz , and then to get information on such bounds even in the asymptotic case, that is Δz → 0.
In order not to burden the presentation, we apply this idea to the simplest case in which at each mesh point z i the water content i coincides with the corresponding threshold ̄i when the equilibrium is reached; among other things this implies by (15) that u N−1 = a − u N + a +ūN−2 .
Following the introduction of Sect. 3.1, in this case, we have to ensure that for each i ∈ {1, … , N − 1} , the water content i is not smaller than the residual water content r and not greater than the saturated water content s ; additionally, for each i ∈ {1, … , N − 2} , the inequalities in (16) have to be satisfied; more explicitly, we have the following constraints on u 0 , u N and ū i : and, by (18) For simplicity, we suppose that the root depth is one meter and then Δz = 1 N m ; thus, recalling that a − = (2 − Δz)∕4 and a + = (2 + Δz)∕4 , all the given inequalities can be rewritten in terms of the densities u * 0 , u * N and ū * i , with i ∈ {1, … , N − 2} as where u * r ∶= u r ∕(Δz) and u * s ∶= u s ∕(Δz) ; additionally, by (16), we have the following constraints: Collecting all together the factors of the same powers of 1/N in each inequality, we rewrite that as analogously, as a consequence of (16), we have As N tends to infinity, the terms that determine the validity of each inequality are the factors of the lower power of 1/N which appears in the inequality; as a consequence, also in the asymptotic case Δz → 0 , that is N → ∞ , we can write a system of constraints on the densities u * 0 , u * N and ū * i which guarantees the existence of a realistic equilibrium in this configuration. Firstly, the inequalities in (46) have to be satisfied; moreover, (47) and (48) hold if and only if one of the following conditions is satisfied: 1 3 The fulfillment of conditions (49) and (50) is guaranteed by the following constraints: Finally, the inequality (51) is satisfied if and only if one of the following condition is satisfied: The same approach could be used in all the other configurations in order to get some information about the bounds on the density u * 0 in the asymptotic case Δz → 0 ; however, in order to avoid further technicalities, we omit this discussion.

Minimum Percolation Below the Root Depth with Water Supply Only by Irrigation
In this section, we want to investigate for which choice of boundary condition 0 at the surface it is possible to minimize the percolation below the root depth, assuming that the steady state is reached and that the water supply is obtained exclusively by irrigation. In this framework, some additional conditions have to be satisfied by our equilibrium points i obtained by (17) in order to guarantee that such equilibrium is realistic; indeed, for any i = 1, … N − 1 , the condition r ≤ i ≤ s should be replaced by the inequality r ≤ i ≤ 0 , that is equivalent to therefore, similarly to Sect. 3.1, such estimates can be written as a system of linear inequalities which involves the variables 0 , N ,̄1, … ,̄N −2 . In order not to burden the presentation, we do not explicitly present all these inequalities; we will study more in details such case in Sect. 4 considering a mesh of five points. We remark that if we assume that there is no capillary rise from soil level below the root zone, we cannot expect a non-positive flux q in z N−1 . Indeed, when a steady state is reached, we can presume that N−1 ≥ N ; instead, by (14), the non-positivity of the flux requires u N > u N−1 . Such conjecture is confirmed by the example treated in Sect. 4.1 (see (54)).

An Example
In this section, we apply all the ideas discussed in Sect. 3 to the case of a mesh of five points z 0 , … , z 4 . In particular, in each of the nine possible configurations, we firstly discuss explicitly which conditions on the boundary values 0 , 4 and on the thresholds ̄1 , ̄2 allow to get a realistic equilibrium, in the framework of Sect. 3.2; then, we will evaluate such equilibrium points. Finally, in Sect. 4.1, we will present some useful simulations in order to discuss for a fixed choice of thresholds ̄1 , ̄2 which choice of boundary conditions 0 and 4 allows to reach the equilibrium in the desired configuration with a minimum flux at the bottom of the root level. Let us suppose that the equilibrium is reached. Therefore, by using the notations introduced in Sect. 3, in the general case model, (11) takes the form: Thus, the hypothetical equilibrium points can be expressed as: If there exists i ∈ {1, 2, 3} such that i =̄i , the equilibrium points have a simpler representations with respect to (52); moreover, in each configuration, the flux can be expressed in terms of the boundary conditions and the thresholds ̄i . Let us discuss each case in detail.
Case 1: 1 and 2 are both below the threshold. In such case, it holds S 1 ( 1 ) = S 2 ( 2 ) = 0 . Thus, by (52), the equilibrium points are given by In order to get a realistic equilibrium, we require that i ≤ 0 for any i = 1, 2, 3, 4 ; moreover, we are supposing i <̄i at each depths. As a consequence, since the conductivity function K is increasing, the following inequalities have to be satisfied: (52) Case 2: 1 <̄1 and 2 >̄2 . In this case, we have S 1 ( 1 ) = 0 and S 2 ( 2 ) = 1 . Thus, by (52) we get the following representation for the equilibrium points: Here, a realistic equilibrium is reached if i ≤ 0 for any i = 1, 2, 3, 4 , 1 <̄1 and 2 >̄2 . Thus, due to the monotonicity of K, we require Case 3: i >̄i for both i = 1, 2 . In this case, the root water uptake is always not zero. Thus, if the equilibrium for the water content is reached, then it holds: Here, we require i ≤ 0 for any i = 1, 2, 3, 4 ; additionally, the conditions 1 <̄1 and 2 >̄2 need to be satisfied. This is equivalent to require the following system of inequalities: , Case 4: 1 >̄1 and 2 <̄2 . Here, we are assuming S 1 ( 1 ) = 1 and S 2 ( 2 ) = 0 . Thus, we get Here, we have the conditions 1 >̄1 and 2 <̄2 . Thus, we want Case 5: 1 <̄1 and 2 =̄2 . In this case, we are assuming that the water content in z 2 is known, and it is equal to the threshold ̄2 ; thus, we have S 1 ( 1 ) = 0 and S 2 ( 2 ) which can vary in the interval [0, 1]. Thus, we are able to represent the equilibrium points in a simpler way, that is .
As in the previous cases, we present which conditions have to be satisfied in order to reach a realistic equilibrium; in particular, since 2 =̄2 , the inequality (16) has to be satisfied for i = 2 ; summarizing, we have the following necessary conditions: Case 6: 1 >̄1 and 2 =̄2 . The water content in z 2 is known and it coincides with the threshold ̄2 ; moreover, the root-water uptake in the first level is maximum, equal to 1 . We can represent the equilibrium points as Analogously to the previous case, the following inequalities have to be satisfied: Case 7: 1 =̄1 and 2 =̄2 . Here, we are assuming that when the equilibrium is reached, the water content in the first and second levels of the mesh coincides with the corresponding threshold. Thus, the sink term in both the levels is defined by means of a differential inclusion. The equilibrium is given by We require that i ≤ 0 and that condition (16) is satisfied for i = 2 , that is: Case 8: 1 =̄1 and 2 <̄2 . In this case there is not root uptake in the second level z 2 ; on the other hand, the sink term in z 1 is defined by means of a differential inclusion. The equilibrium is given by A realistic equilibrium is obtained if the following conditions hold: Case 9: 1 =̄1 and 2 >̄2 . In this case, the root uptake in the second level z 2 is maximum, equal to 2 ; on the other hand, the sink term in z 1 is defined by means of a differential inclusion. The equilibrium is given by In order to guarantee that ( 1 , 2 , 3 ) describes a realistic equilibrium, the following conditions have to be satisfied:

Numerical Simulations
In this section, we present some MATLAB simulations which allow to understand more deeply how the water contents, at each level, are related to the corresponding boundary conditions. The model setup in (5) has been used, with hydraulic conductivity function (6) and parameters as in Broadbridge et al. (2017), i.e., m = 10 , r = 0 and s = 0.485 , and K s = 12 cm h −1 . Additionally, we set 1 = 35 ⋅ 10 −4 h −1 and 2 = 27 ⋅ 10 −4 h −1 and we fix the thresholds ̄1 = 0.34 and ̄2 = 0.29 ; finally, we assume = 1∕7 in (4) and Δz = 6cm.
In this framework, Fig. 1 represents which couples of boundary conditions ( 0 , 4 ) allow to reach a realistic equilibrium in the desired configuration, given that the linear inequalities presented in Sect. 4 have to be satisfied.
We can investigate how the water flux below the root zone changes, according to the chosen configuration. By (14) the flux in the third point can be written has As we explained in Sect. 3.2 if the water supply is due only to irrigation, then we cannot expect a non-positive flux below the root level, since this would require u 3 < u 4 . In fact, assuming for instance that 4 ≡ 0.2 we can evaluate in each configuration which is the minimum value of 0 which allows to reach a realistic equilibrium with minimum water flux at the bottom. We have the following correspondence: Then, we note that in all the possible configurations, one has positive flux; in particular, once the root water uptake is nonzero at least at one level, the flux in z N−1 is positive and the minimum flux is q ≈ 0.0111 cm h −1 . Additionally, by Eq. (54), we note that if the top boundary condition 0 varies arbitrarily in [0.3546, 0.4057), then the corresponding steady solution 1 coincides always with the first threshold ̄1 ; this property allows to apply our strategy not only for constant boundary conditions, but also for a time-dependent top boundary condition 0 = 0 (t) if 0min ≤ 0 (t) ≤ 0max for suitable values of 0min and 0max (see Remark 2). Let us suppose to know that the boundary condition at the bottom remains constant 4 = 0.2 ; moreover, we assume that the thresholds are ̄1 = 0.34 and ̄2 = 0.29 . Finally, we suppose that the boundary condition at the top is described by a periodic function; for instance where 0max and 0min represent, respectively, the maximum and the minimum value assumed by the boundary condition 0 . Then, Fig. 2 shows how the equilibrium water contents at each level of the mesh moves in time, choosing 0max = 0.4056 and 0min = 0.3546 as suggested in (54). On the other hand, Fig. 3 shows how the steady solutions i at each level i = 1, 2, 3 behave when 0 changes in [0, s ] for fixed 4 ≡ 0.2.
Let us assume now less constrains on the boundary conditions 0 and 4 ; in particular, at each level z i , we replace the assumption i ≤ 0 with the assumption i ≤ s , as in Sect. 3.1. Then, the representation of the feasible region in each configuration is described by Fig. 4.  We note that in this case, it is possible to investigate for which choice of boundary conditions it is possible to obtain no percolation below the root zone. Indeed, the black line in Fig. 4 represents the set of couples ( 0 , 4 ) such that the flux q in z N−1 is  (14); furthermore, the flux q is negative above such a curve, and positive below it.

Conclusions and Future Work
The introduction of a discontinuous root water uptake function allows us to investigate the steady-state solutions of Richards' equation (3), reduced to the study of an algebraic linear system. In particular, we show which conditions on thresholds and boundary values guarantee that the obtained equilibrium points are consistent with the hydraulic properties of the soil and with the chosen configuration; such conditions can be expressed as a system of linear constraints. The same approach is followed for case of water supplied only by irrigation, which needs additional constraints. In the latter case, a deeper study is presented for a space discretization made up by five points: we show explicitly the discussed constraints and we represent for each configuration the set of admissible boundary conditions; finally, assuming that the bottom water content is known and constant, we estimate the water content value allowing to minimize water percolation below roots level.
In future works, we plan to extend this treatise to more general hydraulic functions, and to make use of some tools from bifurcation theory, or even to apply this approach to nutrient uptake in unsaturated transport framework. Also, we are interested in pursuing some efforts for a deeper mathematical treatise of memory terms in the uptake function, recently investigated, for instance, in Wu et al. (2020). Furthermore, we will investigate the approximation of non-steady solutions to (3) in the frame of Filippov theory, making use also of recent VEM approaches (da Veiga et al. 2021). Finally, a challenge would be to integrate such physically based models into agricultural DSS as Zaza et al. (2018), Friedman et al. (2016). Funding This paper represents a first result of the RIUBSAL project, funded by Regione Puglia under the call "P.S.R. Puglia 2014/2020 -Misura 16 -Cooperazione -Sottomisura 16.2 "Sostegno a progetti pilota e allo sviluppo di nuovi prodotti, pratiche, processi e tecnologie". The first author also acknowledges that this work was partially supported by Gruppo Nazionale per il Calcolo Scientifico (GNCS-INdAM).

Availability of Data and Material (data transparency)
The datasets used and/or analysed during the current study are available from the corresponding author on request.

Code Availability
The MATLAB codes used for the numerical simulations are available from the corresponding author on request.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.