Area-covering postprocessing of ensemble precipitation forecasts using topographical and seasonal conditions

This study focusses on postprocessing of ensemble precipitation predictions from the NWP system COSMO-E (COnsortium of Small-scale MOdelling). At the time of writing, this is the operational probabilistic high-resolution NWP system of MeteoSwiss, the Swiss national weather service. COSMO-E is run at 2x2km resolution for an area centered on the Alps extending from northwest of Paris to the middle of the Adriatic. An ensemble of 21 members is integrated twice daily at $00:00$ and $12:00$ UTC for five days (120 hours) into the future. Boundary conditions for the COSMO-E forecasts are derived from the operational global ensemble prediction system of the European Centre for Medium-range Weather forecasting (ECMWF).

We focus on daily precipitation amounts in Switzerland. As suggested by Messner (2018), the ensemble predictions and the observations of the daily precipitation amount are square-root transformed before the postprocessing. We use observed daily precipitation at MeteoSwiss weather stations to calibrate the ensemble forecasts. This paper relies on two observation datasets, one for the elaboration of the methodology and one for the subsequent assessment of the models of choice. The datasets are presented in Table 1. The first dataset (subsequently referred to as Dataset 1) consists of ensemble forecasts and verifying observations for the daily precipitation amount between January $2016$ and July $2018$ whereby a day starts and ends at $00:00$ UTC. The data is available for $140$ automatic weather stations in Switzerland recording sub-daily precipitation. The second dataset provides the observed daily precipitation amounts of $327$ additional stations (on top of the $140$ ones of Dataset 1) that record only daily sums. For historical reasons, these daily sums are aggregated from $06:00$ UTC to $06:00$ UTC of the following day. For the purpose of uniformity, these daily limits are adopted for all stations in Dataset 2. The larger number of stations from Dataset 2 is only available from June $2016$ to July $2019$. All stations from Dataset 1 are also part of Dataset 2. The stations of both datasets are depicted in Figure 2.

Since COSMO-E makes forecasts for five days into the future, the different intervals between the forecast initialization time and the time for which the forecast is valid have to be taken into account. This is referred to as forecast lead time. The possible lead times for a prediction initialized at $00:00$ UTC are $1$, $2$, $3$, $4$, $5$ days and $1.5$, $2.5$, $3.5$ and $4.5$ days for one initialized at $12:00$ UTC respectively. Figure 3 illustrates the possible lead times for both datasets. For Dataset 2 the forecast lead times increase from 24 hours to 30 hours for the first day due to the different time bounds of aggregation. Also, only four complete forecasts can be derived from each COSMO-E forecast with Dataset 2. We use Dataset 1 reduced to forecast-observation pairs with lead time equals $3$ for the elaboration of the methodology. This selection procedure depends strongly on the dataset bearing the danger of overfitting. To assess this risk, the models of choice will be evaluated with Dataset 2. This is done for all lead times between $1$ and $4$.

In addition to the forecast-observation pairs and the station location (latitude, longitude, and altitude), topographical indices derived from a digital elevation model with 25m horizontal resolution are used. The topographical indices are available on 7 grids with decreasing horizontal resolution from $25$m to $31$km describing the topography from the immediate neighbourhood (at $25$m resolution) to large-scale conditions (at $31$km). The topographical indices include the height above sea level (DEM), a variable denoting if the site is rather in a sink or on a hill (TPI), variables describing aspect and slope of the site and variables describing the derivative of the site in different directions.

In this paper, the observed precipitation amount (at any specified location and time) is denoted as $y$. $y$ is seen as a realization of a non-negative valued random variable $Y$. The $K$ ensemble members are summarized as $\boldsymbol{x}=(x_{1},...,x_{K})$. Predictive distributions for $Y$ are denoted by $F$ and stand either for cumulative distribution functions (CDFs) or probability density functions (PDFs). In literature, a capital $F$ is used to refer indistinguishably to the probability measure or its associated CDF, a loose convention that we follow for simplicity within this paper. A forecast-observation pair written as $(\boldsymbol{x}_{i},y_{i})$ refers to a raw ensemble forecast and the corresponding observation. A pair written as $(F_{i},y_{i})$ generally indicates here, on the other hand, that the forecast is a postprocessed ensemble prediction.

