A three-dimensional primitive equation numerical model was applied to Lake Michigan for the periods 1982--1983 and 1994--1995 to study seasonal and interannual variability of lake-wide circulation and thermal structure in the lake. The model was able to reproduce all of the basic features of the thermal structure in Lake Michigan: spring thermal bar, full stratification, deepening of the thermocline during the fall cooling, and finally, an overturn in the late fall. Large-scale circulation patterns tend to be cyclonic (counterclockwise), with cyclonic circulation within each subbasin. The largest currents and maximum cyclonic vorticity occur in the fall and winter when temperature gradients are low but wind stresses are strongest. The smallest currents and minimum cyclonic vorticity occur in spring and summer when temperature gradients are strong but wind stresses are weakest. All these facts are in agreement with observations. The main shortcoming of the model was that it tended to predict a more diffuse thermocline than was indicated by observations and explained only up to half of the variance observed in horizontal currents at timescales shorter than a day. ¿ 2001 American Geophysical Union