Modelling the frequency distribution of inter-arrival times from daily precipitation time-series in North-West Italy

A robust statistical inference of the alternation of wet spells and dry spells (DS) is still an open issue, despite its fundamental role in several common hydrological applications, including irrigation management (Rallo et al. 2014; Canone et al. 2015), the analysis of precipitation temporal trends (Agnese et al. 2012) and the eco-hydrological modelling of land–atmosphere feedbacks (Agnese et al. 2016).

In most studies, the statistical analysis of the sequences of rainy days, previously denoted wet spells (WS), and that of no-rainy days, DS, was carried out separately. In particular, the observed WS frequencies were commonly modelled by the geometric distribution, consequent to the statement that a rainy day is a ‘success’ outcome of classical Bernoulli trials, with constant rain probability.

As showed by Chatfield (1966), this assumption of constant rain probability is commonly unrealistic for modelling DS frequency distribution; in fact, as stated by the author, ‘the probability that a dry day will be followed by a dry day does increase with the previous number of consecutive dry days’. By using daily rainfall data of London from 1958 to 1965, Chatfield successfully fitted the dry spell frequencies through the logarithmic-series distribution; this distribution allows accounting for the observed persistence of dry states, because of its peculiar pattern having a slowly decaying tail (Kottegoda & Rosso 1997). In a previous study, Agnese et al. (2014) demonstrated that the two-parameter polylogarithmic distribution significantly outperforms the logarithmic distribution, in order to statistically infer the observed DS frequencies.

As mentioned above, here the joint modelling of both WS and DS is performed, instead of treating them separately. In fact, a single occurrence process regulates the alternating sequences of rainy and no-rainy days. A similar approach was carried out by Strupczewski et al. (2013), who derived a two-component model of flood dynamics that combines the duration of high waters and the exceeding of a discharge threshold, thus distinguishing two possible states, namely, flood and no-flood.

A joint modelling of both WS and DS can be obtained by considering the observed frequencies of the so-called ‘inter-arrival times (IT)’ of precipitation events, ITs, which represent the series of times elapsed between two subsequent rainy days, instead of the commonly used number of dry days. Indeed, the probabilities of rainy day sequences (i.e., wet spells, WS) and of dry sequences (i.e., dry spells, DS), can be easily derived from the IT distribution. IT series can be viewed as a realization of a renewal process (Feller 1968; Cox 1970; Buishand 1977), by assuming that ITs are independent and identically distributed (i.i.d.) random variables. Accordingly, the WS geometric distribution directly emerges from the renewal property, whereas the DS probabilities can be easily calculated by rearranging the IT probabilities.

Another important aspect of this work is to use a discrete probability distribution. In the literature, both discrete and continuous probability distributions were used to model the frequency distribution of inter-arrival time-series (Foufoula-Georgiou & Lettenmaier 1986; Rodriguez-Iturbe et al. 1987; de Groen & Savenije 2006; Manfreda & Rodrıguez-Iturbe 2006), although the inter-arrival time is properly a discrete random variable. In fact, Foufoula-Georgiou & Lettenmaier (1986) demonstrated the statistical advantages of treating the rainfall occurrence as a discrete process, rather than a continuous one. In order to describe the rainfall duration, Manfreda & Rodrıguez-Iturbe (2006) considered the exponential distribution, which is the extension to the continuum of the discrete geometric distribution. Eichner et al. (2011) discussed the influence of the persistence of the daily rainfall process on that of the IT by using long-term correlated data sets.

In a previous work, the ITs have already been modelled by Agnese et al. (2014) using the discrete three-parameter Lerch distribution. This distribution is able to account for some peculiar features usually observed in IT samples, which other distributions are not capable of reproducing as well. These features are the high standard deviation and the high positive skewness, the monotonically decreasing pattern of frequencies with a slowly decaying tail and a distinctive ‘drop’ marking the passage from the wet status (IT= 1) to the dry status (IT= 2). The ability of the Lerch distribution in reproducing all of these features allows the use of this distribution for a variety of hydrological applications. In a previous work regarding the period 1975–2005, Agnese et al. (2014) successfully applied the Lerch distribution to the IT frequencies derived from 55 daily rainfall series recorded in Sicily.

