Fuzzy Analytical Solution of Horizontal Diffusion Equation into the Vadose Zone

: The process of how soil moisture proﬁles evolve into the soil and reach the root zone could be estimated by solving the appropriate strong nonlinear Richards’ equation. The nonlinearity of the equation occurs because diffusivity D is generally an exponential function of water content. In this work, the boundary conditions of the physical problem are considered fuzzy for various reasons (e.g., machine impression, human errors, etc.), and the overall problem is encountered with a new approximate fuzzy analytical solution, leading to a system of crisp boundary value problems. According to the results, the proposed fuzzy analytical solution is in close agreement with Philip’s semi-analytical method, which is used as a reference solution, after testing 12 different types of soils. Additionally, possibility theory is applied, enabling the decision-makers to take meaningful actions and gain knowledge of various soil and hydraulic properties (e.g., sorptivity, inﬁltration, etc.) for rational and productive engineering studies (e.g., irrigation systems).


Introduction
Infiltration is a common physical phenomenon of water movement in porous media that is of great interest in many earth and plant sciences.Soil water flow plays an important role in understanding the phenomena of runoff, groundwater recharge, and the transport of contaminants.Especially in the vadose zone, soil moisture strongly influences plants' growing processes, and for irrigation problems, it is important to know how the soil moisture profiles evolve into the soil and reach the root zone.Historically, Buckingham (1907) [1], Gardner and Widsoe (1921) [2], and Richards (1931) [3] are the pioneers in the development of soil water movement, introducing the ideas of capillary potential, capillary conductivity, and diffusion phenomena in the concept of soil water movement, which were completed later by Childs (1936aChilds ( , 1936b) ) [4,5].Klute (1952) [6] presented the equation of flow arising from Darcy's (1856) law [7], while the law of conservation of mass is as follows: where θ = the moisture content (cm 3 /cm 3 ), K = the unsaturated hydraulic conductivity (cm/s), Φ = the total potential (cm) Φ = Ψ − z, Ψ = the pressure potential or capillary potential (cm), and z = the gravitational component (cm).In the case of one-dimensional horizontal flow or in cases where gravity may be neglected (fine-textured soils), in which the moisture gradient's influence is much more important than gravity, the phenomenon is described by the following equation: ) Hydrology 2023, 10, 107.https://doi.org/10.3390/hydrology10050107https://www.mdpi.com/journal/hydrologywhere D = diffusivity (L 2 /T).According to Philip (1957Philip ( , 1969) ) [8,9], this equation represents the water movement in a horizontal column and is called the absorption equation because it describes the wetting up of the column under tension.Many methods exist with their corresponding advantages and disadvantages, such as semi-analytic methods [9], finite difference [10][11][12][13][14][15], finite element methods [16,17], finite control volume methods [18,19], flux concentration methods [20], and approximate analytical methods [21][22][23][24][25]. Henceforth, our research is focused on approximate analytical methods.Tolikas et al. (1984) [26] presented a simple analytical solution to the problem of horizontal absorption under the assumption that diffusivity is an exponential function of soil water content.Brutsaert (1982) [27] presented a group of 14 solutions with different diffusivities for the problem of desorption, considering, as pointed out by Crank (1975) [28], that sorption and desorption are related by the correspondence principle.Specifically, if the solution for sorption with a given D(θ) is θ = f(x,t) in normalized form, then the function f1 = 1 − f(x,t) is the solution for the desorption problem with the diffusivity D(1 − θ).Lisle and Parlange (1993) [29] used a large class of transformations in conjunction with symmetry properties to obtain an analytical reduction of the transformed equation for certain diffusivities, emphasizing that the symmetry leaves the equation invariant.Hang and Zhiqiang (1997) [30] presented an approximated analytical solution with the Boltzmann transformation method for the case of an exponential diffusivity form and a profile of soil water content of a finite extent under the assumption that the wetting front is known.Prevedello et al. (2008) [31] derived a new analytical solution of the Boltzmann transformation ϕ as a function of matrix potential for horizontal water infiltration into a sand soil sample without invoking the concept or use of D(θ).The derivation assumes that a similarity exists between the soil water retention function and the Boltzmann transformation ϕ.Tzimopoulos et al. (2015) [32] presented a new approximate analytical solution of the nonlinear diffusion equation, and the solution is given for the transformed equation through the Boltzmann transformation.It is considering the case of an exponential function of the diffusivity under the following assumptions: (a) the profiles of soil water content have a finite extent; (b) the concentration at the boundaries is constant; and (c) the reduced flux of Philip (1973) [20] is of a special form [33].More recent works have implemented the diffusion equation in applications for volatile organic compound transport [34], mass transfer [35], and Turing pattern formation [36].
The calculation of water flow in the unsaturated zone requires exact knowledge of the initial and boundary conditions and the various soil parameters, and this assumption is based principally on in situ measurements.In general, this assumption is subject to different kinds of uncertainty due to human and machine imprecision.In many cases, the uncertainties were considered in statistical terms as random variables with given mean values, variances, and correlations.However, these methods require exact knowledge of mean values, variances, and correlations and often suffer from an insufficient amount of accurate in situ measurement data.Zadeh (1965) [37] introduced the theory of fuzzy sets and fuzzy logic and covered all these kinds of uncertainty.A significant number of research studies presenting uncertainties were carried out in that field, especially regarding the fuzzy differentiation of functions.Puri and Ralescu (1983) [38] generalized and extended Hukuhara's fundamental study [39] (H-derivative) of a set of values appearing in fuzzy sets.Vorobiev and Seikkala (2002) [40], O'Regan et al. (2003) [41], and Nieto et al. (2006) [42] have worked in the theoretical and applied research fields on fuzzy differential equations with an H-derivative, but their method has presented certain drawbacks since it has led to solutions with increasing support along with increasing time.To overcome this drawback, the generalized derivative gH (gH-derivative) was introduced by [43][44][45][46].The gH-derivative will be used henceforth for a more extensive degree of fuzzy functions than the Hukuhara derivative.In addition, a recent review work by Mazandarani and Xiu (2021) [47] presents a chronological survey on fuzzy differential equations of integer and fractional orders, with the corresponding future perspectives and challenges to be discussed.
In the present article, the problem of fuzzy, unsteady water movement in the vadose zone in the horizontal direction is examined and described by a nonlinear partial diffusion equation with fuzzy conditions.In this context, we developed and propose a new fuzzy analytical solution that has never been presented before in the literature, coupled with the innovative possibility theory (see Section 2.2.1).Specifically, the nonlinear equation was first fuzzified and subsequently translated to a system of second-order crisp boundary value problems called the corresponding system for the fuzzy problem.Subsequently, using the Boltzmann transformation, the crisp problem was transformed to a system of two classical ordinary differential equations [32].With the help of α-cuts, profiles of moisture content θ versus Boltzmann transformation ϕ were estimated, as were the storativity (S) and cumulative moisture content (I).In addition, for evaluation purposes, the present analytical solution and the approximated solution of Hang and Zhiqiang (1997) [30] were tested using the semi-analytical solution of Philip as a reference solution.

