Looking into the resources essential for socio-economic development, water resources are inevitable. Food, feed, and shelter requirements are growing rapidly with the increase in the population of our world and on the other hand water resources are becoming scarcer to fulfill these demands. Therefore, it is highly required to manage these resources in the best possible way and to contribute towards their sustainability.Water can play a key facilitating role to ensure sustainable development if both groundwater and surface water can be managed efficiently and equitably through the application of fully integrated systems. This study represents an integrated methodology by using the physically-based surface– subsurface flow model (HydroGeoSphere) for the estimation of the effects of climate change and groundwater withdrawal on the corresponding groundwater aquifer systems. This study is conducted for one of the aquifers (Leech lake aquifer) of the City of Yorkton which is recognized as one of the largest urban centers in Saskatchewan depending entirely on groundwater. Finite element based fully coupled model, HydroGeoSphere (HGS) has been selected because it acts as a comprehensive and integrated framework accounting for both surface and subsurface flows, solute and energy transport. In addition, simpler either loosely or externally coupled models do not seem to offer the same level of flexibility in assigning properties and boundary conditions. Digital Elevation Model (DEM) for the subsurface layers has been constructed in Surfer 16 by means of the cross-sectional log data obtained from the City of Yorkton whereas the DEM for the surface layer has been extracted from the Geospatial Canada website. This defines the bottom and the top layers, respectively, and in between lies the computational domain. Validation and calibration exercises were conducted to build confidence in the modeling approach and to estimate some parameters. The validation exercises were conducted against groundwater head observations of four monitoring wells in the period from 2002 to 2015. Simulation results were showing reasonably accurate matches with the observations and this accuracy has been certified by some of the frequently used statistical error analysis tools. From both observations and simulations, it has been found that the drawdowns in the existing observation wells generally lie within the range of 1 to 3 m whereas it could go a maximum of 5 m in the location of production wells. Different recharge and withdrawal scenarios have been evaluated for providing a sustainable management plan, starting from 1 year (short-term) to 100 years (long-term). In addition, an evacuation scenario has been taken into consideration under adverse climate change conditions (assuming neither any rainfall nor snowmelt). This specific case is being simulated to know the withdrawal rate of pumping wells which could evacuate the whole reservoir. Determination and interpretation of exchange flow rates at the interface between surface water bodies and its associated subsurface aquifer layers is another significant outcome of this research. That is to say, this fully coupled model advances the illustration of mutually dependent processes like recharge, which is vital in the context of climate change. This modeling case study has been executed to support the viability of this integrated approach in replicating the existing condition which could be further used as a prediction model in the light of quick and unpredictable changes.