However, rainfall regime in Sicily is characterized by a strong seasonality, splitting of the year into two parts, a wet cold season (from October to March) and a dry warm one (from April to September). Consequently, data have been preliminarily partitioned into two groups falling in those two different periods of the year, and then they were fitted separately via the Lerch distribution. For data belonging to the ‘wet’ semester, all records satisfied the chi-square test of the goodness of fit; whereas for data of the ‘dry’ semester, 11 series did not satisfy the test, likely because of the small available sample size.

The primary objective of this work is to test the Lerch distribution in an area ranging from the highest peaks of the Alps to the Tyrrhenian Sea. In order to do that, the IT series were derived from long daily precipitation series (70/90 years duration), collected in 26 stations located throughout the Piedmont and Aosta Valley, in typical Alpine and sub-Alpine areas of North-West Italy. The climate of the studied area is strongly spatially variable (Figure 1) and was classified as sub-litoranean, sub-alpine, pre-alpine and sub-continental, and differs completely from that of Sicily which is typical Mediterranean. The different parameter values derived by the statistical inference were compared to those obtained in the study about the Sicily region and results were discussed. Therefore, the application of the three-parameter Lerch distribution also to the area of NW Italy represents an important test to establish the wide applicability of the proposed approach. An attempt to explore the validity of the hypothesis (usually accepted) that ITs are i.i.d.s was carried out.

Moreover, an objective of this paper is also to use the Lerch distribution in order to describe the sequences of rainy days, interrupted by 1-day DS, which could have a certain interest from a hydrological point of view, representing in a broad sense, a unique perturbed period, not significantly altering the antecedent soil moisture status. The equivalence on using a single IT three-parameter Lerch distribution, in place of using two distributions, namely DS two-parameter polylogarithmic distribution and WS one-parameter geometric distribution, is also argued.

Finally, the expected return periods of dry and wet spells at the regional scale are determined, in order to quantify the spatial variability of the precipitation regimes in the investigated area, and to provide suitable results for the sake of agricultural water management.