Classical (Crisp) Case
The one-dimensional horizontal nonlinear flow of water in unsaturated soil with initial and boundary conditions for absorption is described by the following equation [6]: under constant head condition, with initial and boundary conditions as follows: • Boundary condition where θ s = the saturated soil moisture content at x = 0 and θ r = the residual soil moisture content.The diffusivity D (θ) can be closely expressed by an exponential function [48][49][50] as follows: The non-dimensional variables are now introduced in Equation ( 1), where L = a characteristic length usually equal to the wetting front, Dr = D(θ r ) = the diffusivity at the residual soil moisture content, τ = the non-dimensional time, D* = the non-dimensional diffusivity, and λ 1 = β (θ s − θ r ).The advantage of introducing nondimensional variables is that they allow the generalization of the theory, speeding up the computational processes.Then, Equation (1) with initial and boundary conditions (2) and (3) becomes the following: Subjacent to the following: • Boundary condition The Boltzmann transformation ϕ* = Xτ −1/2 is now introduced in Equation ( 6), and due to conditions (7) and ( 8), Equation ( 6) is reduced to the following ordinary differential equation: The conditions (2) and ( 3) are transformed into the following: where ϕ * f is the wetting front.Philip (1957) [8] introduced the term sorptivity in order to define the capacity of the medium to absorb or desorb water by capillarity.He expresses the term as follows: The non-dimensional Boltzmann transformation ϕ* and the non-dimensional sorptivity S* are related to the dimensional Boltzmann transformation ϕ and the dimensional sorptivity S by the following expressions: Tzimopoulos et al. (2015) [32] proposed the following: (a) an explicit analytical solution, which has the Boltzmann transformation as the dependent variable and the soil water moisture as the independent variable; (b) a simple form for the sorptivity.Both solutions are presented in normalized form.The solutions were tested with 12 soils and showed a close agreement with Philip's semi-analytical method [9], which has been considered a reference method.The solution concerning sorptivity is given as follows: or in dimensional form: The solution concerning Boltzmann transformation is given as follows: or in dimensional form: where E i (λ 1 ), E i (λ 1 Θ) are the exponential integrals [51,52].Equation ( 16) satisfies the boundary conditions (10).Indeed, for the following: Hydrology 2023, 10, 107 According to Ramanujan et al. (2000) [53], the function E i (x) − γ − lnx, γ ≈ 0.5772156649 . . .for small values of x is included in the following interval: If we choose x = 5 × 10 −5 as the minimum value, the above relation gives a value of E i (λ 1 Θ) λ 1 Θ = 5 × 10 −5 = −8.32631.With this value of integral, the value of ϕ is equal to the moisture content front (ϕ front ), which is as follows: Now the boundary condition ϕ * → ∞, Θ → 0, becomes the following: From the above expressions ( 14) and ( 16), it is clear that by knowing the soil moisture contents θ S , θ r , and the diffusivity as an exponential function of Θ, it is possible to find the following:

•
The moisture content front ϕ f ront ;
Table 1 shows the parameter values for the 12 different soils under consideration and a comparison of the wetting fronts and sorptivities between the values obtained by the analytical values in this paper and Philip's semi-analytical solution, which is in close agreement with the experimental values and was taken as a reference profile [32].The relation easily follows the dimensional form of I where I = cumulative absorption, L (L 3 /L 3 ).From this one, it follows the flow velocity at x = 0:

Fuzzy Case
In order to facilitate an understanding for readers unfamiliar with the fuzzy and possibility theories, we describe definitions concerning some preliminaries in the fuzzy and possibility theories and definitions about the differentiability of fuzzy numbers.

Fuzzy Theory
which is called membership function.The value µ U (x), x ∈ X is the membership value of x ∈ X and expresses the membership degree of x ∈ X , i.e., the truth of fuzzy logical proposition: The membership function µ U (x), defined also as U(x), has the following properties: Definition 2. Let X be a Banach space and U be a fuzzy set on X.We define the α-cuts of U as , and for α = 0, we define the following closure: Definition 3. Let Ҡ(X) be the family of all nonempty compact convex subsets of a Banach space.
The space of all compact and convex fuzzy sets on X is denoted as Hydrology 2023, 10, x FOR PEER REVIEW 6 of 17

Fuzzy Case
In order to facilitate an understanding for readers unfamiliar with the fuzzy and possibility theories, we describe definitions concerning some preliminaries in the fuzzy and possibility theories and definitions about the differentiability of fuzzy numbers.

Fuzzy Theory
which is called membership function.The value  (),  ∈  is the membership value of  ∈  and expresses the membership degree of  ∈  , i.e., the truth of fuzzy logical proposition: The membership function  (), defined also as  (), has the following properties: Definition 2. Let X be a Banach space and  be a fuzzy set on X.We define the α-cuts of  as [ ()] = { ∈   () ≥ },  ∈ (0,1], and for α = 0, we define the following closure:

Definition 3. Let Ҡ(X) be the family of all nonempty compact convex subsets of a Banach space. A fuzzy set 𝑈 on X is called compact if
The space of all compact and convex fuzzy sets on X is denoted as Ƒ (X).[61,62], the membership function and the α-cut form of a fuzzy number  are equivalent, and in particular, the α-cuts [ ] = [ ,  ] uniquely represent  , provided that the two functions are monotonic ( increasing,  decreasing) and  ≤    = 1.For every  ,  ∈ Ƒ (R), the metric structure is given by the Hausdorff distance as follows: Notation:

Fuzzy Case
In order to facilitate an understanding for readers unfamiliar with the fuzzy and possibility theories, we describe definitions concerning some preliminaries in the fuzzy and possibility theories and definitions about the differentiability of fuzzy numbers.Definition 2. Let X be a Banach space and  be a fuzzy set on X.We define the α-cuts of  as [ ()] = { ∈   () ≥ },  ∈ (0,1], and for α = 0, we define the following closure:

Definition 3. Let Ҡ(X) be the family of all nonempty compact convex subsets of a Banach space. A fuzzy set 𝑈 on X is called compact if
The space of all compact and convex fuzzy sets on X is denoted as Ƒ (X).[61,62], the membership function and the α-cut form of a fuzzy number  are equivalent, and in particular, the α-cuts [ ] = [ ,  ] uniquely represent  , provided that the two functions are monotonic ( increasing,  decreasing) and  ≤    = 1.For every  ,  ∈ Ƒ (R), the metric structure is given by the Hausdorff distance as follows:

Definition 4. Let [𝑈 ] ∈ Ƒ (R). The α-cuts of 𝑈 are as follows: [𝑈 ] = [𝑈 , 𝑈 ]. According to the representation theorem of
According to the representation theorem of [61,62], the membership function and the α-cut form of a fuzzy number U are equivalent, and in particular, the α-cuts

Fuzzy Case
In order to facilitate an understanding for readers unfamiliar with the fuzzy and possibility theories, we describe definitions concerning some preliminaries in the fuzzy and possibility theories and definitions about the differentiability of fuzzy numbers.Definition 2. Let X be a Banach space and  be a fuzzy set on X.We define the α-cuts of  as [ ()] = { ∈   () ≥ },  ∈ (0,1], and for α = 0, we define the following closure:

Definition 3. Let Ҡ(X) be the family of all nonempty compact convex subsets of a Banach space. A fuzzy set 𝑈 on X is called compact if
The space of all compact and convex fuzzy sets on X is denoted as Ƒ (X).

Definition 4. Let [𝑈 ] ∈ Ƒ (R). The α-cuts of 𝑈 are as follows:
[ ] = [ ,  ].According to the representation theorem of [61,62], the membership function and the α-cut form of a fuzzy number  are equivalent, and in particular, the α-cuts [ ] = [ ,  ] uniquely represent  , provided that the two functions are monotonic ( increasing,  decreasing) and  ≤    = 1.For every  ,  ∈ Ƒ (R), the metric structure is given by the Hausdorff distance as follows:

Fuzzy Case
In order to facilitate an understanding for readers unfamiliar with the fuzzy and possibility theories, we describe definitions concerning some preliminaries in the fuzzy and possibility theories and definitions about the differentiability of fuzzy numbers.Definition 2. Let X be a Banach space and  be a fuzzy set on X.We define the α-cuts of  as [ ()] = { ∈   () ≥ },  ∈ (0,1], and for α = 0, we define the following closure:

Fuzzy Theory
Definition 3. Let Ҡ(X) be the family of all nonempty compact convex subsets of a Banach space.A fuzzy set  on X is called compact if [ ] ∈Ҡ(X), ∀ ∈ [0,1].The space of all compact and convex fuzzy sets on X is denoted as Ƒ (X).

Definition 4. Let [𝑈 ] ∈ Ƒ (R). The α-cuts of 𝑈 are as follows:
[ ] = [ ,  ].According to the representation theorem of [61,62], the membership function and the α-cut form of a fuzzy number  are equivalent, and in particular, the α-cuts [ ] = [ ,  ] uniquely represent  , provided that the two functions are monotonic ( increasing,  decreasing) and  ≤    = 1.For every  ,  ∈ Ƒ (R), the metric structure is given by the Hausdorff distance as follows:

Fuzzy Case
In order to facilitate an understanding for readers unfamiliar with the fuzzy and possibility theories, we describe definitions concerning some preliminaries in the fuzzy and possibility theories and definitions about the differentiability of fuzzy numbers.Definition 2. Let X be a Banach space and  be a fuzzy set on X.We define the α-cuts of  as [ ()] = { ∈   () ≥ },  ∈ (0,1], and for α = 0, we define the following closure: [61,62], the membership function and the α-cut form of a fuzzy number  are equivalent, and in particular, the α-cuts [ ] = [ ,  ] uniquely represent  , provided that the two functions are monotonic ( increasing,  decreasing) and  ≤    = 1.For every  ,  ∈ Ƒ (R), the metric structure is given by the Hausdorff distance as follows:

