Dealing with non-stationarity in sub-daily stochastic rainfall models

Understanding the stationarity properties of rainfall is critical when using stochastic weather generators. Rainfall stationarity means that the statistics being accounted for remain constant over a given period, which is required for both inferring model parameters and simulating synthetic rainfall. Despite its critical importance, the stationarity of precipitation statistics is often regarded as a subjective choice whose examination is left to the judgement of the modeller. It is therefore desirable to establish quantitative and objective criteria for defining stationary rain periods. To this end, we propose a methodology that automatically identifies rain types with homogeneous statistics. It is based on an unsupervised classification of the space–time–intensity structure of weather radar images. The transitions between rain types are interpreted as non-stationarities.

Our method is particularly suited to deal with non-stationarity in the context of sub-daily stochastic rainfall models. Results of a synthetic case study show that the proposed approach is able to reliably identify synthetically generated rain types. The application of rain typing to real data indicates that non-stationarity can be significant within meteorological seasons, and even within a single storm. This highlights the need for a careful examination of the temporal stationarity of precipitation statistics when modelling rainfall at high resolution.

1 Introduction

Stochastic rainfall models are statistical models that aim at simulating realistic random rains. For this purpose, they generate rainfall simulations which reproduce, in a distributional sense, a set of key rainfall statistics derived from an observation dataset (Benoit and Mariethoz, 2017). The practical interest of stochastic rainfall models is notably to complement numerical weather models for the simulation of rainfall heterogeneity at fine scales, and to quantify the uncertainty associated with rainfall reconstructions. Indeed, numerical weather models face challenges for reproducing rainfall heterogeneity in space and time, in particular at fine scales (Bauer et al., 2015; Bony et al., 2015). Some of the main applications of stochastic rainfall models are therefore the fast generation of synthetic rainfall inputs for local impact studies related for instance to hydrology (Caseri et al., 2016; Paschalis et al., 2014) or agronomy (Mavromatis and Hansen, 2001; Qian et al., 2011), and the downscaling of aggregated precipitation products such as rain observations (Allcroft and Glasbey, 2003; Bárdossy and Pegram, 2016) or numerical model outputs (Vaittinada Ayar et al., 2016; Wilks, 2010). In all cases, the target is the transposition of observed rain statistics into synthetic rain simulations.

Recently, considerable attention has been paid to increasing the resolution of stochastic rainfall models so that they can mimic rainfall at sub-daily timescales. Currently, several high-resolution stochastic rainfall models are able to deal with precipitation data at typical resolutions of 1 min to 1 h in time and of 100×100 m² to 1×1 km² in space (see for example Benoit et al., 2018; Leblois and Creutin, 2013; Paschalis et al., 2013). At such scales, not only the marginal distribution of observed rain intensity matters but the space–time dependencies within rain fields are also important features of the rain process (Emmanuel et al., 2012; Marra and Morin, 2018). In particular, the impact of the advection and diffusion of spatial rainfall patterns (e.g. rain cells or rain bands) have to be modelled (Creutin et al., 2015; Lepioufle et al., 2012). As a consequence, most sub-daily stochastic rainfall models consider rainfall as a random space–time process.

An underlying hypothesis in stochastic rainfall modelling is that of stationarity: the statistics of rainfall are supposed to be constant over a given (space–time) modelling domain. This enables (1) the inference of rainfall statistics from an observation dataset, and (2) the reproduction of these statistics in simulations. The definition of stationary domains can be regarded as a modelling choice, often subjective and left to the judgment of modellers (Journel, 1993). It consists of defining pools of data that are considered similar enough (in a statistical sense) to perform model inference. In the case of stochastic rainfall modelling, the identification of stationary datasets or sub-datasets relies on some phenomenological guesses about rainfall, which serve as fuzzy guidelines to delineate stationary domains. Depending on the application and modelling choices, the parametrization of sub-daily stochastic rainfall models is considered as changing at scales ranging from seasons (Bárdossy and Pegram, 2016; Paschalis et al., 2013; Peleg et al., 2017) to single rainstorms (Benoit et al., 2018; Caseri et al., 2016).

One possible approach to delineate pools of homogeneous rain observations in a more quantitative way is to classify them prior to modelling. A set of predefined criteria is used to build a metric of similarity between the observations, and a classification algorithm is applied to the resulting similarity measures in order to define clusters of closely related rain observations. The result of such a classification procedure, often referred to as rain typing, is the identification of a limited number of rain types which gather rain observations that share similar properties. Until recently, rain typing mainly focused on classifications based on rain intensity only, with the aim to assess the physical processes responsible for rain generation (e.g. distinguish convective and stratiform rains) (Biggerstaff and Listemaa, 2000; Llasat, 2001; Rosenfeld and Amitai, 1995). In the last years, the emergence of metrics characterizing rainfall spatial or space–time behaviour (Aghakouchak et al., 2011; Ramirez-Cobo et al., 2010; Vrac et al., 2007; Zick and Matyas, 2016) paved the way for new rain typing methods based on the arrangement of rain fields in space and in time (Lagrange et al., 2018; Leblois, 2012).

In this context, the present paper focuses on the temporal non-stationarity of rainfall space–time statistics in view of sub-daily stochastic rainfall modelling. We intentionally restrict our investigation to temporal non-stationarities and consider stationarity in space (i.e. constant statistics over the whole area of interest) as a prerequisite modelling assumption. The goal is therefore to identify periods of time during which rainfall space–time statistics remain as constant as possible over a given area. The proposed framework relies on the classification of radar images based on their space–time features. The resulting classes are then used to define rain types that group rain fields with similar statistical signatures. Finally, the transition between rain types is interpreted as a break in the temporal stationarity of rainfall statistics.

