A numerical simulation of groundwater flow and contaminant transport in porous media is presented. Emphasis is given to horizontal subsurface flow constructed wetlands and to evaluation of the optimal range of values for the first-order removal coefficient. The numerical procedure is based on the MODFLOW code and applied for the simulation of horizontal subsurface flow wetlands pilot-scale units, constructed and operated in Democritus University of Thrace, Xanthi, Greece. The effects of the inlet and outlet recharge positions to the optimum contaminant removal are also numerically investigated.