How to Use COMSOL Multiphysics for coupled dual-permeability hydrological and slope stability modeling

(1)

Procedia Earth and Planetary Science 9 ( 2014 ) 83 – 90

Selection and peer-review under responsibility of Dipartimento di Ingegneria Civile, Design, Edilizia e Ambiente, Seconda Università di Napoli. doi: 10.1016/j.proeps.2014.06.018

ScienceDirect

The Third Italian Workshop on Landslides

How to Use COMSOL Multiphysics for coupled dual-permeability

hydrological and slope stability modeling

Wei Shao

a,

*, Thom Bogaard

a

, Mark Bakker

a

a_{Water Resources Section, Faculty of Civil Engineering and Geosciences, Delft University of Technology, 2628CN, Delft, Netherlands}

Abstract

Preferential flow paths, such as cracks, macropores, fissures, pipes etc. are common features of highly heterogeneous slopes. During intense rainstorms, preferential flow has a significant influence on subsurface flow and slope stability. Dual-permeability models are widely used to simulate preferential flow, but have not been incorporated into hydro-mechanical models. In this study, the COMSOL Multiphysics software was used to couple a dual-permeability model with a soil mechanics model for landslide stability evaluation on a hillslope scale with examples of up to 100m spatial scale. Detailed information is provided on how to incorporate current hydrological and soil mechanics theories into COMSOL. The model is benchmarked against two existing solutions and is applied to evaluate the effect of preferential flow on slope stability as an example.

Selection and peer-review under responsibility of Dipartimento di Ingegneria Civile, Design, Edilizia e Ambiente, Seconda Università di Napoli.

Keywords: Preferential flow; Dual-permeability model; local factor of safety; slope stability analysis; COMSOL Multiphysics

1. Introduction

Physically-based subsurface hydrological models can be integrated with slope stability analysis methods to evaluate the risk of landslides in unstable slope. Landslide forecasting is difficult because of the complexity of hydrological processes and failure mechanisms in natural slopes1_{. For example, water flow and solute transport in}

natural slopes are controlled by complicated pore structures and geometries2_{. Chains of (partially connected)}

preferential flow paths are commonly found in various types of slopes, including forest soil and semi-arid land2,3_.

* Corresponding author. Tel.: +31-(0)684470660; fax: +31-(0)152781029.

E-mail address: w.shao@tudelft.nl

Selection and peer-review under responsibility of Dipartimento di Ingegneria Civile, Design, Edilizia e Ambiente, Seconda Università di Napoli.

(2)

Preferential flow through macropores, fractures, and other local high-permeability zones will be relatively fast and influence the subsurface processes2,3,4,5_{. Therefore, increasingly sophisticated theories and models have been}

developed to simulate preferential flow in various environmental systems, such as dual-porosity/dual-permeability models, multi-permeability models, and empirical models6,7_{. The dual-permeability model is the most widely used}

preferential flow model because of its clear physical concept and powerful simulating ability. It assumes that the soil consists of two interacting and overlapping pore domains. In each representative element volume, micropore and preferential flow paths are evenly distributed, which means that the preferential flow domain does not have a specific structure or direction. The fast-flow domain represents the preferential flow paths such as macropores, fractures, fissures, and cracks; the matrix domain represents the soil micropores with relatively low permeability. Water is exchanged between the two domains continuously6,7,8_.

The accuracy of using hydrological-stability models for rainfall-induced landslide forecasting relies on the identification of realistic landslide triggering mechanisms and the correct mathematical description of these mechanisms9,10_{. Preferential flow in the unsaturated soil through cracks, macropores, fissures, pipes etc., has a}

significant influence on pore water pressure distribution and slope stability10,11,12,13_{. There are numerous examples}

and applications of hydrological phenomena that cannot be described with the single permeability Richards’ and Darcian flow concept7,14,15_{. In contrast, most hydro-geotechnical software packages that are used for coupled}