Fuzzy Case
In order to facilitate an understanding for readers unfamiliar with the sibility theories, we describe definitions concerning some preliminaries i possibility theories and definitions about the differentiability of fuzzy nu 2.2.1.Fuzzy Theory Definition 1.A fuzzy set  on a universe set X is a mapping  :  → [0,1] which is called membership function.The value  (),  ∈  is the membersh and expresses the membership degree of  ∈  , i.e., the truth of fuzzy log ( | = ) =  ().

Definition 4. Let [𝑈 ] ∈ Ƒ (R). The α-cuts of 𝑈 are as follows:
[ ] = [ ,  the representation theorem of [61,62], the membership function and the α-cut form ber  are equivalent, and in particular, the α-cuts [ ] = [ ,  ] uniquely vided that the two functions are monotonic ( increasing,  decreasing) and  1.For every  ,  ∈ Ƒ (R), the metric structure is given by the Hausdorff distan . Suppose that the functions U − α and U + α are real-valued functions, differentiable w.r.t.(with respect to) x, uniformly w.r.t.α ∈ [0, 1].Then the function U(x) is gH-differentiable at a fixed x ∈ [a, b] only if one of the following two cases holds: (a) (U − α ) (x) is increasing, (U + α ) (x) is decreasing as functions of α, and U − 1 is increasing as functions of α, and U ∂x .In both of the above cases, the U α (x) derivative is a fuzzy number.Definition 2. Let X be a Banach space and  be a fuzzy set on X.We define the α-[ ()] = { ∈   () ≥ },  ∈ (0,1], and for α = 0, we define the following closur Definition 3. Let Ҡ(X) be the family of all nonempty compact convex subsets of a B A fuzzy set  on X is called compact if [ ] ∈Ҡ(X), ∀ ∈ [0,1].The space of all c convex fuzzy sets on X is denoted as Ƒ (X).[61,62], the membership function and the α-cut form of a ber  are equivalent, and in particular, the α-cuts [ ] = [ ,  ] uniquely repre vided that the two functions are monotonic ( increasing,  decreasing) and  ≤ 1.For every  ,  ∈ Ƒ (R), the metric structure is given by the Hausdorff distance as  .In both of the above cases derivative is a fuzzy number.

Definition 3. Let Ҡ(X) be the family of all nonempty compact convex subsets of a Banach s A fuzzy set 𝑈 on X is called compact if [𝑈 ] ∈Ҡ(X), ∀𝛼 ∈ [0,1]. The space of all compact convex fuzzy sets on X is denoted as Ƒ (X). Definition 4. Let [𝑈 ] ∈ Ƒ (R). The α-cuts of 𝑈 are as follows: [𝑈 ] = [𝑈 , 𝑈 ].
Accordin the representation theorem of [61,62], the membership function and the α-cut form of a fuzzy n ber  are equivalent, and in particular, the α-cuts [ ] = [ ,  ] uniquely represent  , vided that the two functions are monotonic ( increasing,  decreasing) and  ≤   1.For every  ,  ∈ Ƒ (R), the metric structure is given by the Hausdorff distance as follows: .In both of the above cases, the  derivative is a fuzzy number.Definition 6. gH-Differentiable at x0.Let  : [, ] → Ƒ (R) and  ∈ [, ] with  ()  () both being differentiable at x0.We say that [63]: are differentiable real-valued functions with respect to x, uniformly for α ∈ [0, 1], then f (x) is g-differentiable and we have the following [63]:

Fuzzy Case
In order to facilitate an understanding for readers unfamiliar with the fuzzy and possibility theories, we describe definitions concerning some preliminaries in the fuzzy and possibility theories and definitions about the differentiability of fuzzy numbers.Let X be a Banach space and  be a fuzzy set on X.We define the α-cuts of  as [ ()] = { ∈   () ≥ },  ∈ (0,1], and for α = 0, we define the following closure: The space of all compact and convex fuzzy sets on X is denoted as Ƒ (X).[61,62], the membership function and the α-cut form of a fuzzy number  are equivalent, and in particular, the α-cuts [ ] = [ ,  ] uniquely represent  , provided that the two functions are monotonic ( increasing,  decreasing) and  ≤    = 1.For every  ,  ∈ Ƒ (R), the metric structure is given by the Hausdorff distance as follows: .In both of the above cases, the  () derivative is a fuzzy number.Definition 6. gH-Differentiable at x0.Let  : [, ] → Ƒ (R) and  ∈ [, ] with  () and  () both being differentiable at x0.We say that [63]: , (x 0 , t 0 ) ∈ D and U − α (x 0 , t 0 ), U + α (x 0 , t 0 ) are real-valued functions and partially differentiable w.r.t.x.We say that [46,64]: x at (x 0 , t 0 ) if: The same is valid for

Fuzzy Case
In order to facilitate an understanding for readers unfamiliar with the fuzzy and possibility theories, we describe definitions concerning some preliminaries in the fuzzy and possibility theories and definitions about the differentiability of fuzzy numbers.[61,62], the membership function and the α-cut form of a fuzzy number  are equivalent, and in particular, the α-cuts [ ] = [ ,  ] uniquely represent  , provided that the two functions are monotonic ( increasing,  decreasing) and  ≤    = 1.For every  ,  ∈ Ƒ (R), the metric structure is given by the Hausdorff distance as follows: {| −  |, | −  |} Definition 5. gH-Differentiability [63] is the same, then • [(ii)-p]-differentiable w.r.t.x, if the type of [gH-p]-differentiability of both U α (x 0 , t 0 ) and ∂ U α (x 0 ,t 0 ) ∂x p.gH is different, then

Possibility Theory
Definition 11.A possibility measure Π on a set X (e.g., a set of reels) is characterized by a possibility distribution π : X → [0, 1] and is defined by the following: ∀A ⊆ X, Π(A) = sup{π(x), x ∈ A}.
Definition 12.A degree of necessity Ness X on a set X (e.g., a set of reals) is characterized by the non-possibility (one minus possibility) of A complement (A C ).
Definition 14.Two possibility distributions π X , π X are consistent with the probability distribution p X .The π X distribution is more specific than π X if it is π X < π X .
A possibility distribution π * x consistent with the probability distribution p X is called maximal specificity if it is more specific that each other possibility distribution π X :π * x (x) < π X (x), ∀x.
).In this case, the fuzzy number Y * is the fuzzy estimator of Y and verifies the following: so that the probability of the possibility α-cut is greater than 1 − α.