The remainder of this paper is structured as follows: Sect. 2 gives a general overview of rainfall space–time patterns visible in radar images. Section 3 describes a rain typing method based on the previously identified patterns, and explains how the resulting rain types can be used to identify stationary periods. Then, Sect. 4 assesses the performance of this method for both synthetic and real case studies. Next, Sect. 5 discusses the dependence of the proposed rain typing method to the stochastic model in use, as well as the implications of the observed patterns of non-stationarity for sub-daily stochastic rainfall modelling. Finally, Sect. 6 gives general conclusions.

2 Overview of rainfall space–time patterns observed in radar images

Prior to the design of a quantitative method to identify non-stationarities in rainfall statistics, the current section seeks to illustrate with some typical examples the diversity of space–time patterns that can be observed in rain fields, and to give an overview of their temporal evolutions.

We illustrate this study with data collected over the Vaud Alps, Switzerland (Fig. 1a). The area of interest encompasses a network of high-resolution rain gauges used later for validation and covers an area of 60×60 km². For reasons of data availability, we focus in this paper mostly on summer rains observed from 1 July to 31 August 2017. During this period, only periods corresponding to rain events are considered. A rain event is defined as a rainy period isolated by at least 30 min of dry conditions before and after, and resulting in at least 2 mm of cumulated rain height (in average over the area of interest). The dataset is comprised of 17 rain events causing around 250 mm of cumulative rain height.

Figure 1b, c display different rain fields observed by the Swiss weather radar network operated by the MeteoSwiss weather agency (Germann et al., 2006). Weather radars are remote sensing devices providing comprehensive images of rain fields at the regional scale (coverage up to about 200 km from the radar device), with a high resolution (in the present case, 10 min in time and 1 km in space). The resulting rain rate estimates are known to be biased (Berndt et al., 2014), but at the same time radar images are currently the most reliable and exhaustive source of information about the spatial organization of rain fields and their temporal evolution (Emmanuel et al., 2012; Marra and Morin, 2018; Thorndahl et al., 2017): this is why radar images are used in the current section to illustrate rainfall space–time behaviour, and in the following to extract rainfall statistics.

Figure 1 Examples of rain fields over Vaud Alps, Switzerland. (a) Situation map. The area of interest is delineated by the blue square. The red star denotes the location of the rain gauge network used for validation (Sect. 4). (b) Examples of rainstorms with different space–time behaviours. (c) One example of a rain field with a temporally changing behaviour. In (b) and (c) the red arrows displayed in the last panels denote the advection of the rainstorm in 20 min.

A visual inspection of the rain fields displayed in Fig. 1 allows understanding of the benefit of jointly characterizing the spatial and the temporal behaviour of high-resolution rain fields. At the scale of interest, all rain fields shown in Fig. 1b, c have in common strong space–time interactions structured by the following:

• a distribution of intensities that is often skewed (Vrac and Naveau, 2007), with a variable amount of zero values due to within-storm rain intermittency (Mascaro et al., 2013; Schleiss et al., 2011);

• well-defined spatial patterns (Emmanuel et al., 2012; Guillot, 1999; Marra and Morin, 2018), which can be linked to the processes responsible for rainfall production such as rain cells, rain bands or rainstorms;

• a temporal behaviour shaped by the advection and diffusion of spatial patterns over time (Creutin et al., 2015; Lepioufle et al., 2012).

Despite common space–time attributes, the three rain events of Fig. 1b look very different. For example, the Type A event is characterized by a strong spatial intermittency (i.e. dry and wet locations coexist in the same radar image) combined with well-defined areas generating intense rainfall, while the Type C event shows a lower fraction of dry locations, is less spatially structured and generates low rain rates. Hence, it should be possible to find statistical metrics able to distinguish these three rain types based on their spatio-temporal characteristics.

It is worth noting that not only the space–time statistics of rainfall change between rain events, as shown in Fig. 1b, but also these statistics can change within a single rain event, as illustrated in Fig. 1c. In this case, a widespread and spatially continuous rain field (first two images) is replaced by disconnected rain aggregates (last two images). This change in the space–time features is very rapid and takes place in less than 30 min. Such abrupt changes in the space–time behaviour within a rain event are relatively common in our dataset, as discussed later.

Starting from this example, this paper investigates how to detect non-stationarities in rainfall space–time statistics using radar images as primary information. We adopt a Eulerian approach and investigate the temporal variability of rainfall statistics over a given area of interest, as perceived by an Earth-fixed observer.

3 Assessing rain statistics stationarity from radar images

3.1 Extracting space–time information from radar images

To assess the stationarity of rainfall space–time statistics, we propose to start by extracting information on the rainfall space–time behaviour from radar images. To this end, 10 statistical metrics are derived for every radar image (Fig. 2), which are split in three categories that reflect the three main characteristics of rain fields identified in Sect. 2:

• Intensity indices (IIs). These relate to the probability density function (histogram) of the rain intensities measured in a given radar image. The following indices are used:

  • II.1: fraction of the image covered by rainy pixels (informs the intra-storm rain intermittency).

  • II.2: mean rain intensity computed over all rainy pixels.

  • II.3: an 80 % quantile of rain intensities characterizing heavy rain pixels.

