In irrigated semi-arid watersheds, over-fertilization often leads to nitrogen (N) and phosphorus (P) contamination in aquifers and river systems. Modelling tools often are used to evaluate the effect of management practices on nutrient contamination levels. In this study, we assess a suite of best management practices (BMPs) for N and P in a regional irrigated stream-aquifer system using a new numerical model. As it is necessary to consider the interaction of flow in the stream and aquifer, SWAT-MODFLOW model that couples surface and groundwater flow was used. The reactive processes were modelled using RT3D that simulate the reactive groundwater transport. The model is based on the coupled flow model SWAT-MODFLOW, with the groundwater reactive transport model RT3D included to simulate the reactive groundwater transport of NO3 and soluble P and their interactions within the soil-aquifer-stream system. The assessment is performed for a highly managed 732 km2 region in the Lower Arkansas River Valley (south-eastern Colorado), with model results compared to observed groundwater nutrient concentration, in-stream nutrient concentration and loading, and crop yield. A total of 28 BMPs scenarios are evaluated regarding their impact on NO3 and soluble P contamination in the aquifer and river system. The most effective individual BMP in most areas is to decrease fertilizer by 30%, resulting in NO3 and soluble P reduced by 20% and 2% for groundwater concentrations, 25% and 10% for river concentrations, and 27% and 6% for mass loadings into surface water, respectively. Combinations of using 30% irrigation reduction, 30% fertilization reduction, 60% canal seepage reduction, and conservation tillage yield the greatest overall impact to lower NO3 and soluble P concentrations up to 41% and 8% in groundwater, 52% and 40% in river, and 63% and 49% for mass loadings. Targeting BMPs on localized problem areas shows great promise in reducing contamination while maintaining region-wide crop yield. The study demonstrates the SWAT-MODFLOW-RT3D modelling code is a useful tool to examine NO3 and P transport and quantify BMP effects in groundwater-driven watersheds.