To provide insightful information on water quality management, it is crucial to improve the understanding of the complex biogeochemical cycles of Chesapeake Bay (CB), so a three-dimensional unstructured grid-based water quality model (ICM based on the finite-volume coastal ocean model (FVCOM)) was configured for CB. To fully accommodate the CB study, the water quality simulations were evaluated by using different horizontal and vertical model resolutions, various wind sources and other hydrodynamic and boundary settings. It was found that sufficient horizontal and vertical resolution favored simulating material transport efficiently and that winds from North American Regional Reanalysis (NARR) generated stronger mixing and higher model skill for dissolved oxygen simulation relative to observed winds. Additionally, simulated turbulent mixing was more influential on water quality dynamics than that of bottom friction: the former considerably influenced the summer oxygen ventilation and new primary production, while the latter was found to have little effect on the vertical oxygen exchange. Finally, uncertainties in riverine loading led to larger deviation in nutrient and phytoplankton simulation than that of benthic flux, open boundary loading and predation. Considering these factors, the model showed reasonable skill in simulating water quality dynamics in a 10-year (2003–2012) period and captured the seasonal chlorophyll-a distribution patterns. Overall, this coupled modeling system could be utilized to analyze the spatiotemporal variation of water quality dynamics and to predict their key biophysical drivers in the future.