• Spatial indices (SIs). These characterize the spatial arrangement of patterns within rain fields. They are selected among the indices proposed by Aghakouchak et al. (2011). They are computed based on binary images representing rain masks (Fig. 2). Such binary images are obtained by thresholding radar images at a rain intensity of 0.1 mm h⁻¹, and by assigning a 0 value to the pixels under the threshold and a 1 value otherwise. Then, connected components, hereafter referred to as rain aggregates, are identified in every binary image. Their morphological properties are used to derive the following indices:

  • SI.1: fraction of rainy area covered by the largest rain aggregate in the image. This is a first indication of how the rain field is split into aggregates. Let N_p be the total number of rainy pixels in the binary image and N_m be the number of pixels of the largest aggregate; then $\mathrm{SI}.1=\frac{N_m}{N_{\mathrm{p}}}$.

  • SI.2: connectivity index. It is equal to 1 if the rain field is fully connected (one single rain aggregate) and tends to zero if the rain field is split into many disconnected aggregates. Let N_c be the number of rain aggregates in the binary image, then $\mathrm{SI}.2=\mathrm{1}-\frac{N_c-\mathrm{1}}{\sqrt{N_{\mathrm{p}}}+N_c}$.

  • SI.3: perimeter index, characterizing the sinuosity of the contours of rainy areas. It is equal to 1 if all rain aggregates are squares and tends to 0 if the rain aggregates are very sinuous. Let p be the total perimeter of rain aggregates, i.e. the sum of the perimeters of all rain aggregates, then

    $$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{SI}.3={\displaystyle \frac{\mathrm{4}\times \sqrt{N_{\mathrm{p}}}}{\mathrm{p}}}\mathrm{if}⌊\sqrt{N_{\mathrm{p}}}⌋=\sqrt{N_{\mathrm{p}}},\\ {\displaystyle }& {\displaystyle }\mathrm{SI}.3={\displaystyle \frac{\mathrm{2}\times \left(⌊\mathrm{2}\times \sqrt{N_{\mathrm{p}}}⌋+\mathrm{1}\right)}{\mathrm{p}}}\mathrm{if}⌊\sqrt{N_{\mathrm{p}}}⌋\ne \sqrt{N_{\mathrm{p}}}.\end{array}$$

  • SI.4: area index, characterizing the spread of the rain aggregates. It is equal to 1 if the radar image contains one single aggregate, and tends to zero if the rainy pixels are only in the corners of the image. Let A_convex be the area of the convex hull encompassing all the rain aggregates, then: $\mathrm{SI}.4=\frac{N_{\mathrm{p}}}{A_{\mathrm{convex}}}$.

• Temporal indices (TIs). These characterize the temporal evolution of the rain fields. They assess the advection of rainstorms over the ground as well as the temporal deformation of spatial rain patterns. Let I_t and I_t+1 be two subsequent images. In addition, let r_i, j be the cross-correlation between I_t+1 and I_t translated by a vector $\mathit{D}=i.\mathit{E}+j.\mathit{N}$ of coordinates i and j along the eastward and northward directions respectively. Finally, let r^max be the maximum correlation and ${\mathit{D}}^{\max }=i^{\max }.\mathit{E}+j^{\max }.\mathit{N}$ the corresponding displacement vector. Then the TIs are defined by the following:

  • TI.1: eastward component of the displacement vector, i.e. i^max. This index corresponds to the advection of the rainstorm along the west–east direction between times t and t+1.

  • TI.2: northward component of the displacement vector, i.e. j^max. This index corresponds to the advection of the rainstorm along the south–north direction between times t and t+1.

  • TI.3: correlation coefficient between the two radar images I_t and I_t+1 after removing advection effects, i.e. r^max. This index equals one if the spatial rain patterns remain identical between two subsequent radar images (up to a translation), and tends to zero if the images are completely different.

Figure 2 Computation of indices characterizing rainfall space–time statistics for a single radar image. This procedure is repeated for each image with >10 % rainy pixels. Note that the temporally subsequent image is required to compute the time indices (TIs).

3.2 Classification of radar images based on rainfall space–time statistics

The 10 indices defined above are used to classify the radar images in order to obtain a limited number of rain types. To ensure the reliability of these indices, only images with a significant proportion of rainy pixels are used for classification. Indeed, if the number of rainy pixels is low, the SIs are not meaningful and the TIs cannot be computed because the image correlation procedure fails. We therefore only classify the images with more than 10 % rainy pixels, the remaining rainy images (rain fraction <10 %) being typed afterwards.

To define rain types, we adopt an approach based on a Gaussian mixture model (GMM) (Fraley and Raftery, 2002; Pernin et al., 2016) classifier. This classifier has been selected because it allows for an automatic selection of the number of classes, and because it does not require any a priori information about the joint distribution of the rain indices. The idea of the GMM is to approximate the joint distribution of the 10 statistical indices. This approximation being a combination of Gaussian functions, those can be considered as representing several discrete categories. In the GMM, the joint distribution $\widehat{p}\left(\mathit{x}\right)$ of the space–time indices forming a vector x∈ℝ¹⁰ is approximated by a weighted sum of K multivariate Gaussian distributions $\mathcal{N}\left(\mathit{x}\right|{\mathit{\mu }}_k,{\mathrm{\Sigma }}_k),k=\mathrm{1},\mathrm{\dots },K$, with respective mean vector μ_k and covariance matrix Σ_k:

$$\begin{array}{}\text{(1)}& \widehat{p}\left(\mathit{x}\right)=\sum _{k=\mathrm{1}}^K{\mathit{\pi }}_k\times \mathcal{N}\left(\mathit{x}\right|{\mathit{\mu }}_k,{\mathrm{\Sigma }}_k).\end{array}$$

