Determination of Matric Suction and Saturation Degree for Unsaturated Soils, Comparative Study - Numerical Method versus Analytical Method

Matric suction is a soil parameter which influences the behaviour of unsaturated soils in both terms of shear strength and permeability. It is a necessary aspect to know the variation of matric suction in unsaturated soil zone for solving geotechnical issues like unsaturated soil slopes stability or bearing capacity for unsaturated foundation ground. Mathematical expression of the dependency between soil moisture content and it’s matric suction (soil water characteristic curve) has a powerful character of nonlinearity. This paper presents two methods to determine the variation of matric suction along the depth included between groundwater level and soil level. First method is an analytical approach to emphasize one direction steady state unsaturated infiltration phenomenon that occurs between the groundwater level and the soil level. There were simulated three different situations in terms of border conditions: precipitations (inflow conditions on ground surface), evaporation (outflow conditions on ground surface), and perfect equilibrium (no flow on ground surface). Numerical method is finite element method used for steady state, two-dimensional, unsaturated infiltration calculus. Regarding boundary conditions there were simulated identical situations as in analytical approach. For both methods, was adopted the equation proposed by van Genuchten-Mualen (1980) for mathematical expression of soil water characteristic curve. Also for the unsaturated soil permeability prediction model was adopted the equation proposed by van Genuchten-Mualen. The fitting parameters of these models were adopted according to RETC 6.02 software in function of soil type. The analyses were performed in both methods for three major soil types: clay, silt and sand. For each soil type were concluded analyses for three situations in terms of border conditions applied on soil surface: inflow, outflow, and no flow. The obtained results are presented in order to highlight the differences/similarities between the methods and the advantages / disadvantages of each one.


