Land subsidence is one of th e growing, crucial challenges faced in many great cities located in areas containing fine grained deposits. This phenomenon causes problems such as damaging buildings, roads and underground sewerage system. Recently, many areas affected by land subsidence have been observed in some arid regions in Iran including Tehran, Rafsanjan and the Khorasan plains. The main cause leading to this phenomenon is the excessive extraction of groundwater, which leads to drop in groundwater head. The process of land subsiden ce due to change in groundwater level can be analyzed using the theory of consolidation. This means that a drop in the groundwater level results in an increase in the vertical effective stresses which, in turn, causes land subsidence. Moreover, changing groundwater level and increase in effective stress causes changes in some soil properties. In some soil layers, saturated soils turn into unsaturated soils and exhibit different behavior. In this study, the behavior of unsaturated soils is used to study and analyze land subsidence. Observations and measurements on land subsidence in the plains west of the city of Tehran are first studied and then compared with results of numerical modeling using the PLAXIS 3D software. These comparisons show reasonable agreem ent of predicted and measured subsidence. The calibrated numerical model is then used to predict future land subsidence in the same areas. Keywords: Land subsidence, unsaturated soil, numerical simulation, groundwater withdrawal 1 INTRODUCTION Land subsidence can be defined as a gradual settling or a rapid sinking of discrete segments of earth surface (Galloway and Burbey 2011) . This phenomenon may occur due to various processes. Geological processes such as solubility, thawing of ice, vibration, and densification of deposits are some of the natural occurrence s that can lead to land subsidence (Waltham 1989) . Human activities such as undergrou nd mining and withdrawal of gas, oil, and groundwater are other cause s of land subsidence (ZHU et al. 2005) . Land subsidence can lead to increased susceptibility to flooding, structural damages (buildings, roads, and highwa ys, water pipes) and ground fractures (Budhu and Adiyaman 2010) . Recently, the growth in population and increas e in the demand for water resources has led to excessive exploitation of groundwater (Li, Qian, and Wu 2014) . Therefore , groundwater withdrawal has b ecome the primary cause of land subsidence in most areas of the world. Land subsidence due to pumping groundwater has been reported in Jakarta and Samarang, Indonesia (Chaussard et al. 2013) , Shanghai, China (Hu et al. 2004) , Beijing, China (Ng et al. 2012) , and Mexico City, Mexico (Yan et al. 2012) . In Iran, land subsidence has been reported in Rafsanjan (Mousavi et al. 2001) , Tehran (Mahmoudpour et al. 2016) and Mashhad (Motagh et al. 2007) The process of land subsidence due to cha nge in groundwater level can be analyzed using the theory of consolidation , first formulated by Biot (1941). In this respect, the drop in groundwater level increases effective stress which, in turn, causes land subsidence . In this paper, the behav ior of un saturated soil s is used to study land subsidence occurred in the Southwester n plain of Tehran due to the withdrawal of groundwater from October 1984 to October 2011. The results indicate that considering the behavior of unsaturated soils in simulating the process of land subsidence results in a higher amount of settlement compared to considering only the behavior of saturated soils. Then, three possible scenarios are discussed for predicting land subsiden ce of the plain in October 2030: 1) the case that gro undwater level will be constant till 2030, 2) the case that the groundwater level continues to decline with the same trend after October 2011, 3) the case that the ground water will decline at a lower rate after October 2011. 2 SITE DESCRIPTION The Tehran basin covers an approximate area of 2250 km2 and is located in the Southwestern part of the Tehran province (“Fig.1”). Land subsidence has been reported mainly in the southwestern part of the Tehran basin and is attributed to the withdrawal of gro undwater (Mahmoudpour et al. 2016) . In the study area, results of six boreholes are available , and the y indicate that the stratigraphy of the Tehran southwestern plain consists of three aquifers and three aquitards. Also, in this area , there are some piezometric wells which can help in the determination of land subs idence and change in the ground water level for discrete areas of the plain. According to “Fig.1” , the location of piezometric well no. “P13” is approximately close to borehole “BH2”. Therefor e, based on the observation results obtained from “P13”, groundwater level change s are assumed in the location of “BH2”. Soil properties and thickness of each layer in “BH2” are presented in Table 1 . (Mahmoudpour et al. 2016) . Change s in groundwater l evel over a period of 27 years for the Tehran plain from October 1984 to October 2011 are shown in “Fig.2”. According to “Fig.2”, the groundwater level in the plain has decrease d about 11.36 m from October 1994 to October 2011 (Mahmoudpour et al. 2016). Fig. 1. Site location map. Table 1. Properties of soil at “ BH2 ” Hydrological units Lithology Thickness (m) Unit weight, γ (kN/m3) Void ratio, e Elastic modulus (KN/m2) Compression index, C c Recompression Index, Cr Vertical hydraulic conductivity, Kv (m/s) Horizontal hydraulic conductivity, Kh (m/s) Aquitard 1 Silt and clay 16.5 18 0.66 - 0.305 0.101 3.4×10-10 5.1×10-10 Aquifer 1 Silty sand 3.2 16.5 0.28 344.7×103 - - 3.86×10-6 4.83×10-6 Aquitard 2 Silt and clay 17.3 16.55 0.6 - 0.25 0.08 2.2×10-10 3.3×10-10 Aquifer 2 Clayey sand 14 19.8 0.29 2.79×106 - - 9.36×10-6 1.17×10--5 Aquitard 3 Silt and clay 41.5 20.52 0.73 - 0.33 0.11 1.4×10-10 2.09×10-10 Aquifer 3 Silty sand 8 21.77 0.29 5.61×106 - - 2.08×10--5 2.61×10-5 3 FEM MODEL Numerical modeling i s performed for evaluating land subsidence due to the withdrawal of groundwater level in the Southwestern of Tehran plain. The commercial finite element code PLAXIS 3D 2017 is used for numerical analyses. PLAXIS 3D , by implementing fully -coupled deformatio n analysis , helps simulate Land subsidence in saturated and unsaturated soils. For modeling of soil layers, Tetrahedral iso-parametric 10 -node volumetric elements with 4 Gaussian points were used. In the simulation , the aquitards were represented by the hardening soil model and the aquifers by the linear elastic model. For a nalyzing land subsidence in the plain , the mechanics of both saturated and unsaturated soil should be taken into account. In the layers in which the soils are saturated, the amount of land subsidence is determined using the theory of consolidation. According to Terzaghi, the relation ship between the vertical effective stress (  ), the pore pressure ( p ) in the soils , and the total stress ( ) is as follows (Terzaghi 1951) : p (1) The decli ne in ground water level causes decrease in the pore pressure, and therefore, increas e in the vertical effective stress in the soil. For many types of soil, t he void ratio will decrease due to increas e in the vertical effective stress. Based on the theory of soil consolidation, t he relation ship between vo id ratio and effective stress can be written as follows : 10logcC e for (max)zz zz (2) 10logrC e for (max)zz zz (3) where C c is the compression and Cr is the recompression ind ices, and σ'zz(max) is the preconsolidation stress. The comp ression in the soil layers can be related to the change in the void ratio using equation (4) (Zeitoun and Wakshal 2013) . eebb  00 1 (4) where Δb is the change in layer thickness, b 0 is the initial thickness of the layer , and e0 is the initial void ratio of the soil. For validating the results obtained from the PLAXIS software for the condition in which the groundwater level changes , the unsaturated portion of the soils was first ignored , and a consolidation example was solve d analytically using the analytical relationships of the theory of consolidation , and the result was compare d with that obtained from the PLAXIS software. In the example , the ground was assumed to consist of a sandy layer at the top and a clayey layer at the bottom . For the sandy layer , it was assumed that : height of layer= 10 m, γunsat = 17 kN/m3 and γsat =19 kN/m3, and for the clayey layer: height of the layer= 10 m, γsat =18 kN/m3, OCR = 1.1, Cc=0.3, Cr= 0.03 and e0=0.6. In the initial condition, the groundwater level is assumed to be at the depth of 1m below ground surface. T hen, the groundwater level drops 3m in 50 days . The consolidation settlement of the clayey layer after 5years is calculated using analytical relationships and the results are compared with that obtained from numerical modeling using the PLAXIS software. Analytical solution of the problem resulted in a consolidation settlement of 0.063 m for the clayey layer . “Fig. 3” presents the consolidation settlement of the soil (U z) versus depth ( Z), obtained from PLAXIS , which predicts a settlement of about 0.07 05 m at 10 m depth , at the surface of the clay layer. Therefore , results of PLAXIS are considered to be in favorable agreement with those from the analytical relationships. The reason of the difference between results derived from analytical solving and nu merical modeling is that in Analytical solving, for calculating the settlement of consolidation, Terzaghi’s theory is used w hich has simplified assumption., w hile for analyzing settlement by numerical modeling, fully -coupled deformation is implemented whic h leads to accurate results. Fig. 2. Ground water change in Tehran Plain. Fig. 3. Subsidence of the soil derived from PLAXIS 3.1 Unsaturated soils In areas where water level drops with time, suction above the phreatic level increases, and degre e of saturation of the soil decreases. When the degree of saturation of the soil drops to values below 100 %, the soil will be “unsaturated ,” and its behavior is analyzed in the PLAXIS 3D software using the Bishop effective stress and a suitable hydraulic m odel. 3.2 Bishop effective stress The governing equations of consolidation implemented in PLAXIS follow Biot’s theory (Biot 1941) . The formulation is based on small strain theory and Darcy’s law for liquid flow , in which Bishop’s effect ive stress theory is used. (Bishop and Blight 1963)   a w p pm   1 (5) where is the vector of total stresses,  is the effective stress, m is a vector containing unity terms for normal stress components and zero terms for the shear stress ,  is the matric suction coefficient which varies from 0 to 1 , pw and p a are the pore water pre ssure and the pore air pressure (Galavi 2010) . For practical application s, the pore pressure is small enough to be neglected . The matric suction coefficient is derived experimentally and depends on the degree of saturation, porosity and the matric suction (pw-pa)(Bishop and Blight 1963) . In PL AXIS, the matric suction coefficient is assumed to be equal to the effective saturation defined later in Equation (8) . Therefore , the effective stress relationship is simplified to: aepSm (6) in which eS is the effective saturation (Galavi 2010) . 3.3 Hydraulic model In PLAXIS, five type s of hydraulic models have been utilized, namely Van Genuchten, Mual em, linearized Van Genuch ten, spline and fully saturated (Galavi 2010) . Van Genuchten (Van Genuchten 1980) function relates the saturation to the suction pore pressure head p and describes the hydraulic behavior of unsaturated soils.  cngg p a res sar res p g S S S S    1) ( (7) where gP ww p Sres is the degree of saturation at states in which the water remains in the soil even at high suction s, Ssat is the degree of saturation when the -whole of the voids are filled with water. ga, gn, gc are parame ters derived experimentally. The effective saturation is defined as: res satres eS SSSS (8) The relative permeability according to Mualem – Van Genechten is: 2 1) 1(1)()(      egg g e relnn lS S S K (9) gl is an empirical parameter (Mualem 1976) . 3.4 Model cond ition In order to calculate the land subsidence in the plain, the geological properties of the soil observed at “BH2” are selected, and an area of 40m×40m around the borehole is simulated. (“Fig.4 ”). Geological properties of each layer in the model is shown in Table1 . In PLXIS model, change in the groundwater level was considered from October 1984 to October 2011 by defining a flow function table. “Fig.2” presents the variation s of ground water table in the Tehran plain. Based on recorded data, it is concluded that initial water level in October 1984 was at the depth of 7.36 m below ground surface and it dropped by about 11m by October 2011. In the current numerical simulation, the Van Genuchten hydraulic model is used for analyzing the behavior of the unsaturated soil. This model was assigned to layers which are subjected to ground water changes and developed an unsaturated zone. In this study, it is intended to examine the effe cts of unsaturated zone of the soil on the amount of subsidence and, therefore, the thickness of layers is assumed to be very small in the Van Genuchten hydraulic model such that maximum effects of the unsaturated soil can be determined . Fully -coupled def ormation is implemented in the model as an analyzing method for calculating land subsidence due to change in ground water level. Fully -coupled deformation analysis is used when it is important to consider the simultaneous development of deformations and po re pressure changes in the saturated and unsaturated soils. Fig. 4. PLAXIS model of soil at P13 4 RESULTS AND DISCUSSI ON The model was used for the analy sis of land subsidence in the Southwestern of the Tehran plain, investigating the effects of con sidering u nsaturated soils behavior on the results, and predicting land subsidence in the plain during the time period of 1984 -2030. According to “Fig. 5”, in the period from October 1984 to October 2011 the water level mostly declined. However, the trend o f ground water changes in the mentioned duration is complex and consists of many ups and downs. During 1989 -1991 the groundwater level declined, consequently there is a rap id subsidence during the 1990 -1992 period . From 1991 -1994 , the groundwater level gra dually increased, leading to relative reduc tion of land subsidence from 1992 -1996. Land subsidence accelerated again during 1999 -2005 because of the rapid decline in ground water level during 1998 -2004. Also, “Fig5” shows that there is a time delay of about 3 years in occurring land subsidence , which indicates that subsidence will continue to happen a long time after the drop in groundwater level . In order to evaluate effect s of the trend of ground water level change on the results , land subsidence is analyzed for the condition in which water level change s linearly from 1984 -2011 and also the condition in which groundwater level change s according to measurements in the same period. “Fig.5” compares the results of these two conditions . It can be noticed from “Fig.5” that there is a difference between the amount s Fig. 5. Ground water changes and land subsidence at P13 during 1984 -2011 . Fig. 6. Measured Subsidence rate in the plain and location of “BH2” at the end of 2004. of land subsidence result ing from these conditions during 1984 -2011. Land subsidence due to measured ground water change in October 2011 is about 0.82 m and the land subsidence due to linear decline in ground water is 0.8m. Therefore, the amount of land subsidence obtained from th e two condition, at the end of the time period, is approximately the same. 4.1 Comparison with measurements Land subsidence in the Southwestern of Tehran plain is measured using leveling and GPS data, which were carried out by the National Cartography C enter of Iran, and InSar data (Mahmoudpour et al. 2016) . “Fig.6” shows the map of Tehran plain which contains the contours of subsidence rates at the end of 2004. The location of “BH2” is determined in “Fig.6”, according to which the yearly rate of lan d subsidence in the location of the borehole at the end of 2004 is approximately 8 cm/year. In the PLAXIS model, subsidence in Southwestern Tehran is calculated from October 1984 to October 2011. “ Fig.5” shows the results of the model in this duration . For deriving the rate of subsidence at a specific time, the slope of subsidence curve should be calculated. Based on “Fig.6”, the subsidence at the end of 2004 is about 7.3 cm/year. The comparison of the simulated results to observations shows a good mat ch to measured subsidence rate. 4.2 Model predictions of unsaturated soil effects The model was analyzed under two different conditions. In one condition, the suction is ignored in our calculation and amount of subsidence is evaluated without consider ing unsaturated soil. In another condition, the behavior of unsaturated soil is considered in the analyzing of land subsidence. According to “Fig. 7”, considering the behavior of unsaturated soil in the model generally leads to increase in the amount of sub sidence. For the period of October 1984 to October 2011, when measured groundwater changes was available, the amount of land subsidence for the condition that suction was ignored, was about 0.81m, and by considering suction, it was around 0.94m. The land s ubsidence in Southwestern Tehran is predicted for the period of 1984 -2030 by considering three scenarios. These scenarios regarding ground water changes in the future are presented in “Fig. 7”. In the first scenario, it is assumed that ground water level will be constant after October 2011. “Fig.7” shows that the amount of land subsidence will slightly increase after October 2011. Land subsidence increased after October 2011 and by October 2018, it would be constant. The rate of subsidence after 2018 will be close to zero. In this scenario, the amount of land subsidence in the duration from October 1984 to October 2030, with the consideration of suction change effects, will be about 1.06 m and without considering suction will be 0.91m. Therefore, by conside ring the effect of unsaturated soils, the results will be increased 16% at October 2030. Fig. 7. Predicting land subsidence in for duration of 1984 -2030 . In the second scenario, groundwater level declined with the same trend after October 2011. In th is case, land subsidence will occur continuously until 2030. “Fig.7” shows that the amount of land subsidence at the October 2030, with the consideration of suction will be 1.44 m and without considering unsaturated soils will be 1.31m. Therefore, by consi dering suction, the amounts of measured subsidence will be increased 10%. In the third scenario, the water level will decline with the lower rate of trend after October 2011. Consequently, the amount of subsidence that will happen will be lower in comp arison with the second scenario. In this situation, the amount of land subsidence in October 2030, considering unsaturated soils will be 1.25m and with no suction will be 1.06 m. Therefore, the amounts of measured subsidence at the end of 2030 will be incr eased 17% by considering the effect of unsaturated soils. 5 CONCLUSIONS Land subsidence is one of the growing challenges all over the world. Withdrawal of groundwater is the prime cause of land subsidence. In the Tehran plain, declining of water level has been intensified because of the excessive exploitation of ground water. Therefore, the Tehran plain is faced with a relatively large amount of land subsidence. The PLAXIS software was utilized for modeling an area in the Southwestern Tehran plain. The amo unt of land subsidence derived from simulation showed good correlation with measurement. The effects of unsaturated soil on subsidence were analyzed in the model. The amount of land subsidence in the model that was analyzed with considering the behavior of unsaturated soils, was higher than the model in which suction was ignored. The PLAXIS model was used to predict the amount of land subsidence in the time period of October 1984 to October 2030. Three scenarios were assumed for the change in ground wate r level during this period. According to the first scenario, by considering the effect of unsaturated soils, the results will be increased 16% at October 2030. Based on second scenario and third scenario, by considering the effect of unsaturated soils, the results will be increased 10% and 17%, respectively .