BIBECHANA Vol. 22, No. 1, April 2025, 30-40 ISSN 2091-0762 (Print), 2382-5340 (Online) Journal homepage: http://nepjol.info/index.php/BIBECHANA Publisher:Dept. of Phys., Mahendra Morang A. M. Campus (Tribhuvan University)Biratnagar Richards equation for the assessment of landslide hazards Ramesh Chandra Timsina Department of Mathematics, Patan Multiple Campus, Tribhuvan University, Kathmandu, Nepal ∗Corresponding author. Email: rameshtimsina028@gmail.com Abstract In this study, we investigate numerical simulation models for water flow in variably sat- urated (unsaturated) soils. These models are crucial for addressing soil-related challenges and analyzing water-related risks, particularly in the context of water resource management, soil water-induced disasters, and the agricultural impacts of global environmental changes. The Richards equation is one of the most widely used models for simulating water flow in porous media, especially in unsaturated soils. However, as a highly nonlinear parabolic par- tial differential equation (PDE), it has limited analytic solutions, which often lack precision in practical scenarios. This necessitates the development of innovative and robust numer- ical methods for accurate simulations. We introduce a numerical procedure that linearizes the Richards equation using the Kirchhoff integral transformation, followed by discretization with various time-stepping schemes. This approach enables efficient and accurate modeling of water flow. To extend its application, we integrate the numerical solution with a hydrological infinite slope stability model to evaluate landslide hazards. Specifically, we calculate the factor of safety index based on an axisymmetric form of the Richards equation, which helps identify potential landslidenprone areas. Furthermore, our model provides a framework for predicting landslides by considering the interplay between water flow and the physical, geological, and topographical characteristics of a landscape. This integrated approach offers valuable insights for geohazard assessment and the mitigation of water-induced risks. Keywords Richards equation; Moisture content; infinite sloe model; Landslide hazards; Factor of safety. . Article information Manuscript received: May 6, 2024; Revised: January 4, 2025; Accepted: January 13, 2025 DOI https://doi.org/10.3126/bibechana.v22i1.65508 This work is licensed under the Creative Commons CC BY-NC License. https://creativecommons. org/licenses/by-nc/4.0/ 1 Introduction Landslides occur in many regions worldwide, partic- ularly in mountainous areas. Geologically, the Hi- malayan mountain chain, stretching approximately 2,400 kilometers, is among the most tectonically ac- tive mountain ranges on Earth [1]. This vast terrain is home to millions of people and diverse ecosys- tems, spanning countries such as Nepal, Bhutan, northern Pakistan, northern India, and other parts of Asia. The Himalayan region faces significant land sliding problems due to a combination of factors: uneven topography, dynamic geological structures, fragile and fractured rock formations, 30 http://nepjol.info/index.php/BIBECHANA rameshtimsina028@gmail.com https://doi.org/10.3126/bibechana.v22i1.65508 https://creativecommons.org/licenses/by-nc/4.0/ https://creativecommons.org/licenses/by-nc/4.0/ Ramesh Chandra Timsina/ BIBECHANA 22 (2025) 30-40 31 and heavy, concentrated rainfall during the mon- soon season. These elements collectively contribute to frequent and severe landslides, posing substan- tial risks to lives, infrastructure, and the envi- ronment in this region. The increasing frequency of severe and concentrated rainfall events, largely driven by climate change, has significantly exacer- bated landslide risks worldwide, including in the Himalayan region. Heavy rainfall, a primary trig- ger for infiltration-induced landslides, is a common occurrence across this mountainous terrain. These landslides lead to extensive destruction of natural infrastructure, damage to roads, houses, and other man-made constructions, and result in substantial annual economic losses and fatalities. In unsatu- rated soils, infiltration-induced landslides typically occur close to the ground surface, where heavy rain- fall triggers shallow landslides. These shallow land- slides commonly develop within the vadose zone, the unsaturated region above the water table, often destabilizing surface layers in foothill areas. Slope failure, which governs the stability and instability of soil surfaces, is closely linked to the infiltration process. Fluctuations in matric potential within unsaturated soils play a crucial role in this pro- cess. As rainfall infiltrates the soil, the resulting changes in matric potential reduce mechanical sta- bility, triggering landslides. Understanding these dynamics is critical for developing strategies to mit- igate landslide hazards in regions prone to heavy rainfall. Infiltration-induced landslides in the Hi- malayas cause severe damage to lives, property, in- frastructure, and the environment, especially dur- ing the monsoon season. To address these chal- lenges, it is essential to examine the factors, under- lying issues, and environmental impacts associated with such landslides. The primary cause of these landslides is surface failure, which is influenced by the interconnectedness of slope instabilities. Ana- lyzing slope stability in unsaturated soils is crucial, as it provides insights into the mechanical behav- ior of slope materials. A detailed study of these materials enables the calculation of the factor of safety (SF ), a conventional measure for assessing the stability of sliding surfaces along failure planes. Certain weather conditions, including continuous and heavy rainfall, rapid snow melt, and abrupt weather changes, are key drivers of slope instability. These phenomena are closely tied to fluctuations in pore water pressure within the soil. During infil- tration and redistribution processes, such as rain- fall and post rainfall periods, pore water pressure exhibits abnormal variations. Increased pore satu- ration leads to a rise in pore water pressure, which, when repeated cyclically, reduces the soil’s effective stress and triggers slope failure. The Richards equa- tion, a nonlinear partial differential equation, effec- tively governs these processes and can be applied for quantitative landslide hazard assessment. By cou- pling the Richards equation with a slope stability model, it is possible to develop a robust framework for analyzing slope failure in relation to weather data. This integrated approach provides a power- ful tool for capturing the dynamics of slope stability under various environmental conditions [2]. The ac- curacy and reliability of this method depend heav- ily on the quality and availability of soil parame- ters, weather data, and the physical characteristics of the landscape. By incorporating these factors, a dependable procedure for landslide hazard assess- ment can be tailored to specific regions, enabling precise and effective mitigation strategies [3]. In this work, we propose an approach to mod- eling infiltration induced landslides, focusing on the critical role of infiltration processes in predict- ing landslide hazards. To achieve this, the Kirch- hoff transformed Richards equation [4] is integrated with the infinite slope stability model [5], provid- ing a robust framework for landslide prediction. This approach incorporates various hydraulic con- ductivity functions, which are experimentally ver- ified and describe the movement of water through pore spaces. These functions depend on soil and fluid properties, such as intrinsic permeability, den- sity, viscosity, degree of saturation, shear stress, and strength. The relationships between hydraulic conductivity and moisture content are captured us- ing established empirical formulations [6], which are critical for understanding soil behavior under dif- ferent conditions. To account for slope stability, a hydrological slope stability model is coupled with the nonlinear Richards equation. This integration enables the development of a comprehensive formu- lation for assessing landslide hazards. The new for- mulation emphasizes infiltration-induced landslides in unsaturated soils and captures the complex in- terplay between soil moisture dynamics and slope stability. The validity of the proposed approach will be tested by comparing its predictions with inter- polated exact data. Applications and examples will also be conducted to illustrate the practical utility of the method, demonstrating its effectiveness in accurately predicting landslide hazards in diverse scenarios. 2 Background This work aims to provide an integrated under- standing of hydrogeology and the mathematical framework necessary for studying water flow in un- saturated soils and infiltration- induced landslide hazards. A key focus is on deriving water flow rela- tions for unsaturated soils from fundamental equa- tions at the pore scale, offering a foundation for developing governing equations for modeling wa- ter flow [4]. The study of water flow in unsatu- Ramesh Chandra Timsina/ BIBECHANA 22 (2025) 30-40 32 rated porous media requires a well-defined formu- lation of governing equations, supported by appro- priate constitutive relationships [6], and accompa- nied by specific boundary and initial conditions to complete the model. Typically, water flow through unsaturated soil is described using the Richards equation, which is based on semi-empirical foun- dations first developed by Buckingham (1907) [7] and later refined by Richards (1931) [4]. Despite its limitations and drawbacks, the Richards equa- tion remains the most widely used model for simu- lating water flow through both saturated and unsat- urated soils. Its long nstanding history of applica- tion highlights its significance, particularly in infil- tration studies across various disciplines. However, solving the Richards equation is challenging due to the nonlinearities in constitutive relationships, es- pecially those involving hydraulic properties. While analytical solutions exist, they are rare and limited to specific cases requiring simplifications or substi- tutions. Consequently, numerical solutions have be- come the dominant approach, with numerous tech- niques developed to address these challenges. How- ever, achieving high accuracy, robustness, and cost effectiveness remains an ongoing effort. This work emphasizes the critical role of solving the Richards equation in understanding complex flow phenom- ena in unsaturated soils. The discussion extends to the numerical challenges and computational ad- vancements required for robust solutions. To ad- dress these issues, a qualitative analysis is presented to better understand the behavior of the Richards equation and its implications for modeling and sim- ulation. Lots of numerical solution techniques were employed efficiently to solve the Richards equation with the use of Newton–Raphson or Picard iter- ative techniques together with finite difference or finite element or finite volume approximations. Fi- nite difference methods were used by Celia and Bouloutas [8] or Woodward and Dawson [9] or Do- gan and Motz [10]. Finite element methods are found in Celia and Bouloutas [8], Lehmann and Ackerer [ [11], Simunek [12], Forsyth et al. [13], Katvetsi et al. [14]. Mixed finite elements were used by Diersch and Perrochet [15], Knabner and Schneid [16], Chavent and Roberts [17]or Fahs et al. [18]. Manzini and Ferraris [19], Orgogozo et al. [20] or Caviedes–Voulli‘eme et al. [21] or C.E. Zambra et al. [22] were used finite volume method. Primary variable switching between the pressure head and water content; upstream weighting of the relative permeability; methods of lines and trans- formation techniques were used by Celia et al. [8]; Kirkland et al. [23]; Tocci et al. [24]; Williams et al. [25]; Mathews et al. [26]; Ji et al. [27]; Li et al. [28]; Zha et al. [29]; Zhang et al. [30]. Gener- ally, numerical methods are adopted to the pressure head–based form of Richards equation in saturated, unsaturated and for layered soils. However, there is a poor mass balance for unsaturated soils with unacceptable time stepping limitations. The meth- ods which are used in pressure head–based form are also applied to mixed form of Richards equa- tion by Celia et al. [8]with a strategy to evaluate the change in water content over one time step di- rectly from the change of the pressure head. It was found that in this strategy, there is a good mass bal- ance. But according to Kirkland et al. [23] under very dry initial conditions, there was a CPU effi- ciency problem in this method. Numerical methods of the moisture–based form of Richards equation has a significantly improved performance when the infiltration phenomenon into very dry soil ( Hills et al. (1989) [31], Kirkland et al. (1992) [23]; Forsyth et al. (1995) [13]; Zha et al. (2013) [32]). But the topmost drawback of this method is that it can be used for simulating water flow in unsatu- rated region only. To capture and retain the benefi- ciary of the moisture based method Wu and Forsyth (2001) [13]; Kirkland et al. [23] and Zeng et al. (2018) [33] were improved by combining these two. Kirkland (1992) [23] was introduced a new variable having different characteristic at unsaturated and saturated region in which it may be a linear func- tion of moisture content for unsaturated soil and a linear function of pressure head for saturated or nearly saturated soil. For relatively dry initial con- ditions this method has a very good CPU efficiency found by Williams et al. (2000) [25]. Moreover, there should be defined a point in which these two methods are joined which may directly affect the accuracy and robustness of the numerical results. Haverkamp et al. (1977) [6]; Campbell (1985) [34]; Vauclin et al. (1979) [35]; Ross and Bristow (1990) [36]; Williams et al. (2000) [25]; Ji et al. (2008) [27]; Berninger et al(2011) [37]; Zhang et al. (2015, 2016, 2018) [30]; Li et al. (2016) [28] used Kirchhoff integral transformation approach to solve the non-linear Richards equation. Pan and Wierenga (1995) [38] laid out a special numerical method to overcome the difficulties of evaluating the integral in Kirchhoff transformed approach. A second integral transform formula was intro- duced by Williams et al. (2000) [25], interpolating the linear combination of moisture content with an arbitrary constant and the Kirchhoff integral transform. Indeed, this method had really a good efficiency and robustness as compared to others. To solve unsaturated flow equation with Kirchhoff transformation Zhang et al. (2015) [30] introduced a new numerical method as a good name Finite Analytic Method (FAM). In this method, a local analytic solution is used to form a set of algebraic equation on the Kirchhoff transformed variable. This method established an accurate and efficient solution. Zhang et al. (2016) [30], applied FAM on Ramesh Chandra Timsina/ BIBECHANA 22 (2025) 30-40 33 mixed form of Richards equation. Also Zhang et al. (2018) [30], extended it into the two–dimensional Richards equation. However, its good achievement, the Finite Analytic Method is limited only for ho- mogeneous porous media soil. Ji et al. (2008) [27], used the transient pressure head–based method and compared a steady state solution obtained from Kirchhoff transformation to a steady–state solution obtained from time marching scheme. They used Gardner’s constitutive relationship [39] and their result showed that the computational cost in unsat- urated flow simulation can greatly be reduced and the use of Kirchhoff transformation has great impli- cation for subsurface water problems. But they only considered pseudo–heterogeneous layered porous media and limited it only to steady–state condi- tions. Although, Kirchhoff transformation method has some limitations, it seems to be very promising method for simulation of Richards equation in un- saturated porous media region. Following into this method, numerical errors will be reduced effectively because variations of Kirchhoff head in its integral nature are much smaller than those in pressure head. Indeed, when the Kirchhoff transformation is employed along with some specific constitutive relations to solve nonlinear Richards equation, a full linearization of the unsaturated Richards equa- tion is permitted. This linearization has an ability of analytic and semi–analytic solutions. Further- more, Numerical Methods for solving nonlinear Richards equation are pertaining to specific consti- tutive relationships. The popular widespread con- stitutive relationships were developed by Gardener, van Genuchten, Brooks and Corey and Haverkamp et al. [6]. Correspondingly, they gave the well measured relationship between hydraulic conduc- tivity, pressure and moisture content. For example, Heejun Suk et al. (2019) [40], obtained a numeri- cal solution of the Kirchhoff transformed Richards equation in variably saturated flow in heteroge- neous layered soil. This method has followed the Gardener’s relationships and is used to solve the linearize partial differential equation. Egidi et al. (2018) [41], used van Genuchten consecutive relationship for a numerical solution of Richards equation providing a simple method adaptable in parallel computing. A constitutive relationship de- veloped experimentally by Haverkamp et al. were used by Liu Fengnan et al. (2020) [42]to a lin- earized finite difference scheme for the Richards equation under variable–flux boundary conditions. Among these, a most robust numerical method an explicit stabilized Runge–Kutta–Legendre Super Time–Steeping scheme is developed [42]to solve mixed form of nonlinear Richards equation with no source and sink terms and with sink terms as evap- otranspiration [43] numerically adopting the con- stitutive relationship developed by Haverkamp et al. [6]. The scheme is based on the natural phenom- ena taking into account the hydraulic conductivity as a function of water pressure head as a real situ- ation with experimental data. Finally, Landslides in the Himalayan region are driven by a variety of interconnected processes and causes. The primary factors contributing to landslide hazards include geological, morphological, physical, and human- induced causes. Among these, the foothills of the Himalayas are particularly vulnerable to landslides triggered by unusual, heavy, and concentrated pre- cipitation during the monsoon season. Rainfall is the most significant triggering factor, with the on- set of the monsoon frequently leading to widespread landslide occurrences and disasters. These events pose severe risks to people, infrastructure, and ecosystems, affecting countless lives and livelihoods on both large and small scales throughout the re- gion. This paper is organized as follows : in section 3, we present the methodology along with Kirch- hoff transformed axi-symmetrical form of Richards equation, including the numerical methods based on finite difference schemes with different time - stepping and the slope failure models for landslide stability in unsaturated soils. In section 4, we show the results of the numerical procedure implemented in python. In section 5, we discuss the results and the future development of this work. 3 Methodology The flow of water in unsaturated soil follows the Darcy law and mass conservation equation [4]. Thus, the flow phenomena can be expressed by Richards equation, a highly non linear PDE. How- ever, the analytical solution of this equation is im- possible, we established a robust numerical solution of this equation in (1+1)D and (2+1)D by adopt- ing different time - stepping schemes and compared their results and used the appropriate one to find the factor of safety for the prediction landslide in the potential landslide area. We have used the (1+1)D solution for bench marking and obtained the solution for axi-symmetrical form of this equa- tion. 3.1 Richards equation The axi-symmetrical form of Richards equation can be written as ∂θ ∂t = 1 r ∂ ∂r ( rK(ψ) ∂ψ ∂r ) + ∂ ∂z ( K(ψ) ∂ψ ∂z ) +W, (1) where θ(ψ) and K(ψ) represent the moisture con- tent θ(ψ) and hydraulic conductivity K(ψ), respec- Ramesh Chandra Timsina/ BIBECHANA 22 (2025) 30-40 34 tively, and W is the rate of precipitation. The function θ(ψ) and K(ψ) are taken from empirical formulas, the Van Genuchten and Haverkamp et al. models [6, 44, 45] are the most widely used in scientific computations due to their formulation by smooth function. Commonly used versions of these models can be written as: K(ψ) = Ks ( 1− (αψ) n−1 )2 / (1 + (αψ)n) m/2 , (2) θ(ψ) = θr + θs − θr (1 + (αψ)n) m , (3) where α is the reciprocal value of ψ0, i.e., the air entry point [46], θs and θr are the saturated water content and residual water content respectively. Ks is the value of the permeability when the soil is sat- urated. m, n are empirical parameters depending on the soil, also m and n have the relation: m = 1− 1 n . (4) The various hydraulic parameters, e.g., α, θr, θs, n, Ks, are related to the pore size distribution and pore geometry, so they ultimately depend on the soil type and can be determined experimentally by specific laboratory tests. For implementation and to improve the results, we used the Van Genuchten model. 3.2 Linearization Techniques Applying the Kirchhoff integral transformation in equation 1, we let h = ψ − z and define: ϕ(h) = ∫ h 0 K(h′)dh′. (5) Since K(h) > 0 from equation (2), the function ϕ(h) is strictly increasing with K(h) = K(ψ). Tak- ing derivatives of both sides of the transformation with respect to r and z, we obtain: ∂ϕ ∂r = K(ψ) ∂ψ ∂r , ∂ϕ ∂z = K(ψ) ∂ψ ∂z . (6) Again, taking the derivative of equation (6) with respect to r: ∂2ϕ ∂r2 = K(ψ) ∂2ψ ∂r2 + ∂K ∂ψ ( ∂ψ ∂r )2 , (7) and ∂φ ∂z = ∂φ ∂h ∂h ∂z = K(h) ∂(ψ − z) ∂z = K(ψ) ( ∂ψ ∂z − 1 ) = K(ψ) ∂ψ ∂z −K(ψ) (8) Again differentiating of equation (8), ∂2φ ∂z2 = ∂ ∂z ( K(ψ) ∂ψ ∂z ) − ∂ ∂z (K(ψ)) (9) Using the equations (6),(8) and (9), the Richards equation (1), takes the form ∂θ ∂t = ∂2φ ∂r2 + 1 r ( ∂φ ∂r ) + ∂2φ ∂z2 +W (10) with θ(φ) = θ(h). The corresponding initial and boundary conditions for the transformed equation (10) takes the following form .θ(r, z, 0) = φ0(r, z), Rin ≤ r ≤ Rout, 0 ≤ z ≤ Ztop ∂φ ∂z = q(t), r > 0, z = 0, t > 0 φ(r, Zbot, t) = β(t), r > 0, t > 0 ∂φ ∂r = 0, r = Rout , z > 0, t > 0 .φ(Rin,Z, t) = β1(t), t > 0 (11) The Kirchhoff transformation transformed the nonlinear equation (1) to a nonlinear parabolic equation (10). Also we note that the Kirchhoff transformation preserves the uniqueness of the so- lution for the transformed problem. 3.3 Discretization techniques and imple- mentation We have the transformed equation is in two differ- ent state variables. To solve the transformed equa- tion (10) numerically with the prescribed initial and boundary conditions (11), it is feasible to have a sin- gle state variable. For this, and are assumed as single valued continuous functions of one another, and arranging these variables as ∂θ ∂t = ∂θ ∂φ ∂φ ∂t = ( 1 ∂φ ∂θ ) ∂φ ∂t , ∂φ ∂θ = ∂φ ∂h ∂h ∂θ , (12) Differentiating (3) with respect to h, we get ∂θ ∂h = αn(θs−θr)(−m)(1+ |αh|n)−m−1n|h|n−1 ∂φ ∂h = K(h) = K(ψ), (13) Using (12) and (13), the transformed Richards equation (10) takes the form c(φ) ∂φ ∂t = ∂2φ ∂r2 + 1 r ( ∂φ ∂r ) + ∂2φ ∂z2 (14) Ramesh Chandra Timsina/ BIBECHANA 22 (2025) 30-40 35 where the functional coefficient c depends on through h as c(φ(h)) = −α nmn(θs − θr)|h|n−1 K(h)(1 + |αh|n)m+1 (15) 3.4 Finite Difference Scheme We set up a two dimensional (r, z ) uniform grid for an axi symmetric problem in the cylinder geom- etry by subdividing the radial length [Rn, Rout ] into Mr subintervals of width ∆r = Rout−Rin Mr and the height [0, Ztop] in to Mz subintervals of width ∆z = Ztop Mz . We construct a grid (ri, zj , tn) with ri = i∆r, i = 0, 1, 2, ....,Mr, zj = j∆z, j = 0, 1, 2, 3, ...Mz, and tn = n∆t, n = 1, 2, 3..., N . Let φni,j denote φ(ri, zj , tn). The partial differential equation (14) can be approximated using forward difference in time and central difference in space as ∂φ ∂t |(rizj ,tn) ≈ φn+1i,j − φi,jn ∆t , ∂φ ∂r |(rizj ,tn) ≈ φni+1,j − φi−1,jn 2∆r ∂2φ ∂r2 |(rizj ,tn) ≈ φni−1,j − 2φi,jn + φni+1,j ∆r2 , ∂2φ ∂z2 |(rizj ,tn) ≈ φni,j−1 − 2φi,jn + φni,j+1 ∆z2 (16) Using a weighted average of the derivatives( ∂φ ∂r , ∂2φ ∂r2 , ∂2φ ∂z2 ) at the two time levels, tn and tn+1 and adopting Crank-Nicolson(CN) scheme , equa- tion (14) can be discretized as φn+1 i,j − φn i,j ∆t = 1 2cni,j(∆r) 2 [ φn+1 i−1,j − 2φn+1 i,j + φn+1 i+1,j +φn i−1,j−2φn i,j+φ n i+1,j+ 1 2cni,j(∆z) 2 [ φn+1 i,j−1 − 2φn+1 i,j +φn+1 i,j+1+φ n i,j−1−2φn i,j+φ n i,j+1+ 1 4ricni,j(∆r) [ φn+1 i+1,j − φn+1 i−1,j + φn i+1,j + φn i−1,j ] (17) We collect the unknowns on the left hand side : − ( Fr − F ri ) φn+1 i−1,j+(1+2Fr+2Fz)φ n+1 i,j − ( Fr + F ri ) φn+1 i+1,j − Fz(φ n+1 i,j−1 + φn+1 i,j+1) = ( Fr − F ri ) φn i−1,j+(1−2Fr−2Fz)φ n i,j+ ( Fr + F ri ) φn i+1,j − Fz(φ n i,j−1 + φn i,j+1) (18) where Fr = ∆t 2cni,j(∆r)2 , Fz = ∆t 2cni,j(∆z)2 , F = ∆t 4cni,j(∆r) . The equation (18) is coupled at the new time level n + 1. That is, we must solve a system of (linear) algebraic equations, which we will write as Ax = B, where A is the coefficient matrix, x is the vector of unknowns, B is the right hand-side. To solve the system of linear equations, we need to form a matrix system Ax = B, where the solution vector x must have one index. For this, we need a numbering of the unknowns with one in- dex, not two as used in the mesh. We introduce a mapping position (i, j ) = v(i, j ) from a mesh point with indices (i, j ) to the corresponding unknown p in the equation system. p = v(i, j ) = j (Vr + 1) + i, for i = 0, 1, 2, · · · , Vr, j = 0, 1, 2, · · · , Vz, With this mapping, we number the points along the radial direction start- ing with z = 0 and then filled one mesh line at a time. In another way p = v(i, j ) = i(Vz + 1) + j, for i = 0, 1, 2, · · · , Vr, j = 0, 1, 2, · · · , Vz. with r = 0 and then filled one mesh line at a time. From this we can get the general feature of the coefficient matrix obtained from the discretized equation (18). Let Ak,l be the value of element (k, l) in the coefficient matrix A, where k and l are the numbering of the unknowns in the equation system. We have Ak,l = 1 for k = l corresponding to the all known boundary values. Let k be v(i, j ), i.e., the single index corresponding to the mesh point (i, j ). Then, for inte- rior mesh along with boundary, we have Av(i,j),v(i,j) = Ak,k = 1 + (Fr + Fz), Ak,v(i−1,j) = Ak,k−1 = −Fr Ak,v(i+1,j) = Ak,k+1 = −Fr Ak,v(i,j−1) = Ak,k−(vz+1) = −Fz Ak,v(i,j+1) = Ak,k+(vz+1) = −Fz The corresponding right hand side vec- tor in the equation system has the en- tries Bk, where k numbers the equations with the given boundary values. The above mention algorithm can be used to update the transformed variable φn i,j to its value in the next time level φn+1 i,j . But we cannot advance the algorithm to the next time level φn+2 i,j without evaluating the function c(φn+1 i,j ) Ramesh Chandra Timsina/ BIBECHANA 22 (2025) 30-40 36 which requires computing the intermediate variable hn+1 i,j . For this , we imply the equation (13) which can be approximated as hn+1 i,j = hni,j + φn+1 i,j − φn i,j K(hni,j) . (19) 3.5 Landslide Stability analysis in unsatu- rated soil This study focuses on the modeling of infiltration- induced landslides in unsaturated soils, specifi- cally analyzing slope stability and instability using the axi-symmetric Kirchhoff-transformed Richards equation. The approach considers the variation of moisture content with water pressure head and vice versa. Infiltration- induced landslides on ground surfaces often resemble shallow landslides triggered by intense rainfall events. The problem is addressed using an infinite slope model to calculate the safety factor (SF), which serves as a key indicator of slope stability. The safety factor is defined as the ratio of resisting forces that prevent slope failure to driving forces that cause collapse. As (SF) greater than one indicates stability, while as (SF) less than one sig- nals instability. This numerical value acts as a land- slide hazard index, enabling predictions about the likelihood of landslides for given topographic con- ditions. The infinite slope model is adopted due to its simplicity and effectiveness, especially for shal- low landslides where the failure plane is parallel to the ground surface and relatively shallow compared to the slope length. For inter connected to the moisture profile obtained from the above solution schemes we preferred the following hydrological in- finite slope model SF = C zYtsinβcosβ + tanφ tanβ + hθY wtanφ zYtsinβcosβ (20) where C is the effective cohesion, is the moistu- rity, Yt is the unit weight of the soil, Yw is the unit weight of the water, is the angle of internal friction, z is the soil thickness, is the slope of the inclined surface. The combination of the infinite slope model (Equation 20) with Richards equation, allows the landslide hazard evaluation directly from weather data. In particular, the solution of Richards equa- tion will be used to obtain information about in the expression. 4 Simulation Results 4.1 Numerical Setup To demonstrate the numerical solution of the pre- scribed model, the numerical process developed in the above section is written in python and ran on a laptop with 2.8 GHz Quad-Core Intel Core i7 processor. We considered a specific infiltration ex- periment in unsaturated soil and observe the cor- responding outcomes. At first the setup consisted on a landslide area having dimension 70m × 100m. The physical constituent of the landslide area is sandy soil. We use the soil parameters and charac- teristics relationship between the soil moisture con- tent () and the hydraulic conductivity K () from the work of Van Genuchten. K(ψ) = Ks [1− |∂ψ|n]− m 2 [ 1− ( ||n 1 + ||n )m]2 (21) θ(ψ) = θr + (θs − θr) (1 + ||n)m (22) The simulation starts with a uniform satura- tion θ = 0.1 cm3/cm3 and a constant water head ψ = −61.5 cm is maintained at the bottom bound- ary z = Zbot. For the upper boundary z = Ztop at the soil surface, a constant flux q(t) = 13.69 cm/hr for t < 0.7hr and zero normal flux condition for t > 0.7hr is maintained. To compute the approximate solution using different time-stepping schemes, we used a uniform spatial step size ∆z = 2 cm in the axial direction and a step size of ∆z = 2 cm in the radial direction. The simulation was run for 0.75 hr. For the infiltration-induced landslide, we inter- polated the moisture variation obtained from the above experiment to the hydrological model for sta- bility analysis. Data was obtained from the Depart- ment of Hydrology and Meteorology, Government of Nepal, and field observations in the landslide area located in Yangbarak VDC, Ward No. 6, Shripung, Panchthar district. Rainfall data from the region was applied in our simulation model (Equation 20). 5 Results and Discussion The axi symmetrical model was employed to sim- ulate vertical wa- ter infiltration into an unsatu- rated homogeneous porous medium and to ana- lyze the landslide hazards, hydrological model was used. Numerical experiments were carried out by using sand column (proposed by Van Genuchten.) Ramesh Chandra Timsina/ BIBECHANA 22 (2025) 30-40 37 with finite difference schemes(FTCS, BTCS, CN and RKL). Since we do not have exact solution, we have used numerical solution obtain from above (RKL) finite difference scheme for one dimensional as the reference solution Timsina, R.C et al. [42] for numerical experiment of two dimensional wa- ter infiltration. Fig. (2) (a) shows the variation of volumetric moisture content in sand for one di- mension, and Fig. 2 (b) shows the comparison of different time-steeping schemes as discussed in [42]. In this numerical experiment, the function used for the hydraulic conductivity and the water content as in one dimensional case were taken from Van Genuchten [46]. We set the axi-symmetrical form of Richards equation in an unsaturated flow into a region of sandy soil. For this, we suppose the domain is spec- ified to be an annulus, rinner ≤ r ≤ router with rinner strictly greater than zero, other than a disk. The annular domain consists of 70 cm in length with 100 cm of annular radius. We suppose it vertically downward where the z-axis is taken as downward positive with a constant water head ψ = −61.5 cm at the left (Dirichlet) boundary condition r = Rin as dummy. A constant flux boundary condition (Neumann) on the right r = Rout is taken and a constant water head ψ = −61.5 cm at the bottom boundary z = Zbottom. At the upper boundary z = 0 (the soil surface), a constant flux q(t) = 13.69 cm/hr for t < 0.7 hr and zero normal flux condition for t > 0.7 hr. The solution domain was meshed us- ing a spatial step size of ∆z = 2 cm on the axial di- rection and a step size of ∆r = 2.5 cm on the radial direction and the simulation was run for 1 hr. The longitudinal water content profiles and radial water content profiles at time 0 hr and some other times are shown in Figure (3) (top). From this figure, it is possible to note that the residual water content at the bottom of the annulus is reached when time t is about 0.50 hr. The moisture content contours are shown in Figure (4) and corresponding surface plots of moisture content are in Figure (4) (left). We observe the flow phenomena and carried out the variation of water content in the region we have prescribed. We have interconnected the above re- sult to the infinite slope model to evaluate the safety factor. We have measured the slope angle and the thickness of the unsaturated soil as 300 and 8m. We use the parameter of effective cohesiveness is 4.4kPa, the frictional angle is 280, the soil unit weight is 23.4kN/m3. Figure 6 depicts the com- puted results of safety factor for different pressure head observed in the experiment. The factor of safety obtained by the above model demonstrates the possibilities of surface failure and operation of landslide hazards. From figure 6 we conclude that the situation of surface failure may increases when the level of ground water is increased by heavy rain- fall meanwhile the pressure head energy increase and the level of unsaturated region decrease and the possibilities of landslide increases. The axisymmetric form of Richards’ equation was set in an unsaturated flow region of sandy soil. The domain was specified as an annulus, rinner ≤ r ≤ router, with rinner > 0. The annu- lar domain is 70 cm in length and 100 cm in radius. A constant water head ψ = −61.5 cm was main- tained at the bottom boundary, while a constant flux q(t) = 13.69 cm/hr for t < 0.7 hr and zero nor- mal flux condition for t > 0.7hr were applied at the soil surface. The solution domain was meshed with a spa- tial step size ∆z = 2 cm in the axial direction and ∆r = 2.5 cm in the radial direction. The simula- tion was run for 1 hr. Figures 4 and 5 show the longitudinal and radial water content profiles and the corresponding surface plots. Ramesh Chandra Timsina/ BIBECHANA 22 (2025) 30-40 38 The safety factor for different pressure heads and corresponding moisture content is shown in Figure 6. The slope angle and thickness of the un- saturated soil were 30◦ and 8m, respectively. The effective cohesion was 4.4 kPa, friction angle 28◦, and soil unit weight 23.4 kN/m3. The results indi- cate an increased likelihood of surface failure with heavy rainfall. Conclusion The flow of water in unsaturated soil is described by axi symmetrical Kirchhoff trans- formed Richards equation (a highly nonlinear degenerate parabolic partial differential equation) and solved it numeri- cally and interpolate the result in the hydrological model for stability analyses for the possibilities of landslide hazards. Also this work is based on il- lustrative numerical examples in cylindrical form, with realistic parameters. We implemented differ- ent finite difference schemes. The work presented here describes and verifies the application and accu- racy of finite difference schemes to simulate flow in unsaturated porous media in axi symmetrical cylin- drical formation of the land surface. However, this work presented a new super time stepping numer- ical scheme which was able to solve the Kirchhoff transformed axi symmetrical Richards equation and the findings are converges to the theoretical anal- ysis and also the findings are in line with one di- mensional approaches [42, 43]. Figures 2 and 3 de- picted that the numerical simulations results for the same infiltration experiment as in one dimension imposing, zero flux boundary conditions in the lat- eral boundary mimic the results from the one di- mensional model. That is the longitudinal profile are in line with that of the one dimensional model. The numerical method is able to handle short du- ration infiltration and relatively easy to implement. The interpolated result shows that the prediction of landslide is possible for the available parameter de- pending on the constituent of the soil and nature of the land. This work can be extended to achieve to accurate prediction of massive landslide to unsat- urated heterogeneous soils with abruptly changing wetness conditions. Acknowledgments The author thanks the Research Division, Univer- sity Grants Commission, Sano Thimi, Bhaktapur, for the small research grant provided for this work (SRDIG-80/81-S&T-10). Data Availability The data supporting the findings of this study are included within the article. Conflicts of Interest The author declares no conflict of interest. References [1] R. K. Dahal and S. Hasegawa. Representative rainfall thresholds in the nepal himalaya. Sci- ence Direct, 100:429–443, 2008. [2] R. C. Timsina. Infiltration-induced landslide: An application of richards equation. Journal of Nepal Mathematical Society, 7:86–99, 2024. [3] C. Yu Ku, C. Yu Liu, and T.C. Tsai. A novel radial basis function approach for infiltration- Ramesh Chandra Timsina/ BIBECHANA 22 (2025) 30-40 39 induced landslides in unsaturated soils. Water, 14:1–20, 2022. [4] L. A. Richards. Capillary conduction of liquids through porous mediums. Physics, 1(5):318– 333, 1931. [5] J. M. Duncan, S. G. Wright, and T.L. Bran- don. Soil strength and slope stability. John Willey and Sons, Hoboken, NJ, 2014. [6] R. Haverkamp, M. Vauclin, J. Touma, J. P. Wierenga, and G. Vachaud. A compari- son of numerical simulation models for one- dimensional infiltration. Soil Sci. Soc. Am. J., 41(2):285–294, 1977. [7] E. Buckingham. Studies on the movement of soil moisture, volume 38. Bulletin Edition, Washington, Govt. Print. Off., 1907. [8] M. A. Celia, E. T. Bouloutas, and R. L. Zarba. A general mass-conservative numerical solu- tion for the unsaturated flow equation. Water Resources Research, 26(7):1483–1496, 1990. [9] C. S. Woodward and C. N. Dawson. Analysis of expanded mixed finite element methods for a non-linear parabolic equation modeling flow into variably saturated porous media. SIAM Journal on Numerical Analysis, 37(3):701–724, 2000. [10] A. Dogan and L. H. Motz. Saturated- unsaturated 3d groundwater model. ii: Veri- fication and application. Journal of Hydrology Engineering, 10(6):505–515, 2005. [11] F. Lehmann and P. Ackerer. Comparison of iterative methods for improved solutions of the fluid flow equation in partially saturated porous media. Transport in Porous Media, 31(3):275–292, 1998. [12] J. Simunek, M. van Genuchten, and M. Se- jna. Development and applications of the hy- drus and stanmod software packages and re- lated codes. Vadose Zone Journal, 7(2):587– 600, 2008. [13] P. A. Forsyth, Y. S. Wu, and K. Pruess. Robust numerical methods for saturated–unsaturated flow with dry initial conditions in heteroge- neous media. Advances in Water Resources, 18(1):25–38, 1985. [14] D. Kavetski, P. Binning, and S. W. Sloan. Adaptive backward euler time stepping with truncation error control for numerical mod- elling of unsaturated fluid flow. International Journal for Numerical Methods in Engineer- ing, 53(6):1301–1322, 2001. [15] H.-J. Diersch, G. Perrochet, and P. et al. On the primary variable switching technique for simulating unsaturated–saturated flows. Ad- vances in Water Resources, 23(3):271–301, 1999. [16] P. Knabner and E. Schneid. Adaptive hybrid mixed finite element discretization of in sta- tionary variably saturated flow in porous me- dia. In Lecture Notes in Computational Science and Engineering, pages 37–44. Springer Berlin Heidelberg, 2002. [17] G. Chavent and J. E. Roberts. A unified phys- ical presentation of mixed, mixed–hybrid finite elements and standard finite difference approx- imations for the determination of velocities in water flow problems. Advances in Water Re- sources, 14(6):329–348, 1991. [18] M. Fahs, A. Younes, and F. Lehmann. An easy and efficient combination of the mixed finite el- ement method and the method of lines for the resolution of richards equation. Environmental Modeling & Software, 24(9):1122–1126, 2009. [19] G. Manzini and S. Ferraris. Mass–conservative finite volume methods on 2-d unstructured grids for the richards equation. Advances in Water Resources, 27(12):1199–1215, 2004. [20] L. Orgogozo, N. Renon, C. Soulaine, F. Hénon, S. K. Tomer, D. Labat, O. S. Pokrovsky, and M. Sekhar et al. An open source massively parallel solver for richards equation: Mecha- nistic modeling of water fluxes at the water- shed scale. Computer Physics Communica- tions, 185(12):3358–3371, 2014. [21] D. Caviedes-Voullième, J. Murillo, P. Garcia- Navarro, and et al. Numerical simulation of ground water surface interactions by external coupling of the 3d richards equation and the full 2d shallow water equations. In XIX Inter- national Conference on Water Resources, 2012. [22] C. E. Zambra, M. Dumbser, E. F. Toro, and N. O. Moraga. A novel numerical method of high-order accuracy for flow in unsaturated porous media. International Journal for Nu- merical Methods in Engineering, 89(2):227– 240, 2012. [23] M. R. Kirkland, R. G. Hills, and P. J. Wierenga. Algorithms for solving richards equation for variably saturated soils. Water Resour. Res., 28(8):2049–2058, 1992. [24] M. D. Tocci, C. T. Kelley, and C. T. Miller. Ac- curate and economical solution of the pressure- head form of richards equation by the method of lines. Advances in Water Resources, 20(1):1– 14, 1997. Ramesh Chandra Timsina/ BIBECHANA 22 (2025) 30-40 40 [25] G. A. Williams, C. T. Miller, and C. T. Kelly. Transformation approaches for simulating flow in variably saturated porous media. Water Re- sources Research, 36:923–934, 2000. [26] C. J. Matthews, R. D. Braddock, and G. Sander. Modeling flow through a one- dimensional multi-layered soil profile using the method of lines. Environmental Modeling and Assessment, 9:103–113, 2004. [27] S. H. Ji, Y. J. Park, E. A. Sudicky, and J. F. Sykes. A generalized transformation approach for simulating steady-state variably-saturated subsurface flow. Advances in Water Resources, 31:313–323, 2008. [28] N. Li, X. Yue, and L. Ren. Numerical ho- mogenization of the richards equation for un- saturated water flow through heterogeneous soils. Water Resources Research, 52:8500– 8525, 2016. [29] Y. Zha, J. Yang, L. Yin, Y. Zhang, W. Zeng, et al. A modified picard iteration scheme for overcoming numerical difficulties of simulating infiltration into dry soil. Journal of Hydrology, 551:56–69, 2017. [30] Z. Zhang, W. Wang, C. Gong, T. J. Yeh, Y. L. Wang, et al. Finite analytic method for mod- eling variably saturated flows. Science of the Total Environment, 621:1151–1162, 2018. [31] D. Hillel. Fundamentals of Soil Physics. Aca- demic San Diego, California, 1980. [32] Y. Zha, J. Yang, L. Shi, and X. Song. Sim- ulating one-dimensional unsaturated flow in heterogeneous soils with water content-based richards equation. Vadose Zone Journal, 12(2):1–12, 2013. [33] J. Zeng, Y. Zha, and J. Yang. Switching the richard equation for modeling soil water move- ment under unfavorable conditions. Journal of Hydrology, 563:942–949, 2018. [34] G. S. Campbell. Soil physics with BASIC. El- sevier, New York, 1985. [35] M. Vauclin, D. Khanji, and G. Vachaud. Ex- perimental and numerical study of a transient, two-dimensional unsaturated–saturated water table recharge problem. Water Resources Re- search, 15:1089–1101, 1979. [36] P. J. Ross and K. L. Bristow. Simulating wa- ter movement in layered and gradational soils using the kirchhoff transform. Soil Science So- ciety of America Journal, 54:1519–1524, 1990. [37] H. Berninger, R. Kornhuber, O. Sander, et al. Fast and robust numerical solution of the richards equation in soil, homogeneous. SIAM Journal on Numerical Analysis, 49(6):2576– 2597, 2011. [38] L. Pan and P. J. Wierenga. A transformed pressure head-based approach to solve richards equation for variably saturated soils. Water Resources Research, 31(4):925–931, 1995. [39] W. R. Gardner. Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Science, 85(4):228–232, 1958. [40] H. Suk and E. Park. Numerical solution of the kirchhoff-transformed richards equation for simulating variably saturated flow in heteroge- neous layered porous media. Journal of Hy- drology, 579:1–13, 2019. [41] N. Egidi, E. Gioia, P. Maponi, and L. Spadoni. A numerical solution of richards equation: a simple method adaptable in parallel comput- ing. International Journal of Computer Math- ematics, 97:1–18, 2018. [42] R. C. Timsina, H. Khanal, and K. N. Uprety. An explicit stabilized runge-kutta-legendre su- per time-stepping scheme for the richards equation. Mathematical Problems in Engineer- ing, 2021, 2021. [43] R. C. Timsina, H. Khanal, A. Ludu, and K. N. Uprety. A numerical solution of water flow in unsaturated soil with evapotranspira- tion. The Nepali Mathematical Sciences Re- port, 38(2):35–45, 2021. [44] L. F. Richardson. Weather Prediction by Nu- merical Process. Cambridge University Press, 1922. [45] R. J. Brooks and A. T. Corey. Hydraulic prop- erties of porous media. Hydrology papers, Col- orado State University, Fort Collins, Colorado, 1964. [46] M. T. van Genuchten. A closed–form equa- tion for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44(5):892–898, 1980. Introduction Background Methodology Richards equation Linearization Techniques Discretization techniques and implementation Finite Difference Scheme Landslide Stability analysis in unsaturated soil Simulation Results Numerical Setup Results and Discussion