We develop a new vertically implicit transport solver, based on two total variation diminishing (TVD) limiters in space and time, inside a 3D unstructured-grid model (SCHISM), and apply it to the Upper Chesapeake Bay (UCB), which has complex geometry and sharp pycnocline. We show that the model is able to accurately and efficiently capture the elevation, velocity, salinity and temperature in both the deep and shallow regions of UCB. Compared with all available CTD casts, the overall model skills have the mean absolute error of 1.08 PSU and 0.85 degrees C, and correlation coefficient of 0.97 and 0.99 for salinity and temperature respectively. More importantly, the new implicit solver better captures the density stratification, which has great implications on biogeochemistry in this estuarine system. The cross-scale capability of the model is demonstrated by extending the high-resolution grids into a tributary (Chester River) and its sub-tributary (Corsica River), with minimal impact on the model efficiency. The model is also able to capture complex 3D structures at the transition zone between the main bay and the tributary, including the three-layered circulation in Baltimore Harbor. As more and more attention is being paid to the productive shallows in the Chesapeake Bay and other estuaries, the model can serve as a very powerful management tool to understand the impact of both local and remote forcing functions. (C) 2016 Elsevier Ltd. All rights reserved.