To assess a postprocessing model, the conformity of its forecasts and the associated observations is rated. In the case of probabilistic forecasts, a predictive distribution has to be compared with single observation points. We follow Gneiting, Balabdoui and Raftery (2007) by aiming for a predictive distribution maximizing the sharpness subject to calibration.

In this section, we present a conventional ensemble postprocessing approach as a starting point for later extensions. We compared several EMOS models and Bayesian Model Averaging in a case study with Dataset 1 whereby a censored inhomogeneous Logistic regression (cNLR; Messner et al. 2014) model turned out to be the most promising approach. The performance of the models has been assessed by cross-validation over the months of Dataset 1. Of the $31$ months available, the months are removed one at a time from the training data set and the model is trained with the remaining $30$ months. Predictive performance is then assessed based on comparison between the observations at each left-out month from the models trained on the remaining months. The basic models are tested in two versions: For the global version, the model is trained with the data of all stations allowing the later application of the model to all stations simultaneously. The local version requires that models are trained individually for each station with the past data pairs of this station. More details on the alternative approaches to the cNLR model and the results from the comparison of approaches can be found in the supplementary material.

The cNLR approach is a distributional regression model: We assume that the random variable $Y$ describing the observed precipitation amount follows a probability distribution whose moments depend on the ensemble forecast. To choose a suitable distribution for $Y$, we take into account that the amount of precipitation is a non-negative quantity that takes any positive real value (if it rains) or the value zero (if it does not rain). These properties are accounted for by appealing to a zero censored distribution. We assume that there is a latent random variable $Y^{*}$ satisfying the following condition (Messner, Mayr and Zeileis 2016):

In this way, the probability of the unobservable random variable $Y^{*}$ being smaller or equal than zero is equal to the probability of $Y$ being exactly zero.

For the choice of the distribution of $Y^{*}$, we have compared different parametric distributions: a Logistic, Gaussian, Student, Generalized Extreme Value and a Shifted Gamma distribution. For the Logistic distribution, which has achieved the best results, $Y^{*}\sim\mathcal{L}(m,s)$ with location $m$ and scale $s$ has probability density function

The expected value and the variance of $Y^{*}$ are given by:

The location $m$ and the scale $s$ of the distribution have to be estimated with the ensemble members. We note that in the COSMO-E ensemble, the first member $x_{1}$ is obtained with the best estimate of the initial conditions whereas the other members are initialized with randomly perturbed initial conditions. The first member is thus not exchangeable with the other members and it is therefore reasonable to consider this member separately. Members $x_{2},x_{3},...,x_{21}$ are exchangeable, meaning that they have no distinct statistical characteristics. Within the location estimation, they can therefore be summarized by the ensemble mean without losing information (Wilks 2018). Taking this into account, we model the location $m$ and the scale $s$ of the censored Logistic distribution as follows:

where $\overline{x}$ is the ensemble mean and $SD(\boldsymbol{x})$ is the ensemble standard deviation. The five regression coefficients are summarized as $\psi=(\beta_{0},\beta_{1},\beta_{2},\gamma_{0}$, $\gamma_{1})$.

A case study with Dataset 1 showed that the performance of the local cNLR model is better than the one of the global cNLR model (see supplementary material). Similar results have been presented for example by Thorarinsdottir and Gneiting (2010) in a study about wind speed predictions. The local model, however, cannot be used for area-covering postprocessing. To improve the performance of the global model, it is enhanced with topographical covariates. The idea is to fit the regression coefficients only with training data from weather stations which are topographically similar to the prediction site.

