It is widely recognized that Great Lakes of North America is particularly sensitive to climate change and play a critical role in regulating regional climate. The lake-ice-atmosphere interactions in the region lead to large fluctuations in lake thermal structure, ice coverage, and lake water levels. In turn, these changes have dramatic affects on regional hydro-climate due to the massive heat capacity and thermal inertia of the Great Lakes. Accurate presentations of lake-ice-atmosphere interactions in climate modeling remain one of the most critical and unresolved issues. To date, the most advanced representations of the Great Lakes region in regional climate modeling (RCM) are performed with a one-dimensional (1-D) lake model. While some progresses have been made in refining the 1-D deep lake model processes, 1-D lake model is fundamentally incapable of resolving realistic physical processes in the Great Lakes. For the first time, we developed a two-way coupled three-dimensional (3-D) climate-lake-ice modeling system aiming at improving large, deep lakes simulation in regional climate modeling and accurately resolving the hydro-climate interactions. Model results are compared to a wide variety of observational data.  Our results demonstrated the unique skill of the coupled 3-D modeling system in reproducing the regional climate trend and climatic variability and in capturing the physical characteristics of the Great Lakes by fully resolving the hydrodynamics of the large and deep lakes. The simulations of climatology and spatiotemporal variability of lake thermal structure, lake ice are significantly improved with remarkable accuracy. Given the reliable model performance, the surface heat and water budgets are estimated at both seasonal and annual time scales. The implications of the model system and its further development are also discussed.