hydrological process and soil mechanics process modeling, such as PLAXIS15_{and FLAC}16_{, are based on a}

single-permeability formulation for subsurface flow modeling. This is appropriate for many practical geotechnical analyses, which are the core use of these commercial models, but is less appropriate for detailed process understanding and quantification, for which software packages need to be used that can handle more complex conceptualizations.

In this study, COMSOL Multiphysics17_{is used to develop a coupled hydrological and slope stability model for}

both single-permeability and dual-permeability subsurface flow. The aim of this paper is to share our experience with developing a dual-permeability model for unsaturated and saturated conditions and couple that with a slope stability model in COMSOL Multiphysics. The theories of the hydro-mechanical framework are introduced first, including the coupled dual-permeability model and local factor of safety (LFS) slope stability analysis method. We present two benchmark examples to demonstrate the accuracy of our COMSOL model. Finally, an example application is presented to quantify the influence of preferential flow on slope stability at the hillslope scale, and the pros and cons of this approach are discussed.

2. Methods

2.1 Introduction of COMSOL Multiphysics software

COMSOL is a Multiphysics modeling tool that solves various coupled physical problems based on Finite Element Analysis and Partial Differential Equations. COMSOL provides a user-friendly interface for mesh generation, equations configuration, and results visualization. The model can be built with user-defined equations or discipline-specific modules (including stress/strain, heat transfer, computational fluid dynamics, and electromagnetics modules). The Subsurface Flow Module and Solid Mechanics Module can be applied for modeling the coupled hydrological and soil mechanical system, respectively. Implementation of a user-defined general mixed (Cauchy) boundary condition in the Subsurface Flow Module provides the possibility to describe variably-saturated flow in a hillslope. The Solid Mechanics Module can be used to describe various types of elastic and plastic stress-strain analyses in soils. Although COMSOL is a potential tool for coupled hydro-mechanical system modeling, the coupling between a dual-permeability model and a soil mechanics model has not been reported in the literature.

2.2 Richards’ Equation Interface

In COMSOL, the Richards’ equation is incorporated in the module for subsurface flow in single-permeability models. The relationships between unsaturated hydraulic conductivity, soil moisture, and pressure head can be described by the Brooks-Corey Model18_{, the van Genuchten Model}19_{, or a user-defined equation. The Richards’}

(3)

pressure head and flux boundary conditions, combined with the user-defined mixed (Cauchy) boundary condition, can be used to describe variably-saturated subsurface flow under the influence of rainfall and evaporation20.

2.3 Dual-permeability model

The dual-permeability model assumes that the matrix flow domain and preferential flow domain co-exist and overlap in the same representative elementary volume, and thus are homogeneously distributed over the entire modelling network. The coupling between matrix flow and preferential flow in the dual-permeability model is achieved by a water exchange term. The water exchange rate is calculated based on the arithmetic mean of the hydraulic conductivities to gain an exchange rate that is dependent on the moisture and pressure states of both domains21_{. Mathematically, the two Richards’ equations combined with a user-defined water exchange function}

constitute a coupled system of partial differential equations8_:

>

@

( ) ( ) ( ) 1 ( ) ( ) ( ) 1 ( ) ( ) ( ) 2 f w f f e f s f f f f m w m m e m s m m m m f f m m w w f m h C h S h S K h h t w h C h S h S K h h t w K h K h h h D w _ª _º * ª º °¬ ¼ _w ¬ ¼ ° ° _w _* ° ª º ® _w ¬ ¼ ° ° °* °¯ (1)

where the subscript f indicates the preferential flow domain, the subscript m indicates the matrix domain, t is time (T), T is the water content (L3_L-3

), C is the differential water capacity (dT/dh) (L-1_{), h is the pressure head (L), K is}

the unsaturated hydraulic conductivity (LT-1_{), S}

