A three-dimensional numerical model with free water surface was established to investigate flow characteristics of surface buoyant jets and river plumes. Turbulent shear stresses and turbulent buoyancy fluxes were expressed in terms of the eddy viscosities and diffusivities. Stable stratification effects due to density difference between discharged water and receiving ambient water were taken into with empirical formulae. Through a comparison of numerical results with published experimental data the validity of the model was shown and the optimal stratification functions was determined The three-dimensional spreading characteristics were examined and the effects of inlet densimetric Froude number, inlet aspect ratio and water surface elevation on the flow development were discussed.