Hydrologic models provide a comprehensive tool to estimate streamflow response to environmental variables. Yet, an incomplete understanding of physical processes and challenges associated with scaling processes to a river basin, introduces model uncertainty. Here, we apply generalized additive models of location, scale and shape (GAMLSS) to characterize this uncertainty in an Atlantic coastal plain watershed system. Specifically, we describe distributions of residual errors in a two-step procedure that includes model calibration of the soil and water assessment tool (SWAT) using a sequential Bayesian uncertainty algorithm, followed by time-series modeling of residual errors of simulated daily streamflow. SWAT identified dominant hydrological processes, performed best during moderately wet years, and exhibited less skill during times of extreme flow. Application of GAMLSS to model residuals efficiently produced a description of the error distribution parameters (mean, variance, skewness, and kurtosis), differentiating between upstream and downstream outlets of the watershed. Residual error distribution is better described by a non-parametric polynomial loess curve with a smooth transition from a Box–Cox t distribution upstream to a skew t type 3 distribution downstream. Overall, the fitted models show that low flow events more strongly influence the residual probability distribution, and error variance increases with streamflow discharge, indicating correlation and heteroscedasticity of residual errors. These results provide useful insights into the complexity of error behavior and highlight the value of using GAMLSS models to conduct Bayesian inference in the context of a regression model with unknown skewness and/or kurtosis.