Highway networks represent a type of linear land use crossing streams and sensitive water bodies. Stormwater runoff generated from highway surface has long received much attention due to the presence of a variety of contaminants, such as sediments, heavy metals, hydrocarbons, and nutrients. Mandated by the Clean Water Act and U.S. EPA's regulations, state transportation agencies are required to secure their National Pollution Discharge and Elimination System stormwater permits; these permits issued to the agencies shall also incorporate the implementation requirements of established TMDLs. Compliance with the increasing TMDLs has brought a great pressure to the agencies in stormwater management and becomes an emerging issue requiring technically sound watershed-scaled TMDL modeling methodologies to estimate pollutant loads from the highway right-of-way land use. To meet this demand, this research adopts a new approach to develop the watershed-scaled probabilistic volume-to-breakthrough water quality model, i.e., PVbtWQM. It consists of three components: a) a hydrologic connectivity evaluation component; b) a TN (i.e., NO3-N and TKN) loading and TMDL assessment component; and c) a TP loading and TMDL assessment component. Prior to model construction, a comprehensive review is provided to include existing methodologies in modeling highway stormwater runoff quantity and quality, and a watershed TMDL modeling case by using WARMF in the context of highways, Falls Lake Watershed Analysis Risk Management Framework Development. Then, the algorithms of PVbtWQM are derived to evaluate the hydrologic connectivity of highways to receiving streams and quantify highway stormwater runoff pollutant loadings reaching streams as well as their associated uncertainties. They are further programmed in spreadsheets of Microsoft Excel and an 8-step procedure of model manipulation is presented. In addition, an approach to performing a comprehensive terrain analysis with known streams and lakes on a LiDAR-based DEM are also provided to extract the information of highway runoff drainage systems and measure the lengths of diffuse flow pathways by using the Arc Hydro tools and ArcMap. PVbtWQM has been applied to assessing stormwater runoff nutrient TMDLs from the NCDOT road land use for the 26.74-mi2 (17,114-acre) Lake Orange watershed in Orange County of North Carolina. It yields 0.413 ± 0.407 kg/d (1.426 ± 1.408 lb/ac.yr) of TN and 0.090 ± 0.121 kg/d (0.312 ± 0.417 lb/ac.yr) of TP. Without considering the nutrient reduction through diffuse flow paths, it yields 0.918 ± 0.715 kg/d (3.172 ± 2.473 lb/ac.yr) of TN and 0.182 ± 0.208 kg/d (0.628 ± 0.719 lb/ac.yr) of TP source loadings. For the same NCDOT land use, WARMF yields 0.139 kg/d (0.480 lb/ac.yr) of TN and 0.021 kg/d (0.072 lb/ac.yr) of TP source loading. Comparing nutrient loadings of these two models to those estimated by the simple method, PVbtWQM's simulation results fall in the range between the minimal EMCs based and the mean EMCs based estimates of the simple method, lower than and approaching the upper quartile of that range for both TN and TP loadings. If considering the reduction via diffuse flow pathways, PVbtWQM's results are 21% and 11% higher than the minimal EMCs based estimates of the simple method for TN and TP loadings, respectively. WARMF's results are significantly lower than the minimal EMCs based estimates of the simple method by 59% and 79% for TN and TP, respectively, and also lower than the national low values and NC secondary road runoff nutrient loads. This study has shown that PVbtWQM provides more acceptable simulation results and is more reliable than WARMF for estimating nutrient loading from the same highway land uses, and that both TN and TP loadings of stormwater runoff from NCDOT highways have been underestimated in the Falls Lake WARMF Development where the corresponding system coefficients should be further calibrated. Meanwhile, the simulation results of PVbtWQM have shown that the road-to-stream hydrologic connectivity will increase as precipitation increases, and that different types of flow pathways have significant impact on highway stormwater runoff nutrient loadings. In the Lake Orange watershed, the overall hydrologic connectivity of the highway network to receiving streams is 0.32 ± 0.14, ranging from 0.23 to 0.86 during the simulation period between 2004 and 2007; the original road runoff TN and TP loadings are reduced through diffuse flow pathways by 55% and 51%, respectively. Also, it has been shown that the propagation uncertainties in the final results for both TN and TP are quite large, ranging from 78% to 134%, and would be a big concern in highway stormwater TMDL assessment. In short, PVbtWQM has been proven as an informative watershed-scaled highway stormwater TMDL modeling methodology. It provides a new option for state transportation agencies to assess highway stormwater pollutant loads and support their stormwater management decision making.