We assume that the training dataset consists of $n$ forecast-observation pairs $(\boldsymbol{F_{i}},y_{i})$ with $i=1,2,...,n$. The global cNLR model estimates $\psi$ by minimizing the mean CRPS value of all these training data pairs. To select only or favour the training pairs from similar locations, we use a weighted version of the mean CRPS value¹¹1Literature knows another kind of weighted CRPS value: Threshold and quantile weighted versions of the CRPS are used when wishing to emphasize certain parts of the range of the predicted variable. The threshold weighted version of the CRPS is given by $$CRPS_{u}(F,y)=\int_{-\infty}^{\infty}\left(F(x)-\mathbb{I}_{[y,\infty]}(x)\right)^{2}u(x)\enskip dx,$$ (15) where $u$ is a non-negative weight function on the real line (Gneiting and Ranjan 2011, consult their paper for the analogous definition of the quantile weighted version of the CRPS). as the cost function $c$, which is minimized for the fitting of the coefficient vector $\psi$:

$CRPS(F_{i}^{\psi},y_{i})$ refers to the CRPS value of data pair $(F_{i}^{\psi},y_{i})$ where the predictive distribution $F_{i}^{\psi}$ depends on the coefficient vector $\psi$. We use $(F_{i},y_{i})$ as shorthand for $(F_{i}^{\psi},y_{i})$. The weight $w_{i}(s)$ of training data pair $(F_{i},y_{i})$ depends on the similarity between the location it originated and the location $s$ where we want to predict. We set $w_{i}(s)=1$ if training data pair $i$ originated in one of the $L$ closest (be it with respect to the euclidean distance or to some other dissimilarity measure) stations to the prediction site $s$. For the training pairs from the remaining stations, we set $w_{i}(s)=0$. This ensures that the training dataset is restricted to the data from the $L$ stations which are most similar to the prediction site $s$. Consequently, the coefficient vector $\psi^{*}(s)$ which minimizes $c(\psi;s)$ depends on the location $s$.

Following Lerch and Baran (2018), the similarity between locations is quantified with a distance function, which, in our case, is intended to reflect the topography of the respective locations. From the topographical dataset we have about $30$ variables in $8$ resolutions at our disposal. To get an insight regarding which ones to use, we examine the topographical structure of the raw ensemble’s prediction error. We compare the observed daily precipitation amounts with the ensemble means and take the station-wise averages of these differences. These preliminary analyses were made with the first year ($2016$) of Dataset 1. The mean prediction errors per station are depicted in Figure 4(a). The topographical structure of the ensemble prediction error seems to be linked to the mountainous relief of Switzerland.

In a first approach, we define the similarity of two locations via the similarity in their distances to the Alps. It turns out that such an approach depends on numerous parameters (we refer the reader to the supplementary material for more details). The proposed alternative is to focus on the variable DEM describing the height above sea level and use the values provided by a grid with low resolution (here 31km horizontal grid spacing). This ensures, that the small-scale fluctuations in the height are ignored such that the large-scale relief is represented. Figure 4(b) shows the same station-wise means as Figure 4(a), but this time the values of each station are plotted versus their height above the sea level (DEM) in resolution $31km$. The best fit with a polynomial is achieved by modelling the ensemble mean bias as a linear function of the DEM variable; the solid line is depicting this linear function.

The linear dependency appearing in Figure 4(b) motivates the following choice of a distance function to measure the similarity between two locations. Let us define a function $DEM$ which maps a location $s$ to its height above the sea level in the resolution $31$km:

where $\mathcal{D}\subseteq\mathbb{R}^{2}$ is a set with all the coordinate pairs lying in Switzerland. The similarity of locations $s_{1}$ and $s_{2}$ is then measured by the following distance:

Based on this distance, we determine the $L$ stations of the training data set which are most similar to the prediction location, i.e. the stations which have the smallest distances $d_{DEM}$. Let

