A preliminary model for optimal control of moisture content in unsaturated soils

In this paper we introduce an optimal control approach to Richards’ equation in an irrigation framework, aimed at minimizing water consumption while maximizing root water uptake. We first describe the physics of the nonlinear model under consideration, and then develop the first-order necessary optimality conditions of the associated boundary control problem. We show that our model provides a promising framework to support optimized irrigation strategies, thus facing water scarcity in irrigation. The characterization of the optimal control in terms of a suitable relation with the adjoint state of the optimality conditions is then used to develop numerical simulations on different hydrological settings, that support the analytical findings of the paper.


Introduction
More and more often, extreme weather events are accompanied by longer, more intense heat waves and consequent periods of drought, and a forecast global warming is increasing the urgent need of freshwater for human life.In this context, freshwater necessary for agriculture represents almost 70% of the whole amount of freshwater reserve [1].
In this scenario, a wise management of water resources for agricultural purposes is of fundamental importance, even at the irrigation district scale [2].Nevertheless, the vast majority of irrigation models just applies heuristic approaches for determining the amount and the timing of irrigation.Albeit an expert knowledge of agricultural and phenological issue is crucial, very seldom such information is coupled with proper mathematical models for controlling irrigation.In most sophisticated cases, the irrigation is managed mainly by studying soil water infiltration into the root zone, starting from Gardner's pioneering works, reviewed in [3]; several tools have been proposed in this context, generally providing some type of solver for Richards' equation, the advectiondiffusion equation which describes water infiltration in unsaturated porous media accounting also for root water uptake models.In the water resources management or agronomic framework several tools have been proposed in order to benefit from the Richards' equation for irrigation purposes: in [4] a Python code is presented, able to solve the Richards' equation with any type of root uptake model by a transverse method of lines; also, the popular Hydrus software is often used for simulating water flow and root uptake with different crops (e.g apples, as in [5], pecan trees as in [6]).
To the best of our knowledge, control methods are seldom applied to irrigation problems, and always with simplified models.For instance, in [1] a zone model predictive control is designed after defining a linear parameter varying model, aimed at maintaining the soil moisture in the root zone within a certain target interval; in [7] an optimal control is applied to an irrigation problem, modelled by a simpler (with respect to Richards' equation) hydrologic balance law; a simplified optimization method based on the computation of steady solutions of Richards' equation is proposed in [8]; finally, [9] presents an interesting sliding-mode control approach, but only considering a constant diffusion term in Richards' equation.An elegant approach for applying control techniques in a Richards' equation framework is provided in [10], yet with very different applications and tools, i.e. maximizing the amount of absorbed liquid by redistributing the materials, when designing the material properties of a diaper.
In this paper, we propose a model for solving an optimal control problem under the quasi-unsaturated assumptions (see Section 2), which provide a suitable hydrological setting that prevents to reach water moisture saturation in the soil.We derive the appropriate optimality conditions for the boundary control of a class of nonlinear Richards' equations, and implement these results in the development and computation of numerical solutions by a classical Projected Gradient Descent algorithm.Albeit the focus of this paper does not consist in proposing a novel numerical method, and here a standard MATLAB solver is integrated with control tools, some significant advances in the numerical solution of the unsaturated flow model deserve to be reminded; as a matter of fact, the numerics literature on Richards' equation is currently enriching and constantly evolving, since its possibly degenerate and highly nonlinear nature poses several challenges.For instance, the treatment of nonlinearities is a significant issue, and has been faced by different techniques, as Newton methods ( [11,12]), L-scheme or its variants [13,14,15] or Picard iterations [16].The discretization in space has been dealt, for instance, by finite elements or mixed finite element methods [17,18,19], discontinuous Galerkin [20,21], finite volume methods [22,23].A separate mention is deserved by the problem of infiltration in presence of discontinuities, which can be handled by domain decomposition methods [24,25], Filippov approach [26,27].For a more detailed discussion, the interested reader is addressed to the following complete reviews on numerical issues in Richards' equation [28,29].The aforementioned methods can be used to face peculiar issues in the numerical integration of Richards' equation, where MATLAB pdepe is known to blow up in a relatively small integration time (see, e.g., [4]).Thus, they could provide further directions in the development of specific algorithms for solving optimal control problems.
It is worth stressing that there is a vast literature about the problems of existence, uniqueness and regularity of the solution to degenerate parabolic differential equations.In the specific case of Richards' equations, the existence of solutions is tackled in the seminal paper by Alt and Luckhaus [30], whose ideas were subsequently developed by other authors.For example, [31] describes the semigroup approach to determine the existence of weak solutions, further developed in [32].In this paper we adopt the functional framework from [33], that exhaustively describes the minimal assumptions to derive the necessary regularity conditions to develop our analysis, for general classes of hydrological settings.We also refer to [34] for a thorough review of such results.
The paper is organized as follows: In Section 2 we introduce the quasi-unsaturated model of the Richards' equation; in Section 3 we describe the framework to ensure the well-posedness of the model of interest, and then derive first order necessary optimality conditions via the Lagrangian method in Section 4. Finally, we present numerical simulations in Section 5.

The mathematical model
Our work stems from the quasi-unsaturated Richards' model [33], describing a fast diffusion of water in soils.This framework provides a convenient setting to apply optimal control methods and derive optimal irrigation strategies.Indeed, from the mathematical point of view, the quasi-unsaturated diffusive model retains many crucial features of the nonlinear diffusion of water in the soil, while avoiding the special mathematical treatment required to face the limit case of a saturated diffusion.We consider the model of water diffusion in the space domain (0, Z), where Z > 0 is the depth of the domain under consideration.We denote by T ∈ (0, +∞) the time horizon of the interval (0, T ), by Q = (0, Z) × (0, T ) the space-time domain, and by (•, •) and • the scalar product and the norm in L 2 (Ω), respectively.In terms of hydraulic parameters, θ : Q → [θ r , θ S ) is the water content or moisture, where θ r and θ S represent the residual and the saturated water content, respectively.The function β : [θ r , θ S ) → [ , +∞) is the water diffusivity, satisfying the following condition: (H β ) β is locally Lipschitz continuous and monotonically increasing, β(θ) ≥ > 0 for all θ ∈ [θ r , θ S ), and lim The function β * is the primitive of the water diffusivity β that vanishes at 0. Thus, assumption (H β ) implies that β * is differentiable and monotonically increasing on [θ r , θ S ), and satisfies Moreover, the hydraulic conductivity K : [θ r , θ S ] → R is a non-negative, Lipschitz continuous on [θ r , θ S ], and monotonically increasing function.Other hydraulic functions of interest are the liquid pressure head h : Q → (−∞, 0), that is negative for unsaturated porous media, and the specific water capacity C(h) = dθ dh , that practically represents a storage term.The relation between the functions β, K and C is then expressed by β(θ(h)) := K(h) C(h) .With these notations, and assuming the vertical axis with downward positive orientation, the implicit form of the quasi-unsaturated model of the nonhysteretic infiltration of an incompressible fluid into an isotropic, homogeneous, unsaturated porous medium with a constant porosity and truncated diffusivity with non-homogeneous-Dirichlet Boundary Conditions (BCs) is given by the system (2.1) for t ∈ (0, T ) .
The source term f (θ) represents a sink function, that in our model describes the root water uptake.The non-homogeneous Dirichlet BCs are given by two functions v, g : (0, T ) → R, where the BCs at z = 0 is the control input v, that describes the irrigation strategy over the time horizon (0, T ), while the BCs at z = Z is a given function g modeling the interaction with the environment below the root zone Remark 2.1.Let us notice that the first equation in (2.1) is the Richards' equation.
In fact, we can easily compute that that provides the diffusion term in the classical mixed form of the Richards' equation [16, Eq. ( 3)].

Existence of solutions
Following [33, Section 4.5], we first reduce system (2.1) to a problem with homogeneous Dirichlet BCs.For this purpose, we assume the following hypothesis where ω t denotes the time derivative in the sense of distributions from (0, T ) to L 2 (0, Z).Thus, the function ω is defined on the cylinder Q and it attains the boundary values of v and g in z = 0 and z = Z, respectively.Then, introducing the function φ = θ − ω, system (2.1) is equivalent to where φ 0 := θ 0 − ω 0 and, for all φ ∈ V , We now introduce a suitable functional framework for (3.1).Let us consider the Gelfand triple H = L 2 (0, Z), V = H 1 0 (0, Z), and its dual V = H −1 (0, Z), with their usual norms.Then, (3.1) is equivalent to the abstract equation where the operator B(t) : V → V is defined by Notice that, thanks to (H ω ), we have that V ≤ M ω , and thus f B − ω t ∈ L 2 (0, T ; V ).In particular, this implies that the right-hand side of (3.2) is also in L 2 (0, T ; V ).
and, for a.e.t ∈ (0, T ) and for all ψ ∈ V , Existence of solutions to (3.2) with a source term independent of the water content θ is proved in [33,Proposition 5.3].For our purpose, we need to extend such wellposedness result to the case of a nonlinear source term f (θ), as in system (2.1).This can indeed be achieved by following a standard Galerkin approximation approach (see, e.g., [35,Lemma 5.3]).
Moreover, for appropriate initial conditions, we can prove that the solution stays away from the saturation value θ S uniformly in time.To this aim, we introduce the function j : R → (−∞, ∞] defined by and the space Theorem 3.3 (Theorem 5.7, [33]).Assume (H ω ), θ 0 ∈ M j , and that f : [θ r , θ S ] → R is Lipschitz continuous.Moreover, assume that f is non-negative, that is, there exists f m ∈ [0, ∞) such that f m ≤ f , and where Then the solution θ to problem (2.1) satisfies

The optimal control problem
In this section we formally derive the first order necessary optimality conditions for the cost functional where u = v − θ r , v is the control that appears in (2.1), λ > 0 is the coefficient of the control cost, f : [θ r , θ S ] → R describes the normalized root water uptake model as in (2.1), and θ is the solution to (2.1) with f as the sink term.Roughly speaking, the performance index (4.1)optimizes the root water uptake (see, for example, expression (5.2) in Section 5, where f is maximized when f ≡ 1) while minimizing the irrigation cost u.
In this setting, it is natural to consider the following space of admissible control (4.2) Fixing g ∈ L 2 (0, T ), θ 0 ∈ L 2 (0, Z), we introduce the control-to-state operator Λ : ) solution of (2.1).Theorem 3.2 ensures that the mapping Λ is well-posed.We can thus reformulate the minimization of a functional J(θ, u) constrained to the control system (2.1) in terms of the so-called reduced cost functional J : U ad → R defined by J(u) := J(Λ(u), u).We first introduce the Lagrangian functional where p = (p, p 1 , p 2 ) are adjoint variables that will be useful to find a representation of the optimal control.After integration by parts, we can rewrite the Lagrangian functional as Hereafter, we shall assume that the source term f ∈ H 1 (θ r , θ S ) to justify the following computations.In order to derive the first order optimality conditions of problem (4.1)-(2.1)with input constraints (4.2), we enforce the condition D θ L(θ * , u * , p * )θ = 0 for all θ, that determines the equation satisfied by the adjoint variable p; and the condition D u L(θ * , u * , p * ) • (u − u * ) ≥ 0 for all u ∈ U ad , that returns the optimality condition satisfied by any optimal control u * .After direct computations, we get that Thus, we deduce that the adjoint variable p satisfies where On the other hand, since ≥ 0 for all u ∈ U ad .We thus obtain that any optimal solution (θ * , u * , p * ) of problem (2.1)-(4.1)-(4.2) must satisfy the optimality system for t ∈ (0, T ), (4.3a) for all u ∈ U ad , where we recall that v = θ r + u.In the next section, we exploit this optimality system to build suitable algorithms to numerically solve the optimal control problem (4.1)-(2.1).

