Stream flow from snow-dominated mountain ecosystems is an important source of water in the Western US and other parts of the world. Mountain watersheds are complex systems with high spatial variability in elevation, slope, aspect, vegetation, and subsurface properties. Integrated hydrological models such as GEOtop have been designed to simulate water and energy fluxes in these complex environments. GEOtop couples a canopy and surface energy balance to a fully distributed hydrological model that additionally simulates root-water uptake, snow evolution, and soil freezing. The model is used in this study to test an optimization method that efficiently calibrates subsurface properties and accounts for the presence of a significant fraction of stone/rock. The optimization was found to be useful for both point-scale simulations where soil water content was predicted, and catchment-scale simulations where stream flow was the main model output. The optimization method was applied at the Dry Creek Experimental Watershed in the Boise Front, Idaho, and the Libby Creek Experimental Watershed in the Medicine Bow Mountains, Wyoming. In a separate study in the Medicine Bow Mountains, GEOtop predicted above and below-canopy turbulent fluxes were compared against eddy covariance data. The comparison showed significant discrepancies, especially for below-canopy turbulent flux.