Probable maximum flood (PMF) has been used in large hydraulic infrastructure design for decades. However, a complete framework for deriving PMF with uncertainty analysis is lacking. This study investigated from probable maximum precipitation (PMP) with uncertainty to probable maximum storm (PMS), rainfall-runoff (R-R) model with uncertainty, and PMF with uncertainty. A basin-scale model, based on the Hershfield method, was developed for calculating PMP and uncertainty. It showed that PMP and uncertainty from the proposed method were more reliable than other statistical methods because of less restrictive assumptions and the improvement in methodology. The average improvements in terms of the difference percentages for 1-hour, 6-hour and 24hour duration PMPs were 53.84%, 81.04% and 72.60%, respectively. The improvements in uncertainty by the Delta method were 15.54%, 9.71%, 8.93%, respectively. Hydrometeorological Report No. 52 (HMR 52) with updated data was adopted to design the PMS from PMP. The updated storm records indicated a different within/without-storm depth-area relation (WWSDA) than HMR 52. The PMP values with 95% confidence interval (CI) were transferred to PMSs utilizing the DAD maps and the WWSDA relations. With the GLUE methodology, the impacts of different sources of uncertainty on the R-R model performance and uncertainty were investigated. Results showed that the uncertainties were more sensitive to smaller basin area, average area of subbasin, precipitation intensity, and observed peak flow, which means that for extreme events modeling with opposite basin features, the selection of model input and model setup had iii less impact on model performance and uncertainty. PMSs were then routed by the R-R model with uncertainty to obtain the PMF with uncertainty. Results showed that in lower Brazos River basin, the PMF peak flow derived by the average PMP value was 12716 CMS, which was around three times the peak flow from hurricane Harvey (3764 CMS).