1 Introduction

The Richards’ equation, which quantitatively describes the process of fluid infiltrating the porous medium such as soil, concrete, etc., is a nonlinear partial differential equation (PDE) about the saturation of unsaturated porous media in most cases^{[1, 2, 3, 4]}. It is given as

where S is the dimensionless saturation of unsaturated porous media, x and t are the space and time coordinates, K(S) and D(S) are the hydraulic conductivity and diffusivity function [L²/T] of the porous media, respectively. Because the material properties of the porous medium are different, in the present study, several models, such as the Fredlund model^[5], Campbell model^[6], Brooks-Corey model^[7], Gardner model^[8], van Genuchten model^[9], Russo model^[10], etc., have been suggested for determining the basic parameters K(S) and D(S). Among them, the Brooks-Corey model, which had a precise definition and was adopted in Ref. ^[11], is the more commonly used model.

Several analytical and numerical solutions for the Richards’ equation with the Brooks-Corey model were obtained by many researchers, such as Bruce and Klute^[12], Parlange^[13], Prasad and Romkens^[14], Witelski^[15], Lockington et al.^[3], Barry et al.^[16], and Omidvar et al.^[11]. In these studies, some have obtained the implicit solutions of Richards’ equation without the hydraulic conductivity K(S)^{[1, 3, 12, 13, 14]}, and others based on wetting front propagation, have derived the traveling wave approximate solution^[15]. Recently, Omidvar et al.^[11] used the differential transform method (DTM) and the homotopy perturbation method (HPM) to obtain the analytical solution of the Brooks-Corey model with 1-2 orders power series types.

Learning from the studies above, the DTM is found to be robust in finding solutions of Richards’ equation. This paper tries to use the DTM to search the analytical solution of Brooks-Corey model with the high order power series type, which is usually used to describe the unsaturated fluid flow in concrete or soil^{[1, 17]}.t

Nevertheless, there are two problems in the DTM. One is that the convergence speed is very slow as the calculated point is away from the convergence center. For example, some- times it needs about 400 approximations to ensure accuracy in calculating the Brooks-Corey model^[11]. The other is that the DTM, which is based on Taylor’s theorem, is usually hard to solve the discontinuous solution problem, such as discontinuity between initial and boundary conditions^[3].

For the problems of the DTM above, the solution of extension of Heaslet-Alksne technique^[1] is used as an intermediate variable to accelerate the convergence speed of the DTM. For dealing with the discontinuous solution problem, the Boltzmann variable φ=$\frac{x}{{\sqrt t }}$ and a tiny time point △t are introduced in two forms K(S) = 0, and K(S)≠0, respectively, to change the Richards’ equation into an ordinary differential equation (ODE) instead of using the wave variable from Alquran^[18] in this paper. In the end, two numerical examples are given to show the accuracy of the presented approach.

2 Brooks-Corey model and definite conditions

According to Parlange et al.^[1] and Witelski^[4], the functional relationship between D(S), K(S), and S is presented as

where D₀, K₀, m, and n are the empirically-fitted constants.

From the works of Lockington et al.^[3] and Witelski^[4], the initial condition of Eq. (1) could be expressed as a constant moisture concentration

and the boundary conditions are given as in which S₀ is the initial moisture concentration, and S_L (> S₀) is the boundary value. We would consider the one-dimensional movement of water into the initially dry porous medium, which indicates that S₀ is usually a tiny constant in this paper. Here is a discontinuous solution problem about Eq. (1) with Eqs. (4)-(6).

3 DTM with intermediate variable

The DTM was presented in 1986 by Zhou^[19], who used Taylor series for the solution of linear and nonlinear differential equations in the form of polynomial. From the pioneer studies on the DTM conducted by several researchers^{[20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]}, the basic concept of the DTM is that each function S(x) in the governing equation can be expressed as

where n is the derivative order, and The properties of the DTM are given as follows.

Theorem 1 If S(x) = S₁(x) ± S₂(x), then Y (n) = Y₁(n) ± Y₂(n).

Theorem 2 If S(x) = cS₁(x), then Y (n) = cY₁(n), where c is a constant.

Theorem 3 If S(x) = $\frac{{{d^k}{S_1}(x)}}{{d{x^k}}}$, then Y (x) = $\frac{{(n + k)!}}{{n!}}$ Y₁(n + k).

Theorem 4 If S(x) = S₁(x)S₂(x), then Y (n) =$\sum\limits_{{n_1} = 0}^n {} $ Y₁(n1)Y₂(n - n1).

According to Alquran^[18], an intermediate variable ξ is introduced to change the PDE into an ODE. Then, the DTM can be applied to the resulting ODE through the following steps.

(i) Construct an intermediate variable ξ = ξ(x, t) such that