Introduction
The behaviour of unsaturated soils is decisively conditioned by the soil suction, ψ. Soil suction is composed of matric suction and osmotic suction. In geotechnical engineering applications the osmotic suction is usually ignored [1,2,3]. Matric suction is defined as difference between pore air pressure, u a , and pore water pressure, u w . In geotechnical engineering field, pore air pressure can be considered equal with atmospheric pressure and therefore it can be ignored (u a =0), [1]. Determination of matric suction distribution in unsaturated zone represents a major step, necessary for further geotechnical applications regarding specific unsaturated soil behaviour. Matric suction (u a -u w ) is necessary for determination of unsaturated soil shear strength, permeability of unsaturated soils, and effective stress state in unsaturated zone [1,2,4,5], [6]. The dependence between soil moisture and matric suction is described by soil water characteristic curve (SWCC). The equation for moisture-matric suction interdependency has a powerful character of nonlinearity. In the literature are found many SWCC equation proposed by different researchers [1,2,6,7,8]. In this paper is adopted the equation proposed by van Genuchten-Mualen (1980), [1]. According with this approach, the volumetric water content, θ, is defined as a function of matric suction, u a -u w , (1 In equation above, θ(u a -u w ) represents volumetric water content expressed in function of matric suction, θ r represents residual volumetric water content, θ s represents saturated volumetric water content, α and n n are function specific fitting parameters respectively m is another secondary fitting parameter: m=1-(1/n n ). Volumetric water content represents the ratio between water volume from soil pores (V w ) and total volume of the soil (V); θ= V w / V. Saturated state (water fill all the pores: V w =V p ) is defined by saturated volumetric water content (θ s ) which is equal to soil porosity (n) : θ s =V p /V=n.
Another parameter which varies in function of matric suction is soil permeability k w (u a -u w ). The permeability function for unsaturated soils is adopted according to van Genuchten-Mualen, using as constants: saturated permeability; k s , and fitting parameters:α, n n and m, (2).
Saturation degree (S r ), represents the ration between soil water quantity (natural state) and maximum soil water quantity (saturated state): S r =θ/θ s = θ/n. For residual water content it corresponds a residual saturation degree: S r.rez = θ r /n. In unsaturated zone, saturation degree can be defined in function of matric suction and fitting parameters: α, n n , m (3).

Analytical approach (AA)
Analytical calculus approach provides results for one-dimension (vertical) steady state unsaturated water infiltration along unsaturated zone [9]. Unsaturated zone has a lower boundary consisting in ground water table (GWT) and the upper boundary as natural ground level (NGL), according to figure 1. Furthermore, we will focus on matric suction variation, ignoring pore air pressure (u a -u w =-u w ), respectively the variation of saturation degree in unsaturated zone. The variation of these parameters will be concluded in perfect equilibrium condition (v p =0), in inflow condition (precipitation, v p <0) and outflow condition (evaporation/dryness, v p >0). All these boundary conditions are concluded by applying an external flow (v p ) on NGL, (figure 1). The starting point of the model is Darcy law for steady state infiltration on one vertical direction for unsaturated soils (4). In equation 4, v z represents vertical water flow, k w represents unsaturated soil permeability, dH/dz is vertical hydraulic gradient, v p represents external flow applied on NGL, H represent hydraulic head expressed as H=z+u/γ w , [9], where z is the level measured to a datum (GWT), u is pore water pressure and γ w is water volumetric weight.

Figure 1. Conceptual sketch of analytical approach
In the case of unsaturated soil, pore water pressure will be equalized to matric suction and the hydraulic head will become H=z+(u a -u w )/γ w . By introducing these terms in Darcy law, (4) was obtained the differential equation of matric suction in relation to z (the distance of a point from unsaturated zone to GWT).
By substitution of (2) in (5) and by notation the matric suction with S(z) we will obtain the differential equation of matric suction using the van Genuchten-Mualen prediction model for unsaturated soil permeability.
In order to find the value of matric suction along the depth of unsaturated zone (h) from z=0 to z=h we will appeal to finite difference method. The applied algorithm implies establishment of an initial limit z init =0 and a final limit z final =h of the variable interval (figure 1). After that the whole interval will be divided in a finite number, N, of intervals (incremental steps), ∆z, ∆z=(z final -z init )/N. For each interval (∆z) will be assigned an index (i), i=0…N, respectively z i = z init +i  Δz. Border condition will be applied: S 0 =0 (matric suction is equal to 0 on GWT) and furthermore we can find the value of matric suction for any point, i, between z init and z final at Δz intervals (7).
By a progressive application of the algorithm we can find the value of matric suction for any point, p, situated in unsaturated zone at z p distance from GWT if z p is a multiple of Δz. S z S   ( 8 ) In equation 8, S(z p ) represents the matric suction at distance z p from GWT and S zp/Δz represents the matric suction according with equation 7. In a similar manner if we use equation 3 we can find the variation of saturation degree S r (z), defined as a z function for any point, p, included in unsaturated zone and situated at z p distance from G.W.T.
For numerical applications from this paper the height of unsaturated zone was considered 4m (h=4m), total number of intervals is N=500 and Δz=0.008m. The points where has been determined the matric suction and the degree of saturation were disposed at 0.2 m intervals along the height of the unsaturated zone. The previous calculus method was implemented in Mathcad 14 software.

Numerical method -Finite element method for 2D ground water flow (FEM)
As an alternative of analytical approach (AA) there was concluded a simulation using finite element method (FEM) delivered by GEO5-FEM (5.0) software [10]. There was concluded 2D, unsaturated steady state infiltration calculus for an isotropic permeable soil according to constitutive equation [11], [10]. In constitutive equation (10), x, y denotes the direction of the 2D domain and H represents hydraulic head, H= z+(u a -u w )/γ w .
Prediction permeability model for unsaturated soils was adopted Genuchten-Mualen model (2). The geometry of analysis domain is presented in figure 2   impervious (v p =0). In second case on the upper edge of domain is applied a vertical inflow (v p <0) in order to stimulate rainfall condition. In the last case on the upper edge of domain is applied a vertical outflow (v p >0) in order to stimulate evaporation (dryness) condition.

Numerical application and result comparison
In order to reveal the occurred difference between analytical approach (AA) and numerical method (FEM) there were concluded computations for 3 different soil types (sand, silt and clay). The necessary parameters for SWCC function (according to Genuchten-Mualen approach) were determined using RETC 6.02 software (Code for Quantifying the Hydraulic Functions of Unsaturated Soils) and are presented in table 1.

Conclusions and discussion
There was performed analysis for each soil type, according to 3 situations regarding applied boundary conditions. Therefore it can be emphasize certain aspects regarding matric suction and saturation degree variation along unsaturated zone height. For sand, the obtained results with both methods analytical approach (AA) and finite element method (FEM) present similar allure for all applied border condition. Nevertheless the variation of saturation degree is smoother (progressive) in case of FEM analysis (figure 4). For silt, obtained differences between calculus methods, exceed in some cases 100% ( figure 6, 7). There were recorded prevailing differences for the degree of saturation. The biggest differences were obtained for clay (table 4, figure 9 and figure 10). It can be observed there were major differences between AA and FEM in situation of rainfall (inflow boundary condition) and evaporation (outflow boundary condition).
However, certain aspects should be mentioned regarding the base assumptions on which AA was built. GWT is considered constant and it represents the origin of the coordinate system for the determination of the searched parameters (u a -u w , S r ). In real condition, GWT is variable in function of climatic conditions (rainfall, evaporation). FEM captures this aspect. Therefore for silt and clay when are applied inflow and outflow boundary conditions (v p ≠0), the control points situated on GWT level record values of matric suction different from zero. For rainfall condition it can be observed the raising of GWT described by positive values of matric suction recorded on GWT control point. For evaporation condition it can be observed a decrease of GWT described by negative values of matric suction recorded on GWT control point (table 3,4). This physical phenomenon cannot be surprised using AA because the points situated on GWT have an imposed 0 value for matric suction. This assumption is necessary for applying finite difference method on differential equation presented as equation 6. The assumption respect physical reality only in condition of perfect equilibrium (with no inflow or outflow).
After all differences, AA has a major disadvantage consisting in inability of taking into account the transient character of water infiltration. Such approach is necessary because inflow and outflow condition represent in reality climatic events described by a non-constant intensity and a finite duration. This fact is not an obstacle for FEM.