The prediction of moisture and contaminant transport through unsaturated soil to ground water is becoming increasingly important in the fields of hydrology, agriculture and environmental engineering. Computer aided simulation techniques enables one to conduct a series of systematic numerical experiments to analyze flow phenomenon in subsurface hydrology under various physical and chemical processes. The flow movement depends upon medium characteristics, initial and boundary conditions, which reflect, physical processes occurring below the ground. To understand the effects of physical process an efficient and accurate model is needed. Thus the model developed must be able to handle varied initial and boundary conditions. In this regard, infiltration into a very dry soil becomes a very important problem of study. Most of the earlier numerical models developed are concentrated on the development of an efficient algorithm or the modelling of a particular process which govern the flow in unsaturated or saturated-unsaturated homogeneous medium. Not much work has been done on the analysis of variably saturated flow in layered soil medium. Models to simulate unsaturated flow through dry soils, especially through layered soils with varied boundary conditions are very limited. Further, not much studies have been reported in the literature on the prediction of seepage face development and the phreatic surface movement in variably saturated media with layering. These aspects are very important in determining the flow field and the discharge from the domain. A detailed literature review covering above aspects has been made and is reported in this thesis. In the present work, two dimensional numerical models to predict the movement of wetting front in unsaturated domain and the movement of the phreatic surface in homogeneous and layered porous media under various initial and boundary conditions are developed based on finite difference and finite volume techniques. These models can handle flow in both rectangular flow domains and radial flow domains. The initial condition settings include the handling of very dry soil medium without any transformation of the governing equation, handling of infiltration and constant head conditions at the boundaries under steady state as well as transient scenarios. The models are also able to handle various soil moisture characteristics which depict the nonlinear behaviour between hydraulic conductivity, moisture content and pressure head in a soil media. A mixed form of the governing partial differential equation is used in the present study as it leads to better mass conservation. The finite difference model uses a central difference approximation for the space derivatives and an Eulerian backward difference approximation for the time derivative. The fully implicit formulation is solved iteratively using Strongly Implicit Procedure after making Picard approximation for the highly nonlinear coefficients. The process of infiltration into an initially dry soil leads to the development of a steep wetting front. As the finite volume technique is naturally an upwind method, it is expected to play a positive role in modelling such processes accurately. Hence, a finite volume model is also developed by approximating the convective part using a MUSCL approach and a fully implicit central difference method for the diffusive part of the governing equation. The models developed are validated using both experimental data and numerical solutions for problems reported in the literature. The validation problems cover a wide range of physical scenarios such as: infiltration into a very dry soil, infiltration into a dry soil column with gravity drainage, development of water table mound, steady state drainage in a sand filled wedge shaped tank with seepage face development and transient seepage face development in a rectangular domain. Five test problems are used for the validation of the models. The developed models perform very well for the test problems considered, indicating the models' capability in handling such situations. The results obtained by using the present models for simulating flow through highly unsaturated (very dry) soils show that the models perform very well when compared with models which use transformation techniques to handle such problems. The performance of the present models in comparison with the experimental data and numerical models available in the literature show the suitability of the present models in handling such situations. The present models are also used to analyse various types of unsaturated flow problems with varying initial and boundary conditions. The boundary conditions considered are no flow and /or free flow conditions along the left, right and bottom boundaries with infiltration condition along a part of the top boundary. For the various cases considered in the present study, infiltration rate varies from 5 cm/day to 50cm/day through an initially very dry soil of -15000 cm pressure head in homogeneous and layered soils. Different types of soil media considered vary from sandy loam, loam and clay with horizontal and vertical layering of these soils. A total number of 14 cases are analysed. The results are discussed in terms of pressure head variation in the flow domain along with moisture redistribution for all the cases under consideration. It is observed from these studies that the infiltration rate play an important role on the wetting front movement through layered soils depending on the type of layering and the boundary conditions considered. The soil properties of various layers affect the movement of wetting front by changing the direction of movement. Even though the wetting front movement is predominantly vertical, there is a tendency for the wetting front to move in the horizontal direction as it moves from a coarse soil to fine soil. It is also observed that the vertical layering of soils with different hydraulic conductivity helps in redirecting the flow towards the bottom boundary through the neighboring coarser layers. As finite volume method is more suitable for simulating sharp fronts, it is expected to perform better than finite difference method for simulating infiltration into very dry soils. So, a comparison is made between the performance of these two models by using the above test problems. It is observed from these studies that the performance of both the models are same except that the finite volume method takes much more CPU time than the finite difference model. Considering the type of problems tested, it is observed that for modelling unsaturated and saturated-unsaturated flows, finite difference method is better in comparison to finite volume method. It may be due to the predominant diffusive characteristics of the governing equation even while modelling flow through very dry soils. Proper estimation of the seepage height is an important aspect in finding the discharge through the porous medium. It is observed from the literature that the use of a saturated flow model in such situations can lead to an underestimation of the discharge through the porous medium. This effect is more important when dealing with small dimension problems. It is also observed that various parameters which govern the moisture movement through saturated-unsaturated regions affect proper estimation of the seepage face height and there by discharge. Various factors like effect of boundary conditions, type of soil layering, problem dimension and aspect ratio on seepage face development and the associated phreatic surface formation is studied in the present work. It is seen from the present study that the seepage face development is more in rectangular flow domain than in radial flow domain for both homogeneous and layered soils. It is also seen that the seepage face development in rectangular problems are more sensitive than radial flow problems for various factors considered. The seepage height is also influenced by the tail water level. It is seen from the present study that as the tail water level increases the seepage face reduces with no seepage face development for some of the cases studied. This influence is relatively less for radial flow problems. As the length of the domain increases the seepage height decreases. It is seen that for different cases with same horizontal dimension, as the height of the domain increases the seepage face height also increases. This phenomenon is observed for both homogeneous and layered soil medium. The influence of the aspect ratio, which is the ratio of the length to height of the domain indicate that as the aspect ratio increases the seepage height decreases. The type of the soil layering is observed to have a very strong influence on the seepage face development. The study for understanding the effect of soil layering on the development of seepage face and phreatic surface suggest that as the coarseness of the material increases, the phreatic surface become flatter and its steepness increases with the fineness of the soil. The present model is also used for studying the transient phreatic surface movement and the seepage face development. This is studied for homogeneous and layered rectangular soil medium. The present study is used to understand the effect of specific storage on the phreatic surface movement and the seepage face development. The studies indicate that the influence of specific storage on the seepage face development is insignificant in homogeneous soils with only very little effect in the early time for longer domains. It is also observed that the influence of the specific storage is significant in the case of layered soils. This effect depends on the type of layering and the problem dimension and is observed to have influence for relatively longer period. This observation suggests the importance of specific storage on transient seepage face development. When the specific storage effect is considered the drainage of the soil become faster resulting in a faster decline of the phreatic surface with time. The influence of specific storage is also studied considering the problem dimension effect. It is seen that as the aspect ratio increases, the effect of specific storage on the phreatic surface development decreases. The studies with change in the upstream boundary condition from a constant head to a no flow condition indicate that the effect of specific storage has no significant influence on the phreatic surface development for both homogeneous and layered soils.