Climate change affects all the components of the hydrological cycle. Starting from precipitation distribution, climate alterations have direct effects on both surface water and groundwater in terms of their quantity and quality. These effects lead to modifications in water availability for agriculture, ecology and other social uses. Change in rainfall patterns also affects the runoff of natural rivers. For this reason, studying runoff data according to classical hydrological approaches, i.e., statistical inference methods that exploit stationary probability distributions, might result in missing important information relevant to climate change. From this point of view, a new approach has to be found in the study of this type of data that allows for non-stationary analysis. In this study, the statistical framework known as Generalized Additive Models for Location, Scale and Shape (GAMLSS), which can be used to carry out non-stationary statistical analyses, was applied in a non-stationary frequency analysis of runoff data collected by four gauges widely distributed across Sicily (Italy) in the period 1916–1998. A classical stationary frequency analysis of these runoff data was followed by a different non-stationary frequency analysis; while the first was made using annual rainfall as a covariate, with the aim of understanding how certain statistical parameters of runoff distribution vary with changes in rainfall, the second derived information about the temporal variability of runoff frequencies by considering time as a covariate. A comparison between stationary and non-stationary approaches was carried out using the Akaike information criterion as a performance metric. After analyzing four different probability distributions, the non-stationary model with annual rainfall as a covariate was found to be the best among all those examined, and the three-parameter lognormal the most frequently preferred distribution.