This study investigates the role of driving atmospheric forces [winds, net heat flux, and evaporation-precipitation (E-P)] and the possible mechanisms on the mixed layer depth (MLD) spatiotemporal variability in the Bay of Bengal (BoB) using a finer-resolution (similar to 9 km) Regional Ocean Modeling System (ROMS). The model simulation is configured for 2004-2012 with initial and boundary conditions from Mercator Ocean data, atmospheric forcing from ECMWF, and climatological river input. Comparison with in situ observations shows that the model well produces the mean seasonal characteristics of MLD variability. Spatial correlations have been carried out between the atmospheric forces and MLD to determine the relative influence of these forces on the MLD variability. During the southwest monsoon, strong southwesterly wind plays a crucial role in deepening the MLD up to 80 m in the southern BoB. In winter, the decrease in net heat flux and the increase in positive E-P in the upper ocean deepen the MLD in the northern bay to a depth of similar to 60 m. However, the high density stratification in the northern BoB due to the large inflow of freshwater from various rivers limits the seasonality of MLD, particularly in the northern BoB. The positive correlation between E-P and MLD suggests that the advection of the freshwater plume by the southward East India Coastal Current (EICC) forms a thin mixed layer in the western BoB. The model results well captured the contrasting Indian Ocean Dipole (IOD) years and reveal that the interannual variability of MLD seems to be connected to the IOD events in the BoB. The MLD is shallow during the positive phase of IOD (pIOD) in 2006, whereas the mean MLD is more profound (similar to 50 m) during the negative period of IOD (nIOD) in 2010. In pIOD years, the anomalous upwelling coastal Kelvin waves (KWs) propagate, reflect Rossby waves, and trigger upwelling to form a shallow MLD throughout the BoB. The positive net heat flux at the air-sea interface also plays a dominant role in the MLD shoaling in pIOD year, as shortwave radiation increases and exceeds the cooling effect of latent heat flux during this period.