The inference of the model parameters (i.e ${\mathit{\pi }}_k,{\mathit{\mu }}_k,{\mathrm{\Sigma }}_k,k=\mathrm{1},\mathrm{\dots },K$) is performed with an expectation–maximization (EM) algorithm (Fraley and Raftery, 2002). A full covariance model is used for the covariance matrices Σ_k in order to take into account a possible correlation between the indices. The number K of Gaussian mixtures used in the GMM model is selected by minimization of the BIC (Bayesian information criterion) derived from EM fits computed for different numbers K (Schwartz, 1978), with the goal of selecting the GMM model resulting in the best fit while maintaining a parsimonious parametrization. Here the MATLAB Statistics and Machine Learning Toolbox has been used to fit the GMM model, with the function “fitgmdist”.

Once fitted, the GMM can be used to derive a probabilistic classification of any vector x of indices. The probability that a vector belongs to the population G_j whose distribution is the jth mixture component $\mathcal{N}\left(\mathit{x}\right|{\mathit{\mu }}_j,{\mathbf{\Sigma }}_j)$ is given by the following:

$$\begin{array}{}\text{(2)}& \widehat{p}(\mathit{x}\in G_j)={\displaystyle \frac{{\mathit{\pi }}_j\times \mathcal{N}\left(\mathit{x}\right|{\mathit{\mu }}_j,{\mathbf{\Sigma }}_j)}{{\sum }_{k=\mathrm{1}}^K{\mathit{\pi }}_k\times \mathcal{N}\left(\mathit{x}\right|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k)}}.\end{array}$$

A classification of the entire image dataset can thus be obtained by assigning to each image I the class that corresponds to the most probable mixture component (in practice, the MATLAB “cluster” function is used):

$$\begin{array}{}\text{(3)}& G\left(I\right)={\max }_j\left(\widehat{p}\right(\mathit{x}\in G_j\left)\right).\end{array}$$

As mentioned above, the classification procedure can only be applied to radar images with a significant proportion of rainy pixels (>10 %). In addition, it can be desirable to avoid high-frequency successions of rain types, for instance if the targeted stochastic rainfall model requires long-lasting pools of data to perform parameter inference. To do this, we impose a temporal persistence threshold for the rain types. To this end, all images that lead to temporal clusters of classes that do not reach a certain duration are set to unclassified. Here we use a 60 min duration threshold. After cleaning the classification, all the time steps that are unclassified (either because less than 10 % of the pixels are rainy or because the image belongs to a <60 min cluster) receive the type of the classified image that is temporally the closest (i.e. nearest-neighbour interpolation along the time axis). The complete rain typing framework is summarized in Fig. 3.

Figure 3 Rain typing framework.

Since the rain typing method presented above aims at defining groups of rain fields sharing similar statistical signatures, the transitions between rain types can be interpreted as non-stationarities. Similarly, periods with a constant rain type are interpreted as stationary periods.

4 Validation and application

In this section, we validate our rain typing approach in the context of stochastic rainfall modelling. The validation study comprises four steps: first, the proposed approach is tested in a synthetic case in order to determine whether we are able to identify known non-stationarities. Then, real data are used to compare our rain typing strategy with two alternative hypotheses of rainfall stationarity: (1) rainfall statistics are stationary at a seasonal scale and (2) rainfall statistics are stationary at a rainstorm scale. Next, the rain typing method is applied to radar data covering the full year 2017, in order to assess its ability to handle various rainfall situations. Finally, the sensitivity of the classification to the calibration dataset is assessed by comparing rain types computed for the summer of 2017 based on two calibration periods: (1) the summer of 2017 and (2) the full year 2017. Prior to the validation itself, the next subsection describes the stochastic rainfall model used for validation.

4.1 Stochastic rainfall model

The validation of the rain typing approach uses a stochastic rainfall model designed for local-scale (area of a few square kilometres) and high-resolution (up to 1 min) data. This model involves 11 parameters and aims at modelling both the marginal distribution of observed rain intensities and the space–time dependencies that exist within rain fields. It is briefly introduced hereafter; for more details the reader is referred to Benoit et al. (2018). In this model, the marginal distribution of rain rates is accounted for by considering that rain measurements (R) originate from the censoring and power transform (involving parameters $a_{\mathrm{0}},a_{\mathrm{1}},a_{\mathrm{2}}$) of a standardized multivariate Gaussian random field (Z) tainted by an additive measurement noise ($\mathit{ϵ}\sim \mathcal{N}(\mathrm{0},{\mathit{\sigma }}_{\mathit{ϵ}})$) (Eq. 5):

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{ll}R={\left({\displaystyle \frac{Z+\mathit{ϵ}-a_{\mathrm{0}}}{a_{\mathrm{1}}}}\right)}^{\frac{\mathrm{1}}{a_{\mathrm{2}}}}& \mathrm{if}Z+\mathit{ϵ}>a_{\mathrm{0}},\\ R=\mathrm{0}& \mathrm{if}Z+\mathit{ϵ}\le a_{\mathrm{0}},\end{array}\end{array}$$

The multivariate Gaussian latent random field (Z) is characterized by an asymmetric Gneiting space–time covariance (Gneiting, 2002) function ρ which accounts for both the advection and the diffusion of spatial rain patterns (Eq. 6). For two rain observations separated by a spatial lag h and a temporal lag u, the covariance is given by the following:

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathit{\rho }(\mathit{h},u)={\displaystyle \frac{\mathrm{1}}{(u/d{)}^{\mathrm{2}\mathit{\delta }}+\mathrm{1}}}\exp \left({\displaystyle \frac{-\left(\right||\mathit{h}+\mathit{V}.u|\left|\right)/c{)}^{\mathrm{2}\mathit{\gamma }}}{{\left((u/d{)}^{\mathrm{2}\mathit{\delta }}+\mathrm{1}\right)}^{\mathit{\beta }\mathit{\gamma }}}}\right).\end{array}$$

