The application of different multivariate statistical approaches for the interpretation of a large and complex data matrix obtained during a monitoring program of surface waters in Northern Greece is presented in this study. The dataset consists of analytical results from a 3-yr survey conducted in the major river systems (Aliakmon, Axios, Gallikos, Loudias and Strymon) as well as streams, tributaries and ditches. Twenty-seven parameters have been monitored on 25 key sampling sites on monthly basis (total of 22,350 observations). The dataset was treated using cluster analysis (CA), principal component analysis and multiple regression analysis on principal components. CA showed four different groups of similarity between the sampling sites reflecting the different physicochemical characteristics and pollution levels of the studied water systems. Six latent factors were identified as responsible for the data structure explaining 90% of the total variance of the dataset and are conditionally named organic, nutrient, physicochemical, weathering, soil-leaching and toxic-anthropogenic factors. A multivariate receptor model was also applied for source apportionment estimating the contribution of identified sources to the concentration of the physicochemical parameters. This study presents the necessity and usefulness of multivariate statistical assessment of large and complex databases in order to get better information about the quality of surface water, the design of sampling and analytical protocols and the effective pollution control/management of the surface waters.