As mixing plays a dominant role for the physics of an estuary like the Baltic Sea (seasonal heat storage, mixing in channels, deepwater mixing), different mixing parameterizations for use in three-dimensional (3-D) Baltic Sea models are discussed. Within the Swedish regional climate modeling program, SWECLIM, a 3-D coupled ice-ocean model for the Baltic Sea has been coupled with an improved version of the two-equation k-epsilon turbulence model using a corrected dissipation term, flux boundary conditions to include the effect of a turbulence enhanced layer due to breaking surface gravity waves, and a parameterization for breaking internal waves. Results of multiyear simulations are compared with observations. The seasonal thermocline (the main focus of this paper) is simulated satisfactory. During the stagnation period between 1983 and 1993, simulated salinity in the lower layer of the Baltic Sea decreases as observed. Unsolved problems of the k-epsilon model are discussed. To replace the controversial equation for dissipation, the performance of a hierarchy of k models has been tested and compared with the k-epsilon model. In addition, it is shown that the results of the 1-D turbulence submodel depend very much on the dimensionality of the hydrodynamic model. Using the same turbulence parameterization, vertical velocity shear and density gradients are simulated differently in 1-D column models compared to 3-D ocean circulation models. Finally, the impact of two mixing parameterizations on Baltic Sea climate is discussed.