In this model, the advection of rainstorms is assumed to be constant and linear along a vector V defined by its amplitude V_S and direction V_θ. The regularity parameters γ (for space) and δ (for time) control the slopes at the origin of the covariance function and thereby regulate the small-scale variability of the rain fields, and ultimately their smoothness. The scale parameters c (for space) and d (for time), in units of distance and time respectively, control the decorrelation distances of rain patterns. Finally, the separability parameter β controls the space–time interactions. When β=0, the covariance function is space–time separable.

4.2 Detection of rainfall non-stationarity in a controlled setting

The ability of the rain typing method to detect possible non-stationarities is tested by applying it to synthetic time series of radar images. These images are generated using the stochastic rainfall model presented above, with model parameters changing abruptly. This produces (temporal) non-stationary synthetic rain fields. The rain typing method is then applied to the simulated radar-like images (resolution: 1×1 km² in space, 10 min in time; footprint: 60×60 km²) in order to assess if it is able to retrieve the prescribed patterns of temporal non-stationarity.

For generating the synthetic images, we use the stochastic rainfall model described in Sect. 4.1 with model parameters corresponding to three typical rain behaviours identified by visual inspection (Table 1). Figure 4a shows an example of simulated rain field for each rain type.

Table 1 Parameters of the stochastic rainfall model used for the generation of synthetic images.

Figure 4 Identification of non-stationarities in rain statistics for a synthetic case study. (a) Examples of synthetic rain fields simulated for each rain type. (b) Segmentation of the time axis into periods with stationary rain statistics. First row: reference. Rows 2–4: segmentation obtained by rain typing applied to synthetic radar-like images. Dotted black lines represent dry periods.

In Fig. 4b, the series titled “reference” shows the rain types prescribed to the stochastic rainfall model for the generation of the synthetic radar images. Based on stochastic simulations, three sets (realizations) of synthetic images are generated. Each realization is classified to determine whether it is possible to identify the reference rain types. Results show that the proposed method can consistently detect the prescribed rain types and their temporal evolution, with an agreement between the reference and the estimated rain types of 97 %, 96.3 % and 92.6 % for realizations 1 to 3. It also properly estimates the number of rain types prescribed in the reference. The only noticeable difference between the reference and the simulations is the emergence of a very infrequent fourth rain type (accounting for 0 %, 1 % and 1.6 % of the estimated rain types for realizations 1 to 3) at the beginning or at the end of some rain events (in green in Fig. 4b). This is because at these periods, the rain does not cover the whole area of interest, and in certain situations it can produce rain fields with different space–time statistics, which induces this artificial fourth rain type. Except for this fourth rain type, results show that in this synthetic experiment, the proposed method performs well in detecting non-stationarities of rainfall space–time statistics and, in turn, periods during which these statistics remain stationary.

4.3 Impact of rainfall non-stationarity on stochastic modelling of an actual dataset

To further validate our rain typing method, we apply it to a real dataset acquired in the Vaud Alps, Switzerland, during the summer of 2017. In such a real case study, the true succession of rain types is obviously unknown. To assess the performance of the proposed rain typing method, we compare it with two other hypotheses of stationarity that can be found in the literature. We therefore consider three cases, illustrated in Fig. 5:

• Hypothesis H1. The time axis is broken down into rain types interpreted as stationary time periods (the approach proposed in this paper). Applying the rain typing method presented in Sect. 3 to the period of interest leads to six rain types.

• Hypothesis H2. The statistical signature of rainfall is constant over meteorological seasons (Bárdossy and Pegram, 2016; Paschalis et al., 2013; Peleg et al., 2017). This leads to one stationary pool of rain events for the period of interest.

• Hypothesis H3. Rainfall statistics are constant within a single rainstorm but change between storms (Benoit et al., 2018; Caseri et al., 2016). Here 17 rain events are identified following the definition adopted in Sect. 2.

Figure 5 Observation dataset used for validation. (a) Measurement network. Red dots denote rain gauge locations. (b) Picture of a rain gauge with the Vallon de Nant catchment in the background. (c) Observed rain rate averaged over the network. (d) Segmentation of the time axis into stationary periods for the three tested hypotheses. For each hypothesis, segments with the same colour denote periods for which rainfall is expected to have similar space–time statistics.

To compare these three hypotheses, we apply the same stochastic model as above to rain data collected by a dense network of eight high-resolution rain gauges set up in a small (3×6 km², Fig. 5a) alpine catchment called “Vallon de Nant”, situated within the area of interest presented in Sect. 2. Hence, in the following, radar images will be used only to carry out the rain typing presented in Sect. 3 in order to define the hypothesis H1. The remainder of the validation, i.e. stochastic model calibration and simulation under the three tested hypotheses, will be carried out on rain rate time series acquired by rain gauges, and not on radar images. The goal is indeed to keep the following validation as independent as possible from the radar images used for rain typing. By doing so, we seek to prove that the proposed method captures the stationarity of the rainfall process itself, and not only the stationarity of radar images.

