Abstarct: Contamination of the aquifer systems by organic pollutants is a widespread environmental concern. Hence, the management of the groundwater resources has been strongly emphasized in recent years. Accurate methods are required to predict contaminants transport through porous media in order to determine contamination plumes in the subsurface aquifers and control pollution expansion. At present, oil has been revealed as a laxyer on the groundwater surface in the Tehran oil refinery and nearby areas. Field investigations have shown groundwater and soil pollution in a large area. Furthermore, oil has been measured in monitoring and other wells in the study area. The investigations have revealed the presence of more than 1440000 cubic meter of oil throughout the Tehran oil refinery and the Rey industrial area. In this research, comprehensive models have been developed for each zone of the contaminated area. The results of these models can be used for designing an appropriate remediation strategy. To achieve this goal, a conceptual model was first developed to consider effective mechanisms involved in the reactive transport of hydrocarbon pollutants, pollution plume expansion and hydro-chemical reactions. A two-dimensional model has been then presented by modifying a multi-purpose commercial software called PHOENICS to predict multi-species and multi-phases reactive transport of hydrocarbons through the saturated zone of groundwater flow system in the Rey industrial area. PHOENICS solves model equations numerically utilising a finite volume technique. All non-linear terms of the partial differential equations were handled using PHOENICS by coding in PIL (PHOENICS Input Language), and saving in a Q1-file with calls to the GROUND routine for all non-standard computations. The required coding was added through GROUND routine with FORTRAN99 language. The accuracy of the model was verified with analytical solutions to simplified systems and single processes together with the results obtained by different scenarios of published numerical models by other researchers to more complicated cases. The model was then expanded to predict multi-species and multi-phases reactive transport of hydrocarbon pollutants in the study area. The proposed model was calibrated with the field data for a time period of 36 years; showing a close agreement. It was then implemented to predict reactive transport of hydrocarbons for 15 years after calibration period taking different scenarios into account. The most important characteristics of the model presented in this research are that in addition to the physical transport mechanisms, it incorporates variable source, linear adsorption process and biological kinetics models including Monod, Michaelis-Menten, first- and zero- order mathematical exxpressions. A sensitivity analysis was carried out to assess uncertainties of the model presented in this research and to consider the influence of parameters variability on the results of model. For this analysis, the sensitivity of the model to changes in the dispersivity, distribution coefficient, parameters of Monod, Michaelis-Menten, first- and zero- order kinetics modes on the contaminant plumes were examined by performing several simulations. It was found that the model is most sensitive to changes in dispersivities, distribution coefficient, and parameters of Michaelis-Menten, first- and zero- order kinetics model. On the other hand, the model is not sensitive to the parameters of Monod kinetics model. The reason for this insensibility is that degradation is not limited by microbial kinetics in the simulation, but by dispersive mixing. The modeling results show that by removing LNAPL from the surface of groundwater table of the aquifer after 42 years, there will be a significant change in the distribution of dissolved hydrocarbons in the underlying aquifer. However, aquifer remains contaminated during this time. It will be necessary to develop an appropriate remediation strategy baxsed on biological process in the study area. The results obtained from such investigations can be used to monitor the contaminated sites by hydrocarbon pollutants and design methods for controlling pollution and minimizing environmental problems in the groundwater aquifers