be the ordered distances $d_{DEM}(s,s_{j})$ between the $m$ stations from the training data and the prediction location $s$. Let $s_{i}$ be the location of the station where forecast-observation pair $(F_{i},y_{i})$ originated. Then, the weights with respect to the distance $d_{DEM}$ are defined as:

These weights ensure that the data pairs, which originated at one of the $L$ most similar stations, get weight $1$ and the remaining get weight $0$.

Besides this approach, several other topographical extensions of the cNLR model have been tested (with Dataset 1): For their spatial modelling, Khedhaouiria, Mailhot and Favre (2019) propose to vary the postprocessing parameters by expressing them as a function of spatial covariates. We have applied a similar approach and integrated the topographical covariates in the location estimation of the cNLR model. To reduce the number of predictor variables, the topographical variables have been summarized by Principal Components. Additionally, we used the glinternet algorithm of Lim and Hastie (2013) to uncover important additive factors and interactions in the set of topographical variables. A more basic weighted approach has been based on Euclidean distances in the ambient (two and three dimensional) space. All extensions of the cNLR model have been compared with the local and the global fit of this very model.

As the target of this work is to develop an area-covering postprocessing method, the extended and the global models are trained with data where the predicted month and the predicted station are left out. This simulates the situation where postprocessing must be done outside a weather station, i.e. without past local measurements. The training period is set to the last year ($12$ months) before the prediction month. Consequently, forecasting performance can only be assessed with (test) data from $2017$ and $2018$. The case study with Dataset 1 showed that all other topographical approaches are less successful than the DEM approach, more details and the results can be found in the supplementary material (in the section about extension approaches).

In addition to our efforts to quantify similarities between locations, we also aim to investigate ways of further improving postprocessing outside of measurement networks by accounting for seasonal specificities. To examine the seasonal behavior of the local and the global cNLR model, we focus on their monthly mean CRPS values and compare them with the ones of the raw ensemble. Figure 5 shows the monthly skill of the global and the local cNLR model. We use the mean over all the stations from Dataset 1 and depict the months between January $2017$ and July $2018$ such that we have $12$ months of training data in each case. A positive skill of $10\%$ means for example that the mean CRPS value of the raw ensemble is reduced by $10\%$ through postprocessing, a negative skill indicates analogously that the mean CRPS value is higher after the postprocessing. The global model has negative skill in February and March $2017$ and between June and October $2017$. The values are especially low in July and August $2017$. The local model has a positive skill in most months, the postprocessing with this model decreases the forecasting performance only between June and September $2017$. We use these results as a first indication that the postprocessing of ensemble precipitation forecasts is particularly challenging in summer and early autumn.

Next, we are interested in whether there are regional differences in the model performance within a month. The global cLNR model is used as we will extend this model afterwards. We plot the maps with the station-wise means of the skill exemplary for the month with the best skill (February $2018$) and the one with the worst (July $2017$). The maps depicted in Figure 6 show that the skill of the global cNLR model varies between different weather stations. Again, the structure seems to be related to the mountainous relief of Switzerland. We note that for both months the skill in the Alpine region is distinctly higher than in the flat regions.

We use this knowledge to develop an approach which tries firstly to clarify whether the postprocessing in a given month at a given prediction location is worthwhile. The idea is to “pretest” the model with data of similar stations and from similar months by comparing its performance with that of the raw ensemble. For this purpose, the year of training data is first reduced to the data pairs from topographically similar stations, whereby the similarity is measured with the distance $d_{DEM}$ defined in Equation (18). Afterwards, this training dataset is split into two parts: Traintrain and Traintest. The model is adapted a first time with the dataset Traintrain. Afterwards, the performance of this model is assessed with the second part (Traintest) by comparing the mean CRPS of the model with the mean CRPS of the raw ensemble.

The months of the Traintest dataset are selected such that they are seasonally similar to the prediction month. To split the training dataset, three approaches are compared:

• •

  Pretest 1: Pretest with the same month as the prediction month from the year before (Example: January $2017$ for January $2018$)

• •

  Pretest 2: Pretest with the month before the prediction month (Example: December $2017$ for January $2018$)

• •

  Pretest 3: Pretest with both of these months (Example: January $2017$ and December $2017$ for January $2018$).

Let us define the set of indices of training data pairs out of Traintest:

with cardinality $H$. Let further $(\boldsymbol{x_{i}},y_{i})$ be a forecast- observation pair where the forecast is the raw ensemble. $(F_{i},y_{i})$ is a pair with a postprocessed forecast. If

then the pretesting algorithm decides that the raw ensemble is not postprocessed in the given month at the given location. On the contrary if

then the pretesting algorithm decides that the raw ensemble is postprocessed in the given month at the given location, the fit is done a second time with the whole year of training data.

The Pretest approach has been compared with several other seasonal approaches: In a basic approach, we reduce the training period to months from the same season as the prediction month. Another approach uses the sine-transformed prediction month as an additional predictor variable to model the yearly periodicity, an approach comparable to the one of Khedhaouiria, Mailhot and Favre (2019). The third approach reduces the training data to pairs which have a similar prediction situation (quantified with the ensemble mean and the ensemble standard deviation). The methodology for the comparison has been the same as for the topographical extensions introduced in Section 3.2. The Pretest approach turns out to be the most promising method for our case study with Dataset 1, more details and the comparison results can be found in the supplementary material.

For the subsequent evaluation of postprocessing models with Dataset 2, we select a few postprocessing approaches to document the impact of increasing complexity on forecast quality. We will use the raw ensemble and the local such as the global version of the cNLR model as baselines. Further on, we will evaluate the cNLR model extended by the DEM similarity (cNLR DEM). Finally, we will test this same model extended a second time with the pretest approach (cNLR DEM+PT). For the last two models, we have to fix the amount $L$ of similar stations we use for the topographical extension. For the last model we need to fix additionally the pretesting split. To determine the amount of similar stations in use, the numbers which are multiples of ten between $10$ and $80$ have been tested (compare Figure 7). We use the data from August $2017$ to July $2018$ (seasonally balanced) from Dataset 1 and choose the number resulting in the lowest mean CRPS. For the cNLR DEM model (no Pretest) we determine $L=40$. For the cNLR DEM+PT model, we combine the different pretesting splits with the same numbers for $L$. The cNLR DEM+PT model with the lowest mean CRPS value uses Pretest 3 and $L=40$.

We are interested in the area-covering performance of the models. Therefore, we are particularly interested in the performance at locations which cannot be used in model training (as no past data is available). This is the reason why we assess the models only with the $327$ additional stations of the second dataset. None of these stations have been used during the model elaboration in Section 3. When determining $L$ in chapter 3.4 we used a training dataset with $139$ stations ($139$ instead of $140$ as we trained without the past data from the prediction station). For this reason we carry on with using only the $140$ stations of the first dataset to train the models. This rather conservative approach could be opposed by a Cross Validation over all $467$ stations, for which, however, another choice of $L$ would probably be ideal.

The local version of the cNLR model is not able to perform area-covering postprocessing and needs the additional stations in the training from Dataset 2. Despite this, it is fitted and assessed as a benchmark here. We train the models for each of the $327$ stations and each month between June $2017$ and May $2019$. This ensures that we have one year of training data for each month and that we have seasonally balanced results. An individual fitting per station is necessary as the selection of the most similar stations used in the DEM approaches depends on the station topography. The model must also be adapted monthly, as the pretesting procedure (and the training period) depend on the prediction month.