Once the periods of stationarity have been built for each of the three hypotheses, the stochastic model is calibrated for each stationary period. This means that for each hypothesis, a set of model parameters is inferred from observations for each postulated stationary dataset. Then, synthetic rain fields are generated by unconditional simulation under the three hypotheses of stationarity, and in each case 50 realizations (i.e. 50 simulated synthetic rain histories) are compared to actual measurements. To assess the realism of the different scenarios, Fig. 6 shows the simulated cumulative rain heights. Next, Fig. 7 shows quantile–quantile (q-q) plots for four statistics selected to assess the marginal distribution of rain rates and its space–time arrangement: number of rain gauges measuring zero rain at each time step, rain intensity, standard deviation (in time) of rain intensities separated by a time lag of 5 min, and standard deviation (in space) of rain intensities at each time step.

Figure 6 Reproduction of the cumulative rain height (averaged over the whole network) for the three tested hypotheses. (a) Rain is stationary by rain type (H1), (b) rain is stationary during the entire summer period (H2), (c) rain is stationary by rain event (H3). The red line corresponds to the observations and the black lines correspond to simulations.

Figure 7 Reproduction of rainfall statistics for the three tested hypotheses. (a–d) Rain is stationary by rain type (H1), (e–h) rain is stationary during the entire summer period (H2), (i–l) rain is stationary by rain event (H3). (a, e, i) Quantile–quantile (q-q) plot of simulated vs. observed rain intermittency, (b, f, j) q-q plot of simulated vs. observed rain rate, (c, g, k) q-q plot of simulated vs. observed temporal variability of rainfall (at lag 5 min), (d, h, l) q-q plot of simulated vs. observed spatial variability of rainfall. The quantiles used in the q-q plots are centiles. Each centile is denoted by a black cross.

Results show that H1 tends to slightly underestimate the cumulative rain due to an underestimation of very high intensities. This underestimation of heavy rainfall is common to all the three cases and probably originates from the stochastic model itself, which is not designed to handle extreme rainfall due to the simple transform function selected in Eq. (5). This could be improved by adopting a transform function accounting for extreme rainfall (see for example Vrac and Naveau, 2007) but at the price of a more complex parametrization, which is not regarded as essential here because the observed rain rates are mostly low to moderate, and only the 99th centile is affected by rain rate underestimation. Apart from this underestimation of high rain rates, hypothesis H1 allows the other metrics to be correctly reproduced.

Contrary to H1, hypothesis H2 leads to a slight overestimation of the simulated rain height, in particular for the first 30 days (Fig. 6). This is due to the overestimation of moderate rain rates that compensates for the underestimation of extremely high values. This bias in the simulated marginal distribution is due to the lack of flexibility of H2 that imposes a single underlying stochastic model for the entire summer period. This does not allow enough flexibility to capture the diversity of structures emerging from high-resolution data. This is also visible for the simulated variability in space and time, which tends to be overestimated for the low centiles but underestimated for the high centiles.

Under hypothesis H3, simulation results are close to those of H1, with a greater propensity of underestimating heavy rainfall. In addition, the standard deviation in time is not perfectly reproduced for the middle quantiles. This slightly lower performance of H3 compared to H1 can be attributed either to the non-inclusion of intra-storm non-stationarities in this hypothesis, leading to a poorer reproduction of the true rainfall dynamics, or to a poor inference of model parameters in case of short rain events caused by the low amount of observations available for such very short stationary periods.

To sum up, the proposed method consisting of typing rain fields according to their space–time statistical signature derived from radar images (H1) leads to more realistic rainfall simulations than the other approaches, H2 and H3.

4.4 Seasonality of rain type occurrence

The two previous sections have shown that the proposed rain typing method is able to reliably identify rain types when applied to summer rains. To complement the previous findings, the current section investigates the ability of our rain typing framework to classify rainfall from other seasons. To this end, rain typing is performed for the same area of interest as above (Fig 1a), but this time for a complete year of radar observations. Here no comparison data are available since the rain gauge network used for validation in Sect. 4.3 had to be removed after the summer due to harsh local winter conditions. For this reason, this section concentrates on a qualitative analysis of the timing of rain type occurrences throughout the year, with interpretation of some resultant rain types.

The classification framework presented in Sect. 4.3 is applied to radar data covering the entire year 2017. It should be noted that rainfall and snowfall are processed without distinction, and thus the resulting classification produces precipitation types rather than rain types strictly speaking. However, to be coherent with the rest of the paper, we will continue to refer to rain types despite possible mixes between rain and snow during the winter months.

The automatic selection of the number of Gaussian mixtures used in the GMM model leads to 11 rain types for this year. Figure 8 shows the monthly occurrence of these rain types, and Fig. 9 shows the related marginal distributions of the 10 indices used to characterize rainfall space–time behaviour.

Figure 8 Monthly occurrence of rain types during the year 2017.

Figure 9 Marginal distributions of the indices used for rain type classification. Each graph represents a different index. Colours denote rain types.

The monthly occurrence of rain types shows a clear seasonality for most rain types (Fig. 8), with some types occurring mostly during winter months (types 1, 6 and 7), others during summer months (types 2, 3 and 5 ) and one during spring and fall (type 4). By contrast, rain types 8 to 11 seem relatively unseasonal. The seasonality observed in rain type occurrence is in good agreement with the local climatology of rainfall with more snowfall and stratiform rainfall expected during winter, more sleet showers and rain showers expected in spring and autumn, and more thunderstorm induced convective rainfall expected during summer. Of course there is no one-to-one match between the rain types identified by our classification and the physical types listed above, for two reasons: first, the rain typing method is purely statistical and we cannot expect a direct match with the physical processes responsible of rainfall generation. All the attempts at linking the rain types with physical processes made in this section should therefore be regarded as qualitative interpretations instead of well-defined relationships. The second reason for the imperfect match is that one rain generation process (e.g. a thunderstorm generating convective rainfall) can lead to several rain types (e.g. in this case types 2, 3 and 5).