(ii) Change the differential governing equation

into

(iii)φ(ξ) in Eq. (9) could be expressed as

where U(n) is the differential transform of φ(ξ).

(iv) Apply the DTM to Eq. (11), and calculate U(n) as

(v) The approximate solution is

4 DTM with intermediate variable for Richards’ equation

There are two forms for Richards’ equation, namely, the horizontal infiltration form K(S) ≠ 0 and the vertical infiltration form K(S)≠0 in Eq. (1). In this section, we shall discuss them, respectively.

Case 1 K(S) = 0

Due to K(S) = 0, Eq. (1) can be expressed

With the Boltzmann variable φ=$\frac{x}{{\sqrt t }}$ , Eq. (15) can be transformed into^[3] Learning from Bruce and Klute^[12], we integrate Eq. (16) and yield in which When the porous medium is nearly dry, which suggests that S₀ is very close to 0, we have the following approximation of Eq. (17): Substituting the Boltzmann variable φ into the initial condition (4), we obtain Similarly, the boundary conditions (5) and (6) are Combining the same results of Eq. (20) and Eq. (22), we have new boundary conditions

From Eqs. (20)-(24), we can find that the discontinuous problem between the initial condi- tion (4) and the boundary condition (5) can be solved by the Boltzmann variable φ, and the discontinuous solution problem (1) with initial and boundary conditions (4)-(6) is changed into the continuous solution problem (19) with boundary conditions (23) and (24). According to the extension of the Heaslet-Alksne technique^[1] and the step (i) in Section 3, we construct an intermediate variable ξ as follows:

Then, the Boltzmann variable φ can be expressed as in which U(k) ($\frac{{1!}}{{k!}}\frac{{{d^k}\phi }}{{d{\xi ^k}}}$|S=S_L) is the differential transformation of φ.

The new boundary conditions (23) and (24) may be changed into

By introducing the intermediate variable ξ, Eq. (19) is

and boundary conditions (27) and (28) are From Eq. (31), we could find that φ → +∞ is satisfied naturally when S = S₀ in Eq. (26). According to Theorem 4 in Section 3, the transformed version of Eq. (29) is where W(k - 1), Y₁(i), and Y₂(k - 1 - i) (k = 1, 2, 3, · · · ) are the differential transformations of (S_Lⁿ - ξ),$\frac{{d\phi }}{{d\xi }}$, and ∫ξ _SLⁿ-S₀ⁿ $\frac{\phi }{{{{(S_L^n - \xi )}^{\frac{{n - 1}}{n}}}}}$dξ, respectively, and can be expressed as where Y₁(0) and Y₂(0) are initial terms of Y₁(i) and Y₂(k -1-i), respectively, which are given as follows:

Substituting Eq. (26) into Eq. (30), we obtain

Based on the initial terms Y₁(0) and Y₂(0) and the differential transformation formula (38), we could obtain U(k) in Eq. (26) by calculating the parametersW(k-1), Y₁(i), and Y₂(k-1-i) in Eqs. (32)-(35) term by term.

Case 2 K(S) ≠ 0

Because K(S) ≠ 0, the discontinuous problem between the initial condition (4) and the boundary condition (5) can not be solved by the Boltzmann variable φ. At this time, we still suppose that the saturation S could be expressed as

where S_r is the 1 order approximate solution of Case 1 from Eqs. (25)-(38), which could be given as and G(k)(=$\frac{{1!}}{{k!}}\frac{{{d^k}S}}{{dS_r^k}}$|S_r=S_L) is the differential transformation of S.

Instead of solving the differential transformations of the initial condition (4) and the bound- ary conditions (5) and (6) directly, we choose a tiny time point △t → 0, substitute t = △t into Eq. (1), and yield

which is the differential equation about the space variable x. By introducing the intermediate variable S_r in Eq. (40), Eq. (41) can be transformed into where the calculated parameters λ₁ and λ₂ are

According to Eqs. (39) and (40), when t → 0 or x → +∞, we have S_r = S₀ and S = S₀, which indicates that the initial condition (4) and the boundary condition (6) are naturally satisfied. On the other hand, when x = 0, we have S_r = S_L, which means that the boundary condition (5) is

According to Theorem 1 in Section 3, the transformation form of K(S)≠0 in Eq. (42) is

in which According to Theorem 4 in Section 3, W₂(k - 1), W₃(k - 1), and W₄(k - 1) are where where the initial terms F₁(0), F₂(0), F₃(0), F₄(0), and F₅(0) are

Substituting Eq. (39) into Eq. (45), we obtain

Based on the initiazl terms W₁(0), F₁(0), F₂(0), F₃(0), F₄(0), and F₅(0), and Eq. (62), we could obtain G(k) in Eq. (39) (K(S)≠0) by calculating the parameters W₁(k -1), W₂(k -1), W₃(k - 1), and W₄(k - 1) in Eq. (46) term by term.

