Watershed-scale coupled surface water (SW) – groundwater (GW) flow modeling was used to examine changes in streamflow and SW – GW interaction resulting from irrigation and associated SW diversions and GW pumping. The U.S. Geological Survey (USGS) model GSFLOW, an integration of the USGS Precipitation-Runoff Modeling System (PRMS) and the Modular Ground-Water Flow Model (MODFLOW), was utilized for this effort. Processes represented in the model include daily rain, snowfall, snowmelt, streamflow, surface runoff, interflow, infiltration, soil-zone evapotranspiration (ET), and subsurface unsaturated and saturated GW flow and ET. The upper Smith River watershed, an important agricultural and recreational fishing area in west-central Montana, was used to develop the model framework including watershed climate, topography, hydrography, vegetation, and soil properties as well as the scenario used to represent current irrigation practices. The 640 square kilometer modeled watershed area ranged in elevation from 1476 m to 2570 m and was discretized into coincident 200 m by 200 m hydrologic response units (for climate and soil zone flow processes) and grid blocks (for unsaturated zone and GW flow processes) resulting in a grid of 180 rows by 160 columns. The subsurface GW system was discretized into 6 layers (two 50-m thick layers overlying five 100-m thick layers) representing Quaternary alluvium, Tertiary sedimentary rock and bedrock with maximum thicknesses of 50 m, 200 m and 500 m, respectively. The time period from 10/1/2004 to 10/1/2010 was simulated using 2192 daily stress periods. Observed daily temperature maximum, minimum and precipitation records at the White Sulphur Springs 2 weather station were used with temperature and precipitation lapse rates to distribute elevation-dependent temperature and precipitation across the watershed. Initial estimates of parameters controlling soil zone storage, infiltration from the soil zone to GW, and interflow were related to watershed properties such as land-surface slope and soil properties. Potential ET was calculated in GSFLOW using the Jensen-Haise approach. The model ET extinction depth was based on the watershed vegetation distribution, assuming alfalfa was planted in sprinkler irrigated areas, and literature plant root depths. Subsurface saturated and unsaturated zone properties were related to the hydrogeologic units (i.e. alluvium, Tertiary sedimentary rock and bedrock). The stream channel network was constructed by using the Arc Hydro flow accumulation tool and selecting the network of stream segments with flow accumulating from at least 20 km2 (500 grid blocks). Each stream segment of the model was assigned a 12-m wide cross-section that allowed stream width and depth to vary with flow. Model parameters were adjusted to reproduce the general observed streamflow patterns and GW level distributions in the watershed. Model results were used to compare streamflow, GW recharge and SW-GW exchange in the watershed under natural, pre-irrigation conditions (PreIrr scenario); for current irrigation practices involving mainly SW diversion with flood and sprinkler irrigation (IrrCurrent scenario); and for an irrigation scenario based solely on GW pumping for sprinkler irrigation (IrrGW scenario).