The distributions of the space–time indices displayed in Fig. 9 allow refinement of the interpretation of the rain types. As an illustration, one rain type typical of each season will be described hereafter in light of the distribution of indices resulting from the classification:

• Rain type 6 is typical of winter months, and probably corresponds to stratiform rainfall. The rain fields classified in this type are featured by a large fraction of the area of interest covered by rainfall, which leads to a small number of large rain aggregates, and in turn to high connectivity and area indices. Such rainstorms are moving eastward and are well correlated in time, which is typical of stratiform rains over Switzerland. Finally the resulting rain intensities are low.

• Rain type 4 is typical of spring and autumn months, and probably corresponds to sleet or rain showers. The corresponding rain fields are very scattered (low fraction of rain coverage and low area index) and are poorly connected. In addition, it is interesting to note that such rain fields have a low temporal correlation, which reflects a strong variability in time, as is the case for fast-changing mid-season events.

• Rain type 5 is typical of summer months, and probably corresponds to heavy convective rainfall produced by thunderstorms. Indeed, the rain fields classified in this type generate localized (moderate fraction of rain coverage and contribution of main aggregate) but heavy (high mean and Q80 intensity) rainfall. Also, such storms mostly originate from the south-west (eastward and northward advections >0), which is typical for large thunderstorms over the Swiss Alps.

It is important to note that even if different rain types dominate at different months, several rain types occur at each month. In other words, at a given location, there is no single typical summer rain or winter rain, but rather a variable collection of rain types that occur each month. It follows for instance that the rain type 6 is typical of winter months, but it can also occur during summer, although at a lower frequency than the typical summer rain types 2, 3 and 5.

4.5 Sensitivity of the rain typing approach to the size of the calibration dataset

To complete the assessment of our rain typing strategy, we study the sensitivity of the classification method to the dataset used to calibrate the GMM model. To this end, the radar images corresponding to the periods covered by the 17 rain events of interest occurring during summer 2017 are typed based on two different GMM models:

• A first GMM model (model A) is calibrated using the radar images of summer 2017 (July and August). This corresponds to the GMM model used in Sect. 4.3, and this classification results in 6 rain types that all occur during the period of interest.

• A second GMM model (model B) is calibrated using the radar images of the full year 2017. This corresponds to the GMM model used in Sect. 4.4. This model results in 11 rain types, 9 of which occur during the period of interest (rain types 7 and 11 are absent).

Figure 10 Correspondence between rain types derived from two distinct calibration datasets.

Figure 10 shows how the results of rain typing compare when carried out based on both GMM models A and B. Since the number of rain types differs between the models, there is no one-to-one correspondence. This is because the calibration of model B on a larger dataset allows for a more detailed segmentation of the indices into rain types. Indeed, model A focuses only on the summer, and if some typical winter rain types occasionally occur during this period (which is the case as shown in Sect. 4.3), these rain fields are too infrequent to generate their own cluster. They are then assigned to the closest Gaussian mixture. In contrast, the same rain fields have proportionally more relatives in model B, and they can therefore be grouped in a separate rain type. This distinct segmentation of the indices space shows that one cannot unequivocally associate a type derived from model B to a type derived from model A. A rain type in model B can straddle the border between two rain types in model A.

Although imperfect, the links between the two classifications are strong. To assess these links, we first make coincide the two classifications by assigning to each rain type of model B its most common counterpart in model A. In this case, the rain types 1, 2, 3, 4, 5, 6, 8, 9, 10 from model B are paired with respectively rain types 1, 3, 6, 4, 3, 6, 6, 4, 1 in model A. By doing so, we obtain a 73.6 % match between the two classifications. Then, we compare the timing of the intra-event rain type transitions for the two raw classifications (i.e. without the previous pairing). Due to the unequivocal links between the two outputs, model B leads to a more fragmented classification (26 transitions instead of 16 for model A). Despite this, most of the transitions in model A (9 among 16) have a counterpart in model B (i.e. a transition appears in model B during the same hour than the transition in model A). This tends to confirm that the intra-event non-stationarities identified in Sect. 4.3 are well defined, and that their detection is robust.

5 Discussion

5.1 Model dependence in rain typing

The aim of this work is to develop a rain typing strategy able to identify stationary periods for further stochastic rainfall modelling at a sub-daily resolution, with an emphasis on very high resolution models (up to 1 min resolution). Although the main ideas developed throughout this paper can be applied to many different stochastic rainfall models, the detailed settings of the rain typing strategy must be tuned from case to case in order to be compatible with the targeted rainfall model. Indeed, because the final aim is to identify time periods over which stationary statistics can be inferred, the rain typing method necessarily depends on the properties to be modelled and the nature of the model. Therefore, we discuss hereafter how the stochastic rainfall model used in this study has influenced the settings of the rain typing method.

As primary data source, raw radar images have been preferred to combined rain gauge–radar products (Sideris et al., 2014) because their higher temporal resolution (10 min instead of 60 min) is more in agreement with the resolution of the stochastic model used (i.e. 1 min). Radar images acquired at higher temporal resolution (e.g. 5 min) are acknowledged to potentially improve the classification, but such images were unfortunately not available for this study. In counterpart of the high temporal resolution of raw radar products, the observed rain intensities can be impaired by local biases, in particular in mountainous regions such as the one considered in the present application. However, we think that our method is not significantly affected by such biases because radar shadow effects due to topography are constant in time. Therefore, they do not affect the temporal variability of the rain indices used for classification.

