A mesoscale Planetary Boundary Layer (PBL) model with a simple turbulence closure scheme based on the turbulence kinetic energy (TKE) equation and the dissipation (ε) equation is used to simulate atmospheric flow over mesoscale topography. Comparative studies with different parameterizations suggest that with a proper closure assumption for turbulence dissipation, the E-model can simulate the circulation induced by the mesoscale topography with results similar to those obtained using the E-ε model. On the other hand, the first-order closure using O'Brien's cubic interpolation for eddy diffusivities (K) generally produces much larger K profiles in the stable and the unstable regions, which is believed to be due to the overprediction of the height of the PBL. All models with the TKE equation yield quite similar ensemble mean fields, which are found to be little sensitive to the closure assumption for turbulence dissipation, though their predicted magnitudes of TKE and K may differ appreciably. A discussion on the diurnal evolution of the mesoscale topography-induced circulation and the spatial variations of the turbulence fluxes in the surface layer is also given based on the E-ε model results.