Regional Ocean Modeling System (ROMS) is used to develop a new three-dimensional hydrodynamic model for the Chesapeake Bay estuary. Hindcast simulations are conducted for 2 years with markedly different annual river discharges and are compared with time series measurements and high-resolution hydrographic data. The model shows skill in reproducing observed temporal variability in sea level height, salinity, and subtidal current. The agreement with observations is better in the normal runoff year 1997 than in the high runoff year 1996. The model qualitatively reproduces the along-channel and cross-channel salinity distributions during low-to-medium runoff periods. However, during high runoff periods it predicts weaker stratification and a more diffuse halocline than shown by observations. This model/data discrepancy is related to the deficiency of turbulent mixing parameterizations in strong stratification. We have experimented with four turbulence closure schemes (Mellor-Yamada/k-kl, k-$varepsilon$, k-ω, and KPP models) in ROMS but found little difference in the model results. However, vertical stratification shows a strong sensitivity to the background diffusivity. The vertical diffusivity inferred from the model is found to be set by the background diffusivity except in the surface and bottom boundary layers where the turbulence schemes produce similar diffusivity distributions. Among the schemes explored, KPP and k-kl scheme with a background diffusivity of 10-5 or 10-6 m2 s-1 provide the best simulations of the Chesapeake Bay estuary. Both the model sensitivity study and model/data comparison highlight the importance of obtaining a more realistic parameterization for turbulence mixing in a strong pycnocline.