The thesis deals with the numerical modeling of nonreactive and nonlinearly sorbing solutes in groundwater and analysis of the effect of heterogeneity resulting from spatial variation of physical and chemical parameters on the transport of solutes. For this purpose, a hy­brid flux corrected transport (FCT) and central difference method based on operator-split approach is developed for advection-dispersion solute transport equation. The advective transport is solved using the FCT technique, while the dispersive transport is solved using a conventional, fully implicit, finite difference scheme. Three FCT methods are developed and extension to multidimensional cases are discussed. The FCT models developed are anlaysed using test problems possessing analyt­ical solutions for one and two dimensional cases, while analysing advection and dispersion dominated transport situations. Different initial and boundary conditions, which mimic the laboratory and field situations are analysed in order to study numerical dispersion, peak cliping and grid orientation. The developed models are tested to study their relative merits and weaknesses for various grid Peclet and Courant numbers. It is observed from the one dimensional results that all the FCT models perform well for continuous solute sources under varying degrees of Courant number restriction. For sharp solute pulses FCT1 and FCT3 methods fail to simulate the fronts for advection dominated situations even for mod­erate Courant numbers. All the FCT models can be extended to multidimensions using a dimensional-split approach while FCT3 can be made fully multidimensional. It is observed that a dimensional-split approach allows use of higher Courant numbers while tracking the fronts accurately for the cases studied. The capability of the FCT2 model is demonstrated in handling situations where the flow is not aligned along the grid direction. It is observed that FCT2 method is devoid of grid orientation error, which is a common problem for many numerical methods based on Cartesian co-ordinate system. The hybrid FCT2 numerical model which is observed to perform better among the three FCT models is extended to model transport of sorbing solutes. The present study analyses the case of nonlinear sorption with a view to extend the model for any reactive transport situation wherein the chemical reactions are nonlinear in nature. A two step approach is adopted in the present study for coupling the partial differential equation gov­erning the transport and the nonlinear algebraic equation governing the equilibrium sorp­tion. The suitability of explicit-implicit (EI - form) formulation for obtaining accurate solution coupling the transport equation with the nonlinear algebraic equation solved using a Newton-Raphson method is demonstrated. The performance of the numerical model is tested for a range of Peclet numbers for modelling self-sharpening and self-smearing con­centration profiles resulting from nonlinear sorption. It is observed that FCT2 model based on this formulation simulates the fronts quite accurately for both advection and dispersion dominated situations. The delay in the solute mobility and additional dispersion are anal­ysed varying the statistical parameters characterising the heterogeneity namely, correlation scale and variance during the transport of solutes and comparisons are drawn with invariant, cases. The impact of dispersion during the heterogeneous transport is discussed.