Approximate analytical solution for Richards’ equation with finite constant water head Dirichlet boundary conditions

Abstract

The Richards’ equation is widely used to simulate the saturation distributions in porous medium. Due to the high nonlinear of hydraulic diffusivity, it is usually difficult to obtain the analytical solutions, especially for finite boundaries. In this paper, the Richards equation for horizontal infiltration problem with finite Dirichlet boundaries is solved. The whole infiltrating process is divided into three state: the free infiltrating state, the transition state and the steady one, in which the transition state is considered as the key problem for analytical solving. Based on Boltzmann transformation and series expansion technique, an intermediate variable is introduced and an approximate analytical solution is derived for transition state. In addition, the exact solutions for other states are also given in the appendix. The present solutions can be applied for arbitrary nonlinear hydraulic diffusivity in Richards’ equation. Two examples for power law and exponential law diffusivities are shown to confirm the accuracy of present method.

Appendices

Appendix 1

According to Lockington et al. (1999), the free infiltrations in Fig. 2 can be considered as two semi-infinite boundary problems, which mean RE (1) with the definite conditions can be written as:

$$ \left\{ {\begin{array}{*{20}l} {\frac{\partial \theta }{{\partial t}} = \frac{\partial }{\partial x}\left( {K\left( \theta \right)\frac{\partial \psi \left( \theta \right)}{{\partial x}}} \right)} \hfill \\ {\theta {\mid}_{x = 0} = \theta_{L} } \hfill \\ {\theta {\mid}_{x \to + \infty } = \theta_{0} } \hfill \\ {\theta {\mid}_{t = 0} = \theta_{0} ,} \hfill \\ \end{array} } \right. $$

(46)

and:

$$ \left\{ {\begin{array}{*{20}l} {\frac{\partial \theta }{{\partial t}} = \frac{\partial }{\partial x}\left( {K\left( \theta \right)\frac{\partial \psi \left( \theta \right)}{{\partial x}}} \right)} \hfill \\ {\theta \theta {\mid}_{x = L} = \theta_{R} } \hfill \\ {\theta {\mid}_{x \to - \infty } = \theta_{0} } \hfill \\ {\theta {\mid}_{t = 0} = \theta_{0} .} \hfill \\ \end{array} } \right. $$

(47)

According to Parlange et al. (1992), Eq. (46) can be described by the Bruce and Klute equation as:

$$ \left\{ {\begin{array}{*{20}l} {K\left( {\theta + \theta_{0} } \right)\psi^{\prime}\left( {\theta + \theta_{0} } \right) = - \frac{1}{2}\frac{{{\text{d}}\phi }}{{{\text{d}}\theta }}\mathop \int \limits_{0}^{{\theta - \theta_{0} }} \phi {\text{d}}\theta } \hfill \\ {\phi {\mid}_{{\theta = \theta_{L} - \theta_{0} }} = 0} \hfill \\ {\phi {\mid}_{\theta = 0} \to + \infty ,} \hfill \\ \end{array} } \right. $$

(48)

where ϕ (= x/t^½) is Boltzmann variable, ψ' is:

$$ \psi^{\prime} = \frac{\partial \psi }{{\partial \theta }}, $$

(49)

and Eq. (47) can be changed into:

$$ \left\{ {\begin{array}{*{20}l} {K\left( {\theta + \theta_{0} } \right)\psi^{\prime}\left( {\theta + \theta_{0} } \right) = - \frac{1}{2}\frac{{{\text{d}}\phi_{1} }}{{{\text{d}}\theta }}\mathop \int \limits_{0}^{{\theta - \theta_{0} }} \phi_{1} {\text{d}}\theta } \hfill \\ {\phi_{1} {\mid}_{{\theta = \theta_{R} - \theta_{0} }} = 0} \hfill \\ {\phi_{1} {\mid}_{\theta = 0} \to + \infty ,} \hfill \\ \end{array} } \right. $$

(50)

in which ϕ₁ is:

$$ \phi_{1} = \frac{L - x}{{\sqrt t }}. $$

(51)

Chen and Dai (2017) gave the approximate analytical solution of Eq. (45) as:

$$ \phi \approx \mathop \sum \limits_{i = 1}^{n} u_{i} \left( {\mathop \int \limits_{0}^{{\theta - \theta_{0} }} \frac{{K\left( {\theta + \theta_{0} } \right)\psi^{\prime}\left( {\theta + \theta_{0} } \right)}}{\theta }{\text{d}}\theta } \right)^{i} , $$

(52)

where i = 1,2,3,…,n, u_i is the calculation parameter. Substituting Eq. (51) into Eq. (48) and applying the item-by-item derivation to Eq. (48) in $\theta$ = $\theta$ _L-$\theta$ ₀, we have:

$$ \frac{{{\text{d}}^{k} }}{{{\text{d}}\theta^{k} }}\left. {\left[ {K\left( {\theta + \theta_{0} } \right)\psi^{\prime}\left( {\theta + \theta_{0} } \right)} \right]} \right|_{{\theta = \theta_{L} - \theta_{0} }} = - \frac{{{\text{d}}^{k} }}{{{\text{d}}\theta^{k} }}\left. {\left( {\frac{1}{2}\frac{{{\text{d}}\phi }}{{{\text{d}}\theta }}\mathop \int \limits_{0}^{{\theta - \theta_{0} }} \phi {\text{d}}\theta } \right)} \right|_{{\theta = \theta_{L} - \theta_{0} }} , $$

(53)

in which k = 0,1,2,3,…n − 1. Then u_i can be obtained by solved the nonlinear algebraic (52).

Similarly, we obtain ϕ₁ in Eq. (50) as:

$$ \phi_{1} \approx \mathop \sum \limits_{i = 1}^{n} v_{i} \left( {\mathop \int \limits_{0}^{{\theta - \theta_{0} }} \frac{{K\left( {\theta + \theta_{0} } \right)\psi^{\prime}\left( {\theta + \theta_{0} } \right)}}{\theta }{\text{d}}\theta } \right)^{i} , $$

(54)

where v_i is the calculation parameter which can be solved by:

$$ \frac{{{\text{d}}^{k} }}{{{\text{d}}\theta^{k} }}\left. {\left[ {K\left( {\theta + \theta_{0} } \right)\psi^{\prime}\left( {\theta + \theta_{0} } \right)} \right]} \right|_{{\theta = \theta_{R} - \theta_{0} }} = - \frac{{{\text{d}}^{k} }}{{{\text{d}}\theta^{k} }}\left. {\left( {\frac{1}{2}\frac{{{\text{d}}\phi }}{{{\text{d}}\theta }}\mathop \int \limits_{0}^{{\theta - \theta_{0} }} \phi {\text{d}}\theta } \right)} \right|_{{\theta = \theta_{R} - \theta_{0} }} . $$

(55)

Appendix 2

The dash line in Fig. 2 can be described by RE without time t:

$$ \left\{ {\begin{array}{*{20}l} {\frac{\partial }{\partial x}\left( {K\left( \theta \right)\frac{\partial \psi \left( \theta \right)}{{\partial x}}} \right) = 0} \hfill \\ {\theta {\mid}_{x = 0} = \theta_{L} } \hfill \\ {\theta {\mid}_{x = L} = \theta_{R} .} \hfill \\ \end{array} } \right. $$

(56)

Integrated the differential equation in Eq. (56), we have:

$$ K\left( \theta \right)\frac{\partial \psi \left( \theta \right)}{{\partial x}} = C_{1} . $$

(57)

Substituted Eq. (49) into Eq. (57), we obtain:

$$ K\left( \theta \right)\psi^{\prime}\frac{\partial \theta }{{\partial x}} = C_{1} . $$

(58)

Integrating Eq. (58), it is derived:

$$ \mathop \int \limits_{\theta }^{{\theta_{L} }} K\left( \theta \right)\psi^{\prime}\left( \theta \right){\text{d}}\theta = C_{{1}} x + C_{{2}} , $$

(59)

in which C₁ and C₂ is a calculation parameter. Considering the boundary conditions in Eq. (56), the analytical solution for steady state is:

$$ \mathop \int \limits_{\theta }^{{\theta_{L} }} K\left( \theta \right)\psi^{\prime}\left( \theta \right){\text{d}}\theta = \frac{x}{L}\mathop \int \limits_{{\theta_{R} }}^{{\theta_{L} }} K\left( \theta \right)\psi^{\prime}\left( \theta \right){\text{d}}\theta . $$

(60)