e is the effective saturation, Ss is the specific storage, w is the volume fraction of the preferential flow domain or the matrix domain (-), *w is the water exchange term (T-1), and Įw is the

water transfer coefficient (L-2_).

In COMSOL, the dual-permeability model is defined as two coupled Richards’ equations. The coupling is achieved by defining a Mass Source term in the Richards’ Equation Interface, after which the numerical solver recognizes the system of coupled equations during the computations. The distribution of the surface boundary flux between the two domains is defined as the volume ratios of both domains following22_{. Furthermore, a control}

equation is defined to switch from a flux-specified boundary condition to a head-specified condition when ponding occurs.

2.4 Local factor of safety slope stability analysis method

The slope stability analysis is based on the local factor of safety approach. The stress distribution is calculated with a linear elastic model, which is governed by a momentum balance equation:

( ) J T( ) 0

ı b (2) where V is the stress tensor (ML-1_T-2

), b is the unit vector of body forces, and J is the bulk unit weight (ML-2_T-2_),

which is a function of the water content T. The Solid Mechanics Module in COMSOL is used for the stress-strain calculation.

The influence of pore water pressure variation on the suction and effective stress is evaluated with Bishop’s effective stress23,24_:

' ₍ ₎ ₍ ₎

a a w

u u u

V

F

(3) where Vǯ is effective stress, ua is the pore air pressure, uw is the pore water pressure, and F is the matrix suction coefficient which varies from 0 to 1 depending on the saturation. Finally, the local factor of safety (FLFS)is defined as “the ratio of the Coulomb stress at the current state of stress to the Coulomb stress of the potential failure state under the Mohr-Coulomb criterion”25_:

(4)

' ' ' ' 1 3 ' ' ' 1 3 2 cos tan 2 LFS F I c V V I V V ª º « » _¬ _¼ (4) where the cǯ is the cohesion of the slope material (ML-1_T-2

), Iǯ is the friction angle, and V1ǯ and V3ǯ are the first and

third effective stress for the variable saturated soil.

2.5 Coupled dual-permeability and soil mechanics model

The hydrological model and the soil mechanics model are coupled in two ways (Fig.1). First, the body load function depends on the moisture variation, defined in the Body Load Term of the Solid Mechanics Module to manifest the effect of water content variations on the self-weight and stress distribution. Second, the effective stress is determined by the hydrological condition. In the single-permeability model, the pore water pressure result is used directly for effective stress computation and slope stability analysis. Since the shear strength of preferential flow paths is weaker than the surrounding soil matrix, slope failure is determined by the interconnected fissure networks. Therefore, the pore water pressure of the preferential flow domain is used for the coupled dual-permeability and slope stability analysis.

Fig. 1. Structure of coupled dual-permeability model and soil mechanics model

The modelling consists of three sequential steps. First, a stable pore water pressure distribution is obtained from a sufficient warming-up period with a small effective rainfall intensity, and is used as the initial condition for the event rainfall-infiltration study. Second, a variable intensity of event rainfall is used as the input for a transient subsurface flow simulation, while evaporation and transpiration are not being considered. Third, the transient soil moisture and pore water pressure distributions are used for slope stability analysis.

3. Benchmark problems

Two benchmark problems are presented to demonstrate the ability, accuracy and effectiveness of COMSOL software to model subsurface flow and slope stability analysis.

3.1 Furrow irrigation with single-permeability model

The first benchmark is a Furrow irrigation in a single-permeability system with uniform hydraulic conductivity26_{. The study area is 1.5 m deep and 2.0 m wide, with two 20 cm deep furrows. Hydrostatic pressure}

conditions are used as an initial condition. At t=0, the furrows are filled with 10 cm of water. The evolving pressure head profiles at 10 min, 6 hr and 12 hr are shown in Fig. 2.

The top part of Fig. 2 shows the originally published results, and the lower part of Fig. 2 shows the results of the COMSOL model. The results of the two models are very similar, which means COMSOL can be effectively used for subsurface hydrological modeling.