During the postprocessing, we used consistently the square root of the precipitation amount. The CRPS value, which is in the same unit as the observation, refers to this transformation as well. To get an idea of the actual order of magnitude, the values are converted into the original size, in which the precipitation amount is measured in mm. As a first step, $21$ forecasts are drawn from the fitted censored Logistic model. Afterwards, these values and the corresponding observations are squared and the mean CRPS is calculated as for the raw ensemble. The Brier Score, which assesses the ability of the forecaster to predict if a given precipitation accumulation amount is exceeded, is also evaluated for the squared sample of the predictive distribution. The thresholds used within the Brier Score focus on three precipitation accumulations: No rain ($<$ 0.1 mm/d) , moderate rain ($>$ 5 mm/d) and heavy rain ($>$ 20 mm/d).

First of all, let us give an overview of the different models. Figure 8 depicts the mean CRPS values for the different postprocessing approaches and lead times. We refer to Chapter 3.4 for a recap of the model adjustments concerning the DEM and DEM + Pretest model. For lead time $1$, the global cNLR model is able to reduce the mean CRPS value by $2.3\%$. A further improvement is achieved by the local and the DEM model, which show equivalent performances and reduce the mean CRPS value by $4.5\%$ compared to the raw ensemble. Even slightly better results are delivered by the DEM + Pretest model which reduces the mean CRPS value by $4.8\%$. The skill of the global model decreases with increasing lead time. While the skill is still positive for lead times $1$, $1.5$ and $2$, the model performs roughly equally as the raw ensemble for lead time $2.5$. From lead time $3$, the mean CRPS value of the global model is even higher than the one of the raw ensemble. The local and the DEM model perform about the same for lead times between $1$ and $2.5$, for lead times above $3$ the DEM model performs slightly better. The DEM+Pretest model performs best for all lead times. It reduces the mean CRPS value between $4.8\%$ for lead time $1$ and $2.0\%$ for lead time $4$. It is noticeable that the DEM + Pretest model achieves a near constant improvement in the mean CRPS of approx. $0.07$ for all lead times. Relatively - i.e. as a Skill Score - this corresponds to less and less of the total forecast error. We note additionally, that the improvement which is achieved through the extension of the DEM model with the Pretest depends on the lead time. While the Pretest reduces the mean CRPS of the DEM model only by $0.4\%$ for lead time $1$, the obtained reduction corresponds to a proportion of $1.7\%$ for lead time $4$.

Figure 8 summarizes the average performance of the models over all months. To assess the seasonal performance of the different approaches, the monthly means of the Skill Score are plotted in Figure 9. We use the raw ensemble forecast as reference and depict the results for lead time $1$ and lead time $4$. For lead time $1$, we note that the DEM + Pretest model is the only one with non-negative skill in almost all months, implying that this model only rarely degrades the quality of the ensemble prediction. While the model delivers in summer and early autumn equivalent results as the raw ensemble, the monthly mean CRPS value can be reduced by up to $20\%$ in winter months. The same improvement is achieved with the models without Pretest, but they have a slightly worse overall performance as they degrade forecast performance during summer and autumn. For longer lead times (illustrated exemplarily for lead time $4$, right panel of Figure 9), postprocessing is less successful in improving forecast quality with forecasts in summer often deteriorating for all but the DEM + Pretest method. With pretesting the seasonal cycle in quality improvements is much less apparent for lead time $4$ than for lead time $1$. This is likely due to the combination of calibration + pretesting which is performed at individual stations, which guarantees (in expectation) that the quality of postprocessed forecasts is at least as good as that of the direct model output. If, on the other hand, there is considerable miscalibration of forecasts even if only at a few stations, this can be exploited. We also detect noticeable differences in the improvements which are achieved by extending the DEM model with the Pretest between lead times $1$ and $4$: The improvement for lead time $4$ is higher in most months, especially for June to August $2017$ and August to October $2018$.

We have also examined the spatial skill of the models. Therefore, we have compared the station-wise mean CRPS of the models with the one of the raw ensemble. The skill, which is very similar at neighbouring stations, increases in the Alps and is marginal or non-existent in the Swiss plateau. The resulting spatial distribution of the skill looks similar as the mean bias depicted on Figure 4(a), the maps are therefore not shown here.