Fuzzy Model
We write Equation (1) in its fuzzy form according to [69]: with the following: • Boundary conditions The fuzzy parameters were considered to be the following: These parameters are considered symmetrical triangular fuzzy numbers and are illustrated in Figure 1.We can find solutions to the fuzzy problem of Equation ( 25) and the initial and boundary conditions ( 26) and ( 27) by utilizing the theory of [22,46,63,64,70] and translating the above fuzzy problem to a system of second-order crisp boundary value problems, hereafter called the corresponding system for the fuzzy problem.Therefore, many crisp BVP systems are possible for the fuzzy problem, but in the current work, the solution is restricted to the following systems (Cases A and B), which describe exactly the physical problem.
The following fuzzy expression is valid and is referring to α-cuts: In this case, for more convenience, the above expression could be written as follows, related to fuzzy partial differential equations: The main difference between the above two cases is that the operator (−) describes the left side of the α-cut while the (+) describes the right side of the α-cut, constructing in this way a possible range in which the desired value will be.
Case A For Case A, the following expression arises: Following the same procedure described above for the classical case, the non-dimensional and dimensional forms are presented below:

•
Non-dimensional form: 2 ) • Dimensional form: Case B For Case B, the following expression arises: Following the same procedure described above for the classical case, the non-dimensional and dimensional forms are presented below:

