Annual maxima of daily precipitation sums can be typically described well with a stationary generalized extreme value (GEV) distribution. In many regions of the world, such a description does also work well for monthly maxima for a given month of the year. However, the description of seasonal and interannual variations requires the use of non-stationary models. Therefore, in this paper we propose a non-stationary modeling strategy applied to long time series from rain gauges in Germany. Seasonal variations in the GEV parameters are modeled with a series of harmonic functions and interannual variations with higher-order orthogonal polynomials. By including interactions between the terms, we allow for the seasonal cycle to change with time. Frequently, the shape parameter ξ of the GEV is estimated as a constant value also in otherwise instationary models. Here, we allow for seasonal–interannual variations and find that this is beneficial. A suitable model for each time series is selected with a stepwise forward regression method using the Bayesian information criterion (BIC). A cross-validated verification with the quantile skill score (QSS) and its decomposition reveals a performance gain of seasonally–interannually varying return levels with respect to a model allowing for seasonal variations only. Some evidence can be found that the impact of climate change on extreme precipitation in Germany can be detected, whereas changes are regionally very different. In general, an increase in return levels is more prevalent than a decrease. The median of the extreme precipitation distribution (2-year return level) generally increases during spring and autumn and is shifted to later times in the year; heavy precipitation (100-year return level) rises mainly in summer and occurs earlier in the year.