5 Numerical simulation

In this section, we apply the method in Section 4 to study two examples.

Example 1 D(S) = 0.138 243S⁸, K(S) = 0, S₀ = 0, S_L = 1^[3].

According to Eq. (25) in the case 1, the intermediate variable ξ is

and the Boltzmann variable φ(= $\frac{x}{{\sqrt t }}$) has the form which is similar to Eq. (26).

By choosing k = 1, 2, 3, respectively, we obtain U(k) in S ∈ (0, 1] as follows:

k = 1, U(0) = 0, U(1) = 0.197 182;

k = 2, U(0) = 0, U(1) = 0.203 251, U(2) = -0.012 703;

k = 3, U(0) = 0, U(1) = 0.202 618, U(2) = -0.012 663 6, U(3) = 0.001 320 735.

Substituting U(k) into Eq. (64), we obtain 1-3 orders approximation solutions as follows:

1 order approximation solution

2 orders approximation solution

3 orders approximation solution

In Eqs. (65)-(67), the intermediate variable ξ is given in Eq. (63).

Compared with the finite element method (FEM) in the space coordinate x varying from 0 to 0.25mm, which is divided into 280 units and calculated for results with saturation S in t = 1min as the exact solution, the present method shows the results of 1-3 orders approximations of Eqs. (65)-(67) in Fig. 1. The results of saturation S obtained from Eq. (67) are compared with the results from the FEM in Table 1.

From Fig. 1, it can be seen that the saturation S decreases from 1 to 0 as the Boltzmann variable φ increases. Due to the effect of hydraulic diffusivity D(S), the saturation S decreases very fast from 0.7 to 0.1 when φ is changed from 0.18 to 0.19mm·min^{-$\frac{1}{2}$} . The solutions of the present method are closer to the exact solutions as the order k increases from 1 to 3 in Fig. 1.

In Table 1, the results of the 3 orders approximation solution are very close to the results of the FEM. The maximum relative error value of present method is -0.227 56% in S = 0.4 for the 3 orders approximation solution. In the meantime, using the DTM without introducing the intermediate variable ξ, we obtain the 1-3 orders approximate solutions by solving Eq. (19) and boundary conditions (23) and (24) with D(S) = 0.138 243S₈, K(S) = 0, S₀ = 0, and S_L = 1. The corresponding results are shown in Fig. 2.

Compared with Fig. 1, the results of 1-3 orders approximate solutions without the interme- diate variable in Fig. 2 are more difficult to converge to the results obtained from the FEM.

Example 2 D(S) = 12 250.497 546 61S⁶, K(S) = 6 078.024 95S¹¹, S_L = 0.36, and S₀ = 0.042 84^[28].

According to the case 1 and Eq. (40), the 1 order approximation solution of K(S) = 0 is

and the saturation S could be expressed as

From the case 2, supposing the 1 order approximation solution of S_r as the intermediate variable, the tiny time point △t = 1 × 10^-6s and choosing k = 1, 2, 3 in Eq. (69), respectively, we obtain G(k) in S ∈ (0.042 84, 0.36] as follows:

Substituting G(k) into Eq. (69), we obtain that 1-3 orders approximation solutions as follows:

1 order approximation solution

2 orders approximation solution

3 orders approximation solution

According to Eq. (69), the intermediate variable S_r in Eqs. (71)-(73) is

Compared with the FEM in x varying from 0 to 4mm, which is divided into 800 units as the exact solution, the results of 1-3 orders approximations of Eqs. (70)-(72) are shown in Fig. 3. The accuracy of saturation S obtained from the 3 orders approximation of Eq. (72) is compared with that from the FEM and presented in Table 2.

According to Fig. 3, the saturation S decreases from 0.36 to 0.042 84 when the infiltration depth x increases. The solutions of the present method are closer to the exact solutions as the order k increases from 1 to 3 gradually. In Table 2, the differences of results between the present method and the FEM are small. The maximum relative error value of the present method is -3.454 29% when x = 3.42mm in the 3 orders approximation solution.

6 Conclusion

In this paper, we have applied the DTM for solving the Richards’ equation with the Brooks- Corey model. By introducing an intermediate variable ξ into the governing differential equation, we successfully increase the convergence rate of solutions with the DTM and derive the 1-3 or- ders approximate solutions of K(S) = 0. For the form of K(S)≠0, we solve the discontinuous solution problem by introducing a tiny time point △t and do the differential transformation based on the corresponding 1-3 orders approximate solutions of K(S) = 0 at the time t = △t. The presented examples demonstrate the accuracy of the present solution by comparing the present results with the results obtained from the FEM.