This model archive contains model files and a complete repeatable workflow for source water delineation for groundwater supply wells in the Neversink/Rondout Basins of New York. The U.S. Geological Survey, in cooperation with the NYSDEC and NYSDOH, began an investigation in 2019 with the general objectives of (1) improving understanding of the regional groundwater flow system (2) delineating areas contributing recharge to eight priority production wells, and (3) quantifying the uncertainty of these contributing areas in a probabilistic way that can be used to inform decision-making related to priority well source water protection. To accomplish these objectives, a MODFLOW 6 groundwater model (version 6.2.1) was created encompassing eight priority wells and the surrounding flow system, which includes portions of the Neversink and Rondout basins, in Sullivan County and Ulster County, NY. The model was built using python tools (such as flopy, modflow-setup, and sfrmaker) that facilitate transparent and repeatable model development using existing datasets. The model parameters were estimated with a stepwise approach using an iterative ensemble smoother (iES) implementation of the Parameter ESTimation (PEST) software PEST++ (version 5.0.10). We evaluated initial “best guess” parameter bounds with a prior Monte Carlo analysis. Results of the first prior Monte Carlo analysis were used to make informed adjustments to model parameter bounds (typically resulting in expanded bounds), and a second prior Monte Carlo analysis was run to identify improved ranges for model parameters during history matching. The history-matching effort produced an ensemble of parameter values for the groundwater-flow model that spans the range of values within prior uncertainty bounds and is informed by the historical observation data, within a reasonable range of uncertainty on those observations. This history-matched ensemble was used to delineate the areas contributing recharge to priority wells by running the groundwater-flow and particle tracking (MODPATH7) models for each ensemble member. Deterministic contributing areas computed for each ensemble member were aggregated to produce maps showing the probability that a location contributes recharge to priority wells. Finally, this process was repeated for six scenarios, representing a wide range of possible pumping levels, to incorporate uncertainty in pumping rates related to population growth or other future decisions. Increasing pumping rates generally led to larger contributing recharge areas and larger areas of high probability. These maps offer a view of the overall uncertainty of the areas contributing recharge to priority wells in the study area and provide a tool for risk-based decision making related to the source water protection of these wells.