In this work, the daily precipitation data long series (70/90 years) from 26 stations have been used. They belong to the available databases of the National and Regional Observatories (ARPA Piemonte and Centro Funzionale Valle d'Aosta), and of the Italian Meteorological Society (Mercalli & Cat Berro 2008). These stations are located throughout the Piedmont and Aosta Valley area (Italy) and cover an area of 28,000 km², characterized by different climatic regimes, ranging from sub-continental to sub-litoranean ones, in between 113 and 2,170 m altitude. It is an area characterized by very high mountains like Mont Blanc, Matterhorn and Gran Paradiso, quite close to the plains and to the sea. In addition to the mentioned orographic issue, the irrigated agricultural Po Plain plays a major role in the air humidity pattern because of the widespread surface irrigation, which saturates wide areas. This produces a notable variation in the average statistics of precipitation occurrence through the year and through the study area.

The stations were selected in order to have highly reliable and long records, with few gaps in the data series. Figure 1 shows the location of the stations (dots) on a physical map of the area. Table 1 reports the available years of observations, the elevation, the latitude and the longitude (UTM 32 T, m). For the purpose of this study, a common starting year (the year 1919, if available) was considered for all 26 data series. The IT series were obtained from the recorded precipitation observations using the precipitation threshold h* = 1 mm (Figure 2), which is usually adopted in Italy by regional observatories in order to discriminate between rainy and no-rainy days.

A simple analysis of the monthly-averaged frequencies of wet days through the year (Table 2), revealed the presence of a ‘wet’ period in spring (peak in May) and a ‘dry’ period in winter (minimum in January), for all the stations. Furthermore, another ‘wet’ period in autumn is observed in the area located at low–medium altitudes, whereas a secondary peak in August is present at the high altitudes. These facts complicate the splitting of the year into two homogeneous periods valid for the whole area, as for the Sicily case.

In spite of these problems, with almost similar characters for all the sites, a simple partition of the whole year into two contiguous periods was chosen. Accordingly, a first period, denoted as the ‘cold season’, including the four months (from December to March) in which lowest precipitation day frequencies are observed ‘on average’ and a second period, namely, the ‘warm season’, covering the remaining months (from April to November), were selected. This choice is supported by Figure 3, which reports, for all the considered stations, the average number of events per month, <# events/month>, during the ‘cold season’ and the ‘warm season’. As can be observed, all the months of the ‘cold season’ show lower average number of events per month than the months of the ‘warm season’.

Then, for each season, a statistical inference was performed by fitting the IT data frequencies through the three-parameter Lerch (Lch) distribution. For the sake of comparison, the statistical inference was also applied to the whole year data set (Y).

Before performing the statistical inference, the occurrence of a change-point in some average annual IT data series, <IT>, was detected. For this purpose, a number of different tests can be applied (Baiamonte & D'Asaro 2016; Pohlert 2018), such as either the Buishand test or the standard normal homogeneity test. However, for its simplicity and for its wide application in hydrology (Zhang et al. 2009), the Pettitt test (Pettitt 1979) was considered in this study. The Pettitt test localizes the year when most likely a breakpoint (discontinuity or inhomogeneity) occurs. Detecting the change-point allows dividing the entire data series into two intervals, before and after the change-point, forming two homogeneous groups.

For changing in one direction, the statistics K_T⁺ = max {U_t,T} and K_T⁻ = −min {U_t,T} were calculated, with K_T= max {K_T⁺, K_T⁻}, indicating that large K_T⁺ corresponds to a shift down in level, whereas large K_T⁻ corresponds to a shift up in level from the beginning of the series, respectively.

Table 3 reports the statistics K_T⁺, K_T⁻, K_T, the p-value significance level, the year of changing and the shift upward or downward for the considered 26 stations.

The results showed that the null hypothesis was rejected (i.e., no change-point occurs) only for eight stations out of the selected 26 ones. Thus, 18 data series were considered as a whole, whereas in the remaining eight series, the period before the changing year (detected by the Pettitt test) was eliminated. In those cases, only the more recent homogeneous period was considered for the application of the Lerch distribution.

The cases where the null hypothesis was rejected are depicted in Figure 4, showing the <IT> time-series, the location of the change-point, and the disregarded period highlighted by a shaded area. As can be observed, the upward or downward shift agrees with the results of the Pettitt test reported in Table 3. However, it should be emphasized that, as recently observed by Serinaldi et al. (2018), such change-points should be accepted as physically meaningful if they can be related with a predictable process based on theoretical models and/or well identified physical dynamics justifying causality. In our case, neither of the occurrences was possible to find.

Let be a time-series of precipitation data spaced at a uniform time-scale, τ, e.g., day; a day k is a ‘rainy’ day if ⁠, otherwise k is a ‘no-rainy’ (or dry) day; is any precipitation threshold.

From the H series (of size N), the inter-arrival time series, ⁠, representing the succession of the times elapsed between two subsequent rainy days, could be easily derived; IT size, M, equals the number of rainy days. Accordingly, an isolated rainy day is identified by the conditions IT_j > 1 and IT_j+1 > 1, whereas a wet spell, WS, of duration n+ 1 is identified by any sequence, in IT, of n consecutive 1 values. Moreover, by neglecting any IT_k = 1, the DS sequence, can be easily obtained, by DS_j = IT_j− 1; then, if r indicates the frequency of ITs equal to 1, the size L of DS (or, equivalently, of WS) results as: ⁠. Figure 2 reports an example of the terms and sequences defined above.

In order to reproduce the main features of the above-mentioned empirical IT distributions, Agnese et al. (2014) suggested applying the three-parameter Lerch distribution.

In the most general form, the distribution reported in Equation (5), hereafter denoted as Lch, depends on the three parameters θ, s and a (Table 4). Positive s values ensure that the probability mass function (pmf) is a monotonically decreasing function with mode X= 1.

Therefore, the distributions expressed by Equations (8) and (9) are different forms (particular cases) of Lch distribution.

In Table 4, the Lerch family of the probability distributions and the corresponding parameters domain are resumed (ID = 1–7) (Zörnig & Altmann 1995).

In order to perform the statistical inference of the previously introduced pmfs, the well-known maximum likelihood estimation (MLE) method was used. The analytical solution of the MLE method for the general three-parameter Lch implies the equality of the arithmetic, geometric and harmonic expectations to the corresponding sample means (Gupta et al. 2008; Agnese et al. 2014). Obviously, two-parameter and one-parameter distributions of the Lerch family do not fulfil all equality conditions; in Table 4, conditions that are satisfied or not satisfied by each distribution are also indicated.

The classical χ² procedure is not appropriate to test the goodness-of-fit of long-tailed distributions with several small classes (less than five elements) and strong asymmetry, as in our case, as pointed out by Martínez-Rodríguez et al. (2011). Therefore, the evaluation of the goodness-of-fit was performed by operating with artificial IT synthetic samples generated via Lerch distribution and numerically reconstructing the null hypothesis of χ² via Monte Carlo simulation. A set of 2,000 replicates of the IT sample was generated by randomly sampling from the theoretical distribution obtained by inference, and the associated p-value was defined as the fraction of the 2,000 replicates for which reconstructed χ_j² values (with j = 1, 2, …, 2,000) are greater than χ_ref², obtained by the classical procedure (Hope 1968). Accordingly, the null hypothesis was accepted if the p-value was greater than 0.05, the value that is usually assumed as the significance level.

For the whole year data set and for all year records, Table 5 reports the comparison of the statistic, D, calculated for the alternative model (Lch) and that calculated for the null models (ID = 2–7). The critical chi-square statistic, at the significant level, p = 0.05, is also reported. The statistic D appears greater than the corresponding critical chi squared, in all the series and in all probability distributions, justifying the use of the three-parameter Lch.

Table 5 also shows that the extended log (hereafter referred to as Ext log) shows the best performance (lowest D values) among the other distributions, whereas the one-parameter Geo shows the worst performance (greatest D values).

As can be noticed from the comparison of the elevation data reported in Table 1 and from Figure 1, the selected stations include a wide range of altitudes (from 113 m to 2,170 m a.s.l.), which is representative of the whole region. The comparison between the average IT, <IT>, reported in Table 6, highlights substantial differences occurring between the two seasons. However, a variety of conditions are present among the stations, within the same season, as is illustrated by the minimum and maximum (shown in bold in Table 6) of the average IT values, <IT>. Table 6 also reports the starting and the ending year of the observations, which accounts for the results obtained by the application of the Pettitt test, the average cluster size, <WS>, calculated according to Equation (2), which also characterizes each station.

In particular, Table 6 shows a minimum <IT> value of 2.7 days in the ‘warm season’ for Lago Toggia station, one of the more northern stations, and a maximum one of 6.57 days in the ‘cold season’ for Cuneo station, one of the more southern ones. As can be observed, the p-values associated with the χ² test are greater than 0.05, for the ‘cold season’, for all the IT data sets, whereas, for the ‘warm season’ and for the all-years record, only 22 and 18 data sets satisfy the χ² test, respectively.

For the 26 considered stations, Figure 5 illustrates the comparison between the observed frequencies and the theoretical probabilities calculated according to the Lch distribution. As can be observed, the dots are close to the line of perfect agreement, especially for high probabilities, showing evidence of the suitability of Lch in describing empirical ITs. In particular, the goodness-of-fit for IT = 1 also emphasized the Lch capability in describing the geometrical distributed WS (Equations (3), (5b) and (11b)).

Figure 6 also illustrates the fitting of the more common polylogarithmic (Poly) and of the logarithmic distribution (Log), showing that with respect to the Lch distribution, they are not able to reproduce the observed IT frequencies.

The cumulative Lch distributions for the whole year data set, are compared with those corresponding to the ‘cold’ and to the ‘warm’ seasons, in order to demonstrate the seasonality effect, in Figure 7, for Lago Toggia and Torino stations. For both stations and with respect to the whole year Lch distributions, the figure shows greater differences in probability for the ‘cold season’ than for the ‘warm season’.

With reference to the whole year data set, in order to synthetize results of the statistical inference via the Lerch distribution, for the all considered stations, Figure 8(a) compares the observed WS frequency distribution with the geometric distribution derived according to Equation (2), with parameter 1 − p₁ ≡ P{IT >1}. Regarding the mentioned series of Lago Toggia and Torino stations, Figure 8(a) shows that the geometric distribution fits the empirical one well. Analogously, Figure 8(b) illustrates a good agreement of the DS empirical distribution with the theoretical polylogarithmic one derived from IT theoretical distribution, according to Equation (3). Thus, Figure 8(a) and 8(b) together confirm that both WS and DS distributions can be derived from IT distribution fitted with the three-parameter Lerch. The good agreement between the geometric WS distribution and the empirical one, illustrated in Figure 8(a), suggests that the assumption that ITs are independent random variables with identical distribution (i.i.d.) is well founded.

The above described ‘wet day sequences’ can represent in a broad sense a unique perturbed period, even though just 1-day DS occur in between the same sequences. Equation (13) could have a certain interest from a hydrological point of view, if considering that 1 no-rainy day in between 2 rainy days does not significantly alter the root zone soil moisture status, mainly for clayey soils.

Regarding the parameters of the Lch distribution, for the whole year data set, the θ values varies in the interval (0.84 ÷ 0.92) and, as expected, higher values are associated with the ‘drier’ stations and lower values with ‘wetter’ ones. Moreover, the range of θ values associated with the ‘warm season’ (0.82 ÷ 0.91) is close to that of the whole year, whereas for the ‘cold season’ the range of θ values (0.90 ÷ 0.95) is narrower, if compared to that of the whole year data set. In Figure 11(a), the θ values are plotted vs. <IT>, for the all-year data set, for the ‘warm season’ and for the ‘cold season’. As can be observed, ‘θ behaves somewhat like a scale parameter’ (Kemp 2006) and can be also related to a physical parameter, the elevation, which is used as the independent variable in Figure 11(b). In particular, θ values are negatively correlated with elevation, with coefficient of determination equal to 0.86, 0.86 and 0.78 for the ‘cold season’, the ‘warm season’ and for the whole year data sets, respectively.

Figure 12 compare parameters θ, s and a, of the ‘warm season’ with those of the ‘cold season’. It is evident that θ values of the ‘warm season’ show a relationship with the corresponding θ values of the ‘cold season’, whereas s and a parameters do not show a clear behaviour, with a similar and narrow range of variability. Furthermore, no relationship of s and a parameters with other physical/climatic variables have been found. In the same figure, the corresponding values of Lch parameters obtained for the Sicily region in a previous work (Agnese et al. 2014) are reported. Similar considerations are valid for the parameters associated with the Sicily region, although the values are different, and their range of variability is wider, compared to that of the Piedmont and Aosta Valley area, likely indicating the different precipitation regime of the two areas. It is interesting that the warm season Sicily's θ, s and a values are in a completely different range. The θ values cold season also show a dfferent range.

To conclude the analysis of the results, Figure 13 shows the average cluster size <WS>, which according to the geometric distribution (see Equation (2)) is equal to 1/(1 − p₁), as a function of the IT harmonic mean, for the ‘cold season’, for the ‘warm season’ and for the whole year data set. A clear relationship is detected, according to the fact that the harmonic mean enhances low IT values. Moreover, the considered seasons slightly influence this relationship.

An application of the described Lch pdf at the regional scale was carried out, on the basis of the good agreement between empirical and theoretical distribution. In particular, this application utilizes the good agreement in the wide range of probability values, including extremes with low-level probability of exceedance. Thus, the application was performed for a fixed long dry spell (DS = 45 days) and for a fixed wet spell (WS = 5 days), the values of which, and their associated spatial distribution, can play an important role in farm and environmental management. In fact, high duration of DS influences crops' growth, whereas the probability of occurrence of long precipitation spells is associated with environmental risks. Results are represented in terms of return period, T (years) and are illustrated in Figure 14(a) and 14(b). The expected return period was evaluated by using the known average IT values calculated for the whole year record (Y), reported in Table 4, according to ⁠. The probability was calculated for each station by Equation (3). The analysis was carried out for the whole year data set, thus separated ‘cold season’ and ‘warm season’ series (Table 4) were not considered. It is interesting to observe that the high elevation areas (compare Figure 1), are generally characterized by lower return period of WS = 5 days and higher return period of DS = 45 days, if compared with the low elevation areas. It is important to highlight that the expected return period for the above-mentioned was determined under the assumption that both WS and DS are time independent and uncorrelated random variables. In this particular case, the extended return period formulation introduced by Fernández & Salas (1999) and then by Volpi et al. (2015), which differentiated the return period as an average occurrence interval or as an average recurrence interval, according to a not time-dependent or a time-dependent condition, respectively, leads to the same results.

Interestingly, for the case of time-dependent random variables, Volpi et al. (2015) also introduced an equivalent return period (ERP), which resembles the classical definition of return period, but in the case of independence. The ERP, which for this study matches the classical return period formulation considered above, is able to control the probability of the random variable under the time-dependent condition, preserving the virtue of the classical return period formulation that is insensitive to the time-dependent condition.

The statistical analysis performed in this work demonstrated the applicability of the discrete three-parameter Lerch probability distribution to the frequency distribution of the inter-arrival time-series in the NW of Italy, characterized by precipitation regimes ranging between that of hilly areas influenced by the Tyrrhenian Sea, to that of the Po Plain and to that of the Alps. In the past, the same distribution has already been successfully tested in the Mediterranean area (Sicily island) (Agnese et al. 2014), which has very different precipitation regimes from those investigated in this study.

The analysis was successfully applied to 70–90 years' long inter-arrival time series derived from 26 daily precipitation stations records and the i.i.d. hypothesis of the IT random variable was then checked. Because of the climatic spatial variability of these regions, and of Sicily, the application of the three-parameter Lerch distribution has demonstrated its wide applicability. The spatial variability of all the three parameters of the distribution, for both ‘cold’ and ‘warm’ seasons, is smaller than in Sicily. In particular, the differences found in the θ parameter, confirm its intrinsic character of scale parameter, whereas both s and a parameters appear to be purely fitting parameters, but with markedly different ranges during the cold season, with respect to Sicily. In this paper, it is shown that the Lerch distribution successfully describes the sequences of rainy days interrupted by 1-day DS, with useful practical applications from a hydrological point of view. Moreover, the equivalence on using a single IT three-parameter Lerch distribution, in place of two distributions, namely, the DS two-parameter polylogarithmic and WS one-parameter geometric ones, is also tested. Finally, an application at the regional scale describing the probability of exceedance of a long wet spell and a long dry spell, has confirmed the value of this distribution in quantifying the wide spatial variability of the precipitation statistics in the considered area.

Research was supported by PRIN 2015 project granted by MIUR (Italian Ministry for University and Research) No. 2015AKR4HX, by the H2020 project ‘ECOPOTENTIAL: Improving Future Ecosystem Benefits Through Earth Observations’, coordinated by CNR-IGG (http://www.ecopotential-project.eu) (it has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 641762), and by the CNR flag project NEXTDATA.