Once the radar product has been selected, another important choice is the size of the area of interest from which the information will be extracted. Since the study area is very small (3×6 km²), we could in theory have taken a window of the same size to analyse radar images. However, a significantly larger window (60×60 km²) has been chosen because the space indices used for classification rely on the count and on the spatial arrangement of rain aggregates, which requires several aggregates within the window of interest for reliable computations. This implies a window slightly larger than the expected scale of rain aggregates, which ranges from 10 to 30 km in the present case (Benoit et al., 2018). On the other hand, the window extent is restricted by the wish to limit the variety of rain behaviours existing in that window, thereby preserving as much as possible the stationarity of rainfall in space.

Finally and more importantly, the indices selected for classification must be consistent with the statistics embedded in the stochastic rainfall model. Here 10 indices are needed to characterize the evolution of the main features of rainfall considered by the stochastic model, namely the marginal distribution of rain intensity, the spatial arrangement of rain aggregates, and the advection–diffusion of rainstorms. This large number of indices reflects the complexity of the stochastic rainfall model in use. In case of a different model (due to for example a single-site study, a different resolution or another climatology) the set of indices could be modified.

5.2 Consequences of non-stationarity on sub-daily stochastic rainfall modelling

The succession of rain types identified in Sect. 4.3 and 4.4 has two important implications on how stationarity should be regarded in sub-daily stochastic rainfall models.

Due to the variability observed in rainfall statistics, we believe that it is often incorrect to assume rain stationarity over long periods, such as months or seasons, for characterizing of rainfall statistics. In the perspective adopted in this paper, the seasonality observed in monthly rain statistics is attributed to a variation in rain type occurrence rather than to a smooth change in rainfall behaviour (see Sect. 4.4). As such, it would be incorrect to assume constant model parameters over long periods of time because of the risk of mixing distinct rainfall statistics during model calibration. This results in the emergence of artificial rain types whose statistics are an average of the statistics of the rain types that actually occur during the period of interest. In the extreme case, such artificial rain types may not even correspond to any actual rainfall event. This leads in turn to improper space–time dependencies in the simulated rain fields.

A second striking observation is that rainfall statistics can change drastically within a single rain event. As a result, the hypothesis of rain stationarity along entire rain events can be invalidated in some instances. At least in our dataset, such non-stationary events seem relatively frequent (at least 7 non-stationary rain events out of 17 have been identified in our data). This observation is not new, since it may lead to temporal asymmetry (or temporal irreversibility) in rain rate time series, which is discussed by Müller et al. (2017). Our framework offers a way to deal with this phenomenon through the identification of stationary periods prior to the stochastic modelling of rain. Then, stochastic modelling is carried out separately for each rain type, and the results are merged afterwards. This allows the generation of synthetic rain fields presenting a temporal asymmetry, even if the stochastic model itself is only capable of generating symmetric rain fields for a given set of model parameters. The temporal asymmetry is then carried by the temporal arrangement of rain types within a single rainstorm.

6 Conclusions

This paper proposes a quantitative method to identify stationary rainfall periods, that is, periods during which a set of 10 statistics representative of rainfall space–time behaviour at the local scale remains broadly constant. It is based on a classification of radar images into groups of rain fields sharing similar statistics when observed at high resolution. For reasons of data availability, we focused our investigation on summer rains over the Vaud Alps, Switzerland. However, most of the results obtained in this context are expected to be transferable to other seasons, as illustrated in Sect. 4.4, as well as to other mid-latitude areas in cases where extratropical rainstorms significantly influence the precipitation regime.

The application of the proposed rain typing method in the context of sub-daily stochastic rainfall modelling shows that our method is able to (1) identify abrupt changes in rainfall statistics during a controlled synthetic experiment, and to (2) delineate relevant stationary periods when applied to an actual dataset. The succession of rain types obtained for our observation dataset is characterized by the coexistence of several distinct rain types, with switches between types occurring within single rainstorms. In the context of sub-daily stochastic rainfall modelling, this observation highlights the need to delineate stationary periods based on actual observations rather than subjective assumptions about the rain process.

A possible future work would be to use the proposed method to simultaneously type precipitation fields observed by radars and simulated by numerical models. This may provide a new metric to assess the precipitation component of high-resolution numerical weather or climate models. Indeed, the proper reproduction of rain types and rain type successions in model outputs would indicate a correct simulation of the overall space–time behaviour of rainfall by the model.

It could also be interesting to apply the proposed rain typing method to long-term archives of radar images in order to investigate the temporal behaviour of rain type occurrence. The resulting information could be used as the starting point for the design of a statistical model of rain type occurrence, and in turn a stochastic rain type generator. Coupled with already existing high-resolution stochastic rainfall models, this would allow the design of high-resolution stochastic rainfall generators that are able to reproduce local rainfall statistics over long simulation periods under the assumption of a steady climate. An alternative to the development of a stationary rain type generator would be the design of a synoptically conditioned stochastic rainfall generator (Bárdossy and Plate, 1992; Peleg and Morin, 2014) by linking the occurrence of rain types with the state of the atmosphere simulated by one or several climate models. This could be achieved by assessing the statistical relationships between rain types derived from radar observations and meteorological variables derived from climate model reanalyses. If strong dependences are found, then the future evolution of climate variables simulated by general circulation models could provide precious insights into the possible evolution of rain type occurrence in a changing climate.