Construction of the Hang and Zhiqiang Solution
For comparison purposes, the Hang and Zhiqiang (1997) [30] solution was selected mainly due to its simplicity in applications.However, this solution is valid under the following two certain conditions: 1.
Diffusivity is an exponential function of Θ; 2.
The value of the moisture wetting front is a priori known.
The Hang and Zhiqiang analytical solution initially had the following form: ln[e βθ 0 − ( A(e βθ 0 − e βθ r ) − 0.5ϕ 2 In order to make a comparison with Equation ( 6) and with Philip's reference solution, the above equation is transformed into the following expression: The numerical values for ϕ f ront were taken from Equation (20).In both cases, these values are the mean of the values ϕ θ α=1 , θ − α=0 , and θ + α=0 .Figure 4 is considered very informative for drawing conclusions since it illustrates the membership functions of the profiles (Figure 4a,b) and the sorptivity S (Figure 4c,d) for both samples (1, 2), also giving space for the application of the possibility theory and the meaningful interpretation of the fuzzy results.In that regard and by carefully observing the different cases in Figure 4, we can conclude that, as an example, at α-cut = 0.05, we have valuable information regarding the θ and S intervals of confidence for all cases, which, according to the possibility theory, could provide us with confidence for a probability greater than 95%.For example, in Figure 4b, the interval [0.394, 0.438] denotes that the desired value will be in this range with a probability of >95%.To compare the results produced by the three different methods and for both samples 1 and 2, the mean relative error was chosen as a metric.

Results
Table 2 presents the differences in an analytical way to validate the strong agreement of our proposed solution with Philip's reference solution.Specifically, regarding sample 1, the difference is ε mean = 3.7 × 10 −3 , while for sample 2, the difference is ε mean = 3.7 × 10 −3 .As it emerged from the visualized results, the deviation between the Hang and Zhiqiang with Philip's reference solution is confirmed here as well, where for samples 1 and 2, we have ε mean = 2.5 × 10 −2 and ε mean = 7.8 × 10 −2 , respectively.

Conclusions
The proposed fuzzy solution of Richard's equation is based on the Hukuhara theory [39] with the generalized Hukuhara (gH-derivative), as well as its extension to partial differential equations.Today, practical problems useful for engineers can be solved, considering the fuzziness of various sizes.
In the current work, real measurements of 12 different types of soils were used in order to draw reliable conclusions.The application of our proposed solution to real cases of soils has proved that the present analytical solution is in close agreement with Philip's semi-analytical method, which is used as a reference solution that proves the accuracy and reliability of the new fuzzy analytical method.In addition, it was observed to be extremely easy and simple to calculate in comparison to other methods, without affecting the accuracy of the results.
Furthermore, it should be highlighted that according to the possibility theory, for practical cases, such as irrigation, drainage, and water resources projects, the engineers and designers could now have a better understanding of the real physical conditions (including the uncertainties), knowing the confidence intervals of the crisp value of sorptivity and cumulative infiltration with a certain strong probability.

Definition 15 .Definition 16 .
For a number Y with a known and continuous probability distribution function p, the fuzzy number Y, which has a possibility measure Π( Y) = µ Y is the fuzzy estimator of Y and has an α-cut of Π Y (α) = Y[α].This fuzzy number satisfies the consistency principle and verifies P( Y[α]) = Ness Y [α] = 1 − α, so that the probability of the possibility α-cut is equal to 1 − α.The α-cuts Y[α] are the confidence intervals of P, and the confidence level is α.Conjecture [68].For a function Y = Y (X 1 , X 2 , .... X n ) with unknown probability distribution function, a fuzzy number may be constructed Y * = Y( X * 1 , X * 2 , . . .X * n ), and the α-cut is equal to the following: Y (a) The soil moisture contents θ s (saturated soil moisture content); (b) θ r (the residual soil moisture content); (c) The parameter λ 1 (in the exponent of the diffusivity D).

Figures 2 and 3
Figures 2 and 3 illustrate the profiles of this study versus (vs.) the profiles of the Hang and Zhiqiang solution and Philip's semi-analytical method (the reference method).It is found that the mean relative errors ε mean between the present analytical method and the Hang and Zhiqiang solution are as follows: (a) Soil sample 1 ε mean = 2.669 × 10 −2 ; (b) Soil sample 2 ε mean = 8.04 × 10 −2 .In both cases, these values are the mean of the values ϕ θ α=1 , θ − α=0 , and θ + α=0 .

Figure 2 .
Figure 2. Profiles of moisture content θ versus Boltzmann transformation ϕ in the case of soil sample 1.The red color corresponds to θ α=1 , the blue color to θ − α=0 , and the green color to θ + α=0 .

Figure 3 .
Figure 3. Profiles of moisture content θ versus Boltzmann transformation ϕ in the case of soil sample 2. The red color corresponds to θ α=1 , the blue color to θ − α=0 , and the green color to θ + α=0 .

Figure 4 .
Figure 4. Membership function of (a) moisture content θ at ϕ = 1.2 cm/h 1/2 in the case of soil sample 1; (b) moisture content θ at ϕ = 6 cm/h 1/2 in the case of soil sample 2; (c) sorptivity in the case of soil sample 1; and (d) sorptivity in the case of soil sample 2.

Figure 6 .
Figure 6.Membership function of cumulative absorption I.