As proposed by Thorarinsdottir and Schuhen (2018), we use more than one Scoring Rule for the assessment of our postprocessing methods and apply the Brier Score to evaluate forecast quality of specific events. For precipitation forecasts, the ability to predict whether it will rain or not is of particular interest, this is captured by the Brier score for daily rainfall $<$ 0.1 mm (precision of observation measurements). We extend the assessment by also considering forecasts of moderate and heavy precipitation characterized by daily rainfall above 5 mm/d and 20 mm/d. Figure 10 illustrates the skill of the different postprocessing models by comparing the mean Brier Scores of the different models, thresholds and lead times with the ones of the raw ensemble forecasts. The assessment with the Brier Score confirms that the improvement achieved with the DEM + Pretest model is higher than the one with the other models. Only for lead times 1 and 1.5 and moderate or heavy rainfall, the local model outperforms the DEM + Pretest approach. Overall, the skill decreases with increasing threshold and increasing lead time. This is to be expected given that the postprocessing focuses on improving forecasts on average and exceedances of high thresholds are relatively rare (4$\%$ of the observed daily rainfall falls in the heavy rainfall category). Also, we use square-root transformed precipitation in the optimization which further reduces the importance of heavy precipitation events in the training set. The plot confirms further that the global model performs worst, for moderate or heavy rainfall and a lead time above 1.5 respectively 2.5 even worse than the raw ensemble. As in measures of the CRPS, the local and the DEM model score comparable for all thresholds and lead times between 2.5 and 4. For smaller lead times, the local model performs better for all thresholds.

Raw ensemble forecasts are often underdispersed and have a wet bias (Gneiting et al. 2005). This holds for the ensemble precipitation forecasts used in this study as well. Figure 11 (right) shows the verification rank histograms for the raw ensemble. Again, we use the data from June $2017$ to May $2019$ of Dataset 2 for this assessment and depict the results for lead time $1$. We note that the histogram for the raw ensemble has higher bins at the left and right marginal ranks and higher bins for the ranks which lie in the first half of $1,2,...21$. Therefore, it indicates that the ensemble forecasts are underdispersed and tend to have a wet bias. This raises the question whether the results of the DEM + Pretest model are still calibrated.

To be able to evaluate the calibration of the full predictive distribution of the postprocessing models, we do not use the reverse transformation of the square root for this assessment. Additionally, we have to use a randomized version of the PIT as our predictive distribution has a discrete component (Thorarinsdottir and Schuhen 2018):

where $V\sim\mathcal{U}\left([0,1]\right)$ and $y\uparrow Y$ means that $y$ approaches $Y$ from below.

Figure 11 shows the PIT histograms for the DEM and the DEM+Pretest model. As expected, the PIT histogram of the Pretest model lies between the one of the raw ensemble and the one of the DEM model without Pretest. The first and the last bins, which are higher than the other bins, indicate that the Pretest model is underdispersed. However, it seems much less gravely than for the raw ensemble. Since the remaining tested seasonal approaches produce worse results, this slight miscalibration of the DEM + Pretest model is a disadvantage we have to accept for the moment.

Finally, we want to get an idea of the acceptance behaviour of the DEM + Pretest model. Figure 12 shows when and where the DEM + Pretest model does (light point) and does not (dark point) postprocess the raw ensemble. Again, we focus on the results of lead time $1$. The plots show that the months where the model uses the raw ensemble lie mostly in summer and autumn. The postprocessing during these months is accepted at stations which lie in bands of different widths parallel to the Alps. The model postprocesses the raw ensemble at all stations during almost all months in between December and May, the only exception are March and April $2019$. Therefore, it appears that the Pretest approach can address the seasonal difficulties of postprocessing ensemble precipitation forecasts.