A 0-D (box type) mathematical flushing model and a three-dimensional baroclinic numerical model have been presented that are used to simulate transport processes in coastal waters. The numerical model consists of hydrodynamic, transport and turbulence model components. In the hydrodynamic model component, the Navier-Stokes equations are solved with the Boussinesq approximation. The transport model component consists of the pollutant transport model and the water temperature and salinity transport models. In this component, the three dimensional convective diffusion equations are solved for each of the three quantities. In the turbulence model, a two-equation k-epsilon formulation is solved to calculate the kinetic energy of the turbulence and its rate of dissipation, which provides the variable vertical turbulent eddy viscosity. Horizontal eddy viscosity can be simulated by the Smagorinsky algebraic sub-grid scale turbulence model. The, solution method is a composite finite difference-finite element method. In the horizontal plane finite difference approximations and in the vertical plane finite element shape functions are used. The governing equations are solved implicitly in the Cartesian coordinate system. The horizontal mesh sizes can be variable. To increase the vertical resolution, grid clustering can be applied. In the treatment of coastal land boundaries, the flooding and drying processes can be considered. The developed numerical model predictions are verified by the hydraulic model studies conducted on the forced flushing of marinas in enclosed seas.