(5)

Fig. 2. Results of furrow irrigation in a single domain system from Vogel [26] (top row) and our COMSOL model (bottom row) 3.2 Local factor of safety calculation

The second benchmark is the local factor of safety (LFS) distribution for a hillslope25. Consider a slope with an inclination of 30q, friction angle of 30q, Poisson's ratio of 1/3, and cohesion of J+/10 (H is the slope height). The boundaries are far away from the study area, the vertical sides and bottom are roller boundary conditions, and the top side and slope are free boundary conditions. The local factor of safety results are shown in Fig. 3.

Fig. 3. The local factor of safety distribution obtained by Lu et al [25] (left) and with our COMSOL model (right)

Fig.3 shows the original published LFS results (left panel), and the solution from user-defined LFS functions in COMSOL (right panel). The COMSOL solution is smooth and stable, and provides detailed information of the factor of safety distribution. This comparison with the results of Lu et al. shows that our COMSOL model can accurately calculate the local factor of safety distribution by a series of user-defined equations.

4. Example of coupled dual-permeability and slope stability model

In this study, a coupled dual-permeability hydrological and slope stability model is developed in COMSOL to quantify the influence of preferential flow on slope stability by comparing the factor of safety results from a single-permeability model and a dual-single-permeability model. Detailed analysis results are based on a series of numerical experiments. In this paper, one straightforward example of a slope with rainfall infiltration is presented to demonstrate the capability of the COMSOL software in developing such models.

Consider a slope consisting of two layers. The slope is 6 meters high and 15 meters wide. The upper layer is 2 meters in depth and consists of sandy loam. The lower layer is clay (Fig.4). A large computational area was defined to reduce the influence of boundary conditions on hydrological and slope stability results. The total computational area is 57 m wide and 33 m high. The mesh generation interface of COMSOL can generate a mesh with different densities according to user specifications. The upper area is very sensitive to rainfall input, while the pressure head in the lower layer will not change as dynamically. Therefore, the mesh density of the upper layer is approximately

(6)

0.2 m in the vertical direction, and 0.5 m in the horizontal direction. The mesh of the lower part is varied to reduce the computation time. The variable density of the mesh increases the computation efficiency without adversely affecting the numerical accuracy. Fig. 4 shows the computational mesh of the study area and the computational area.

Fig. 4. The geometry of computation area and mesh generation results

Table 1 summarizes the parameter settings. The Brooks-Corey equation is used for the soil water retention curve. The volumetric percentage of the preferential flow domain is 10% of the total volume. For the upper layer, the permeability of the preferential flow domain is 100 times larger than for the matrix domain. In the lower layer, the ratio Ks(f) / Ks(m) is equal to 5. The overall saturated hydraulic conductivity (weighted mean of matrix and preferential flow permeabilities) of the dual-permeability model is equal to the saturated hydraulic conductivity of the single-permeability model. The soil mechanics parameters in the single-permeability model and dual-permeability model are all the same: cohesion is 5 kPa, friction angle is 35° and Poisson’s ratio is 0.35.

Table 1 Parameters for Brooks-Corey model

Soil type Basic parameter Saturated hydraulic conductivity (cm/day) Brooks-Corey fitting parameter Tr Ts Ks Ks(f) Ks(m) aBC(1/m) nBC lBC

Upper layer(Sandy Loam) 0.041 0.412 62.16 570.2 5.702 6.8 0.322 1 Lower layer(Clay) 0.09 0.385 1.44 5.143 1.028 2.68 0.131 1 Note: in table 1,Tr is residual water content, Ts is saturated water content, Ks is saturated hydraulic conductivity of single-permeability model,

Ks(f) and Ks(m) is saturated hydraulic conductivity of preferential flow domain and matrix domain in dual-permeability model, the overall

