The use of the generalized least squares (GLS) technique for estimation of hydrological regression models has become good practice in hydrology. Through a regression model, a simple link between a particular hydrological variable and a set of catchment descriptors can be established. The regression residuals can be treated as the sum of sampling errors in the hydrological variable and errors in the regression model. This article presents a method for exploring and developing a parameterized form for the cross-correlation between the regression model errors. Given an initial GLS analysis, a reweighted set of regression residuals is defined such that the covariance of these residuals is essentially similar to that of the model errors. The cross-products of the reweighted regression residuals, pooled within bins, are then used to identify a structure and to fit a parameterized form for the cross-correlations of the regression errors. These estimated cross-correlations are then used to inform improved GLS and reweighting steps, leading to a recursive procedure. The main advantage of the recursive GLS procedure is that it allows for a simple demonstration of the actual existence of model error correlation as well as for exploring suitable models for the correlation. However, while the procedure for identifying a parametric description for the model error correlation is reasonable, estimation of the parameters within the recursive procedure is affected by the data-binning step. Thus it is suggested that once a structure for the correlation has been decided, further data analysis—such as final decisions about variables included in the regression model and final estimation of parameters—should be undertaken within a maximum likelihood framework. The procedure has been tested on annual maximum flow data from 602 catchments located throughout the United Kingdom. A set of Monte Carlo experiments further confirmed the ability of the recursive GLS procedure to correctly identify and estimate the true model error correlation.