Algorithm and numerical simulations
Our optimization procedure will follow the Projected Gradient Descent (PGD) described in Algorithm 1 (see [35] for a thorough introduction to such optimization algorithms).However, when solving (4.3) with PGD, it could happen that Theorem 3.3 is not satisfied at each iteration, thus incurring numerical difficulties due to the singularity of water diffusivity β at θ = θ S .Therefore, we shall approximate it by truncation: given a small ε > 0, we define (5.1) as shown in the Figure 1.Regularization (5.1) is a standard technique when dealing with Richards' equation to handle singularities in the diffusion term, and it is used in finite difference schemes [30,36] or FEM [37] for both the mathematical and numerical analysis of degenerate, and possibly doubly-degenerate, parabolic equations.In the following simulations, we are then actually computing the numerical solutions to (4.3) after replacing β with β ε as defined in (5.1).
Moreover, we have set a maximum number of iterations equal to 100 before exiting PGD iterations, a tolerance of 10 −5 , and a regularization parameter ε = 10 −3 for water diffusivity in (5.1).We stress that the order of magnitude of ε has been chosen so to be consistent with that of the different θ S values selected in all the simulations that follow.Algorithm 1 Projected Gradient Descent Algorithm applied to Richards' equation for determining optimal irrigation.Intermediate steps for solving direct (line 3) and adjoint (line 4) problems and for computing the optimal descent direction step-size (line 6) are performed using MATLAB pdepe and fmincon functions, respectively.p n ← adjointProblem(θ n ) Solving the adjoint problem 5: s n ← arg min s J(pr end if 13: end while Moreover, we select a root water uptake model of Feddes type (used, for instance, in [4,6]) as source term in (2.1).Its expression is given by with the following values, in cm: h 4 ≈ −820, h 3 ≈ −400, h 2 ≈ −350, h 1 = 0. Also, we set ϕ = 0.1/Z, where Z is the soil depth, and λ = 0.1 in (4.1).
Remark 5.1.Let us notice that the maximum value for f (h) in (5.2) is set to 1 for normalization purposes.In fact, when it comes to practical problems, one uses f (h) properly rescaled according to experimental evidences through the factor ϕ, which is the ratio of the potential transpiration rate and the rooting depth, as explained in [38].
Moreover, we stress that in general one does not necessarily require the source term f (h) to be zero for the values of h corresponding to the boundary of [θ r , θ S ].However, from a physical point of view, it makes sense for a source term to vanish when the soil is either dry or saturated.This is exactly the case of Feddes-type source terms as the one we consider in (5.2) for our numerical simulations.
In Example 5.2 and Example 5.3 below, we simulate a soil described by Haverkamp model [16,39], whose constitutive relations are given by representing water retention curve and hydraulic conductivity, respectively.We first verify that Haverkamp model falls within the quasi-unsaturated model for a suitable choice of the parameters involved in the model setting.In fact, in Haverkamp model we have ( with β 1 , β 2 > 0, that is Lipschitz, monotonically increasing and bounded from below.In order to satisfy assumption (H β ), we shall ensure that lim From (5.4), a straightforward computation provides that this condition is satisfied if and only if β 2 > 1.
We have performed a simulation of T = 3 hours with a maximum depth of Z = 70 cm.
As in (2.1), boundary condition at the top varies in time according to the irrigation strategy: (5.6) θ(0, t) = θ top (t) = u(t) + θ r , while bottom condition has been chosen so to be constant over time: (5.7) θ(Z, t) = θ bottom (t) = 0.9θ r + 0.1θ S , t ∈ [0, T ].Finally, initial condition is linearly varying over time as Our simulations have been produced using MATLAB active-set algorithm; we report that same results are obtained using MATLAB sqp.We have observed that convergence is reached after 3 iterates within the given tolerance, and the numerical solution is locally optimal.Results are in Figure 2. As can be seen, the optimal control framework succeeds in determining an optimal control that optimizes the performance index (4.1),with a reduced water consumption and average water content over time.
Example 5.3.In this second simulation, using the same soil as in Example 5.2, we consider a time-varying bottom condition where θ b1 := 0.9θ r + 0.1θ S , θ b2 := 0.7θ r + 0.3θ S , while top condition is given by Moreover, initial condition is It turns out that MATLAB sqp converges, in 3 iterates, to a local optimal solution, further providing the best results if compared to active-set.Results are displayed in Figure 3.
In Example 5.4 and Example 5.5 that follow, we consider the classical Van Genuchten-Mualem constitutive relations in the unsaturated zone, given by In order to verify under which conditions Van Genuchten-Mualem model satisfies the quasi-unsaturated model, we need to analyze its corresponding function β(θ(h)).Letting ϕ(h) := 1 1 + |αh| n , from (5.12) and exploiting the fact that h < 0, it follows that (5.13) It is an easy computation that β is Lipschitz, monotonically increasing and bounded from below.In order to satisfy assumption (H β ), there needs lim From (5.13), this condition is satisfied if and only if n > 1.This is the case for the simulations reported below.More specifically, we are going to consider a Berino loamy fine sand and a Glendale clay loam, with parameters drawn from [40,26].Here, as in Example 5.3, we consider a time-varying bottom condition where θ b1 := 0.3θ r + 0.7θ S , θ b2 := 0.1θ r + 0.9θ S ; boundary condition at the top of the domain is, again as in Example 5.3, However, now initial condition is a quadratic polynomial function of depth and is set as Results relative to this soil are depicted in Figure 4 and are obtained using MATLAB active-set, where Z = 50 cm and T = 12 hours.
For this experiment, we fix Z = 30 cm and T = 36 hours; bottom boundary condition is given by where θ b1 := 0.5θ r + 0.5θ S , θ b2 := 0.7θ r + 0. Results are depicted in Figure 5. Here, we employed to MATLAB sqp for solving the optimization problem by PGD.

Conclusions
In this paper we introduce an optimal control approach aimed at optimizing the water content provided by irrigation, applying Richards' equation for unsaturated flow.We make use of quasi-unsaturated model introduced in [33], extending the well-posedness results for nonlinear sink terms and deriving suitable optimality conditions for an irrigation performance index of tracking type.We set the model within a MATLAB solver by implementing a properly adapted Projected Gradient Descent method, and provide significant numerical results over a meaningful variety of soils; a deeper analytical treatise of the control system is beyond the scopes of this paper, and it is currently under investigations by the authors.This paper could pave the way to an extensive use of control techniques for optimizing irrigation in real life applications, and this framework could easily be incorporated in existing irrigation software based on Richards' equation solvers.Moreover, it is worth investigating qualitative features of the more general saturated-unsaturated model, for which there is an increasing need of both numerical and analytical results and approaches.In this context, tools from set-valued analysis and discrete control techniques could carry improvements in understanding such problems.

( a )
Optimal water content in the spatio-temporal domain.(b) Optimal water content profiles over time.(c)Average optimal water content and optimal control.

Figure 2 .
Figure 2. Numerical simulations relative to Example 5.2, where initial condition is given by (5.8) and boundary conditions are as in (5.6) and (5.7), respectively.

( a )
Optimal water content in the spatio-temporal domain.(b)Average optimal water content over time.(c)Optimal water content profiles and optimal control.

Figure 3 .
Figure 3. Numerical simulations relative to Example 5.3.For this simulation, initial condition is given by (5.11) and boundary conditions are given by (5.9) and (5.10), respectively.

( a )
Optimal water content in the spatio-temporal domain.(b)Average optimal water content over time.(c)Optimal water content profiles and optimal control.

Figure 4 .
Figure 4. Numerical simulations relative to Berino loamy fine sand in Example 5.4, where we set bottom and top boundary conditions as in (5.14) and (5.15), respectively, whilst cubic initial condition is as in (5.16).

( a )
Optimal water content in the spatio-temporal domain.(b)Average optimal water content over time.(c)Optimal water content profiles and optimal control.

Figure 5 .
Figure 5. Numerical simulations relative to Glendale clay loam in Example 5.5.Bottom condition is given by (5.17), top condition by (5.18); initial condition is a quadratic polynomial as defined in(5.19).
Figure 5. Numerical simulations relative to Glendale clay loam in Example 5.5.Bottom condition is given by (5.17), top condition by (5.18); initial condition is a quadratic polynomial as defined in(5.19).