saturated hydraulic conductivity is equal to w(f) Ks(f) + w(m) Ks(m), aBC, nBC and lBC are fitting parameters of Brooks-Corey model.

In COMSOL, MUMPS (multifrontal massively parallel sparse direct solver) is chosen to obtain the numerical solution. The benefits of MUMPS are that it can save time and memory. The drawback of this solver is that it is more sensitive to the mesh quality. For hydrological computations, the mesh density is high enough to get stable and accurate solutions with MUMPS. The error tolerance when iterating is set to 0.001 to ensure computational efficiency and an accurate solution.

A rainfall of 20 mm/hr is applied for 15 hr. The computation is carried out on a desktop computer with a four Core CPU. The consumed time for the single-permeability model is relatively small, about 5 hr on a desktop computer, while for the dual-permeability model, the computation time is about 30 hr. The initial distribution of pore water pressure and local factor of safety are almost the same for both the single-permeability model and dual-permeability model. After rainfall, the water content and slope stability results are shown in Fig.5 and Fig.6.

(7)

Fig. 5. Water content distribution att=15hr,

Fig. 6. Local factor of safety distribution at t=15hr,

The results show that the single-permeability model can only simulate the regular wetting front which is almost parallel to the slope surface. In contrast, the dual-permeability model with a homogeneous distribution of 10% preferential flow paths in the entire domain can simulate the fast wetting front propagation caused by preferential flow. Water can reach the deep part of the soil matrix by preferential flow, instead of the matrix wetting front parallel to the surface. After a sufficient amount of rainfall, the preferential flow causes a perched groundwater table, which significantly reduces the factor of safety. The slope stability analysis results of the dual-permeability model show that the minimum local factor of safety is smaller, and the failure area is larger compared with the single-permeability model results. Under the same amount of rainfall input, the single-single-permeability model overestimates the slope stability condition as compared to the dual-permeability model.

The hydrological results manifest that the dual-permeability model developed with the COMSOL software has the capability of simulating complex subsurface flow processes, such as non-equilibrium flow and bypass flow. The coupled dual-permeability and slope stability model were used successfully for the slope stability analysis under the influence of preferential flow.

5. Conclusions

The benchmark results show that our COMSOL model is capable of simulating flow governed by Richards’ equation and the local factor of safety analysis in a hillslope profile. From the numerical experiments, the coupled dual-permeability and slope stability model can quantify the influence of preferential flow on slope stability, and provide a deeper understanding of slope failure initiation and evolution. The COMSOL interface for mesh generation, equations configuration, results visualization, and discipline-specific modules are easy to use, the coupling between partial differential equations can be achieved easily, and the numerical solver is efficient and stable. In conclusion, the COMSOL software is a flexible platform to build novel models for coupled physical processes. The planned future work includes the extension of the dual-permeability model to simulate complex hydrological systems under the influence of evaporation, transpiration, interception, surface runoff, etc.

References

1. Tsutsumi D, Fujita M. Relative importance of slope material properties and timing of rainfall for the occurrence of landslides. Int J Erosion

(8)

2. Jarvis NJ. A review of non-equilibrium water flow and solute transport in soil macropores: principles, controlling factors and consequences for water quality. Eur J Soil Sci 2007;58(3):523-546.

3. Uchida T, Kosugi KI, Mizuyama T. Effects of pipeflow on hydrological process and its relation to landslide: a review of pipeflow studies in forested headwater catchments. Hydrol Proc 2001;15(11):2151-2174.

4. Hencher SR. Preferential flow paths through soil and rock and their association with landslides. Hydrol Proc 2010;24(12):1610-1630. 5. Krzeminska D, Bogaard T, Malet J-P, v Beek L. A model of hydrological and mechanical feedbacks of preferential fissure flow in a

slow-moving landslide. Hydrol Earth Syst Sci 2013;17(3):947-959.

6. ŠimĤnek, J, Jarvis NJ, van Genuchten MT, Gärdenäs A. Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone, J Hydrol 2003;272(1–4):14-35.

7. Köhne JM, Köhne S, ŠimĤnek J. A review of model applications for structured soils: a) Water flow and tracer transport. J Cont Hydrol 2009;104(1–4):4-35.

8. Gerke HH, van Genuchten MT. Evaluation of a first-order water transfer term for variably saturated dual-porosity flow models. Water Resour

Res 1993;29(4):1225-1238.

9. Coe JA, Michael JA, Crovelli RA, Savage WZ, Laprade WT, Nashem WD. Probabilistic assessment of precipitation-triggered landslides using historical records of landslide occurrence, Seattle, Washington. Environ Eng Geosci 2004;10(2):103-122.

10. van Asch TWJ, Van Beek LPH, Bogaard TA.(2009), The diversity in hydrological triggering systems of landslides. In: Picarelli L, Tommasi P, Urciuoli G, Versace P, editors. Rainfall-Induced Landslides. Napoli: Officine Grafiche F. Giannini & Figli; 2009. Vol. 1 p. 151-156. 11. Bogaard T. Analysis of hydrological processes in unstable clayey slopes. Koninklijk Nederlands Aardrijkskundig Genootschap: Faculteit

Ruimtelijke Wetenschappen, Universiteit Utrecht, Utrecht; 2001.

12. Gerscovich DMS, Vargas Jr EA, de Campos TMP. On the evaluation of unsaturated flow in a natural slope in Rio de Janeiro, Brazil. Eng Geo 2006;88(1–2):23-40.

13. Krzeminska D, Bogaard TA, van Asch TWJ, van Beek LPH. A conceptual model of the hydrological influence of fissures on landslide activity. Hydrol Earth Syst Sci 2012;16(6):1561-1576.

14. Vargas Jr EA, Velloso RC, de Campos TMP, Costa Filho LM. Saturated-unsaturated analysis of water flow in slopes of Rio de Janeiro, Brazil. Comput Geotech 1990;10(3):247-261.

15. Brinkgreve R, Engin E, Swolfs W, Waterman D, Chesaru A, Bonnier P, Galavi V. PLAXIS 3D. 2010. 16. Itasca F. Fast Lagrangian Analysis of Continua. Version 4.0 User’s Guide. Itasca Consulting Group Inc; 2002. 17. Comsol A. COMSOL multiphysics user’s guide. Version July; 2013

18. Brooks RH, Corey AT. Hydraulic properties of porous media. Hydrology Papers, Colorado State University (March); 1964.

19. Van Genuchten MT. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci Soc Am J 1980;44(5):892-898.

20. Chui T, Freyberg D. Implementing Hydrologic Boundary Conditions in a Multiphysics Model. J Hydrol Eng 2009;14(12):1374-1377. 21. Laine-Kaulio H. Development and analysis of a dual-permeability model for subsurface stormflow and solute transport in a forested

hillslope. Aalto University; 1979

22. Dusek J, Gerke HH, Vogel T. Surface Boundary Conditions in Two-Dimensional Dual-Permeability Modeling of Tile Drain Bromide Leaching. Vadose Zone J 2008;7(4):1287-1301.

23. Bishop AW. The principles of effective stress, Norges Geotekniske Institutt; 1960.

24. Griffiths D, Lu N. Unsaturated slope stability analysis with steady infiltration or evaporation using elasto-plastic finite elements. Int J Num

Anal Meth Geomech 2005;29(3):249-267.

25. Lu N, ùener-Kaya B, Wayllace A, Godt JW. Analysis of rainfall-induced slope instability using a field of local factor of safety. Water Resour

Res 2012;48(9):W09524.

26. Vogel T, Gerke HH, Zhang R, Van Genuchten MT. Modeling flow and transport in a two-dimensional dual-permeability system with spatially variable hydraulic properties. J Hydrol 2000;238(1–2):78-89.