Hazard Assessment under Multivariate Distributional Change-Points: Guidelines and a Flood Case Study

Abstract

One of the ultimate goals of hydrological studies is to assess whether or not the dynamics of the variables of interest are changing. For this purpose, specific statistics are usually adopted: e.g., overall indices, averages, variances, correlations, root-mean-square differences, monthly/annual averages, seasonal patterns, maximum and minimum values, quantiles, trends, etc. In this work, a distributional multivariate approach to the problem is outlined, also accounting for the fact that the variables of interest are often dependent. Here, the Copula Theory, the Failure Probabilities, and suitable non-parametric statistical Change-Point tests are used in order to provide an assessment of the hazard. A hydrological case study is utilized to illustrate the issue and the methodology (viz., assessment of a dam spillway), considering the bivariate dynamics of annual maximum flood peak and volume observed at the Ceppo Morelli dam (located in the Piedmont region, Northern Italy) over a 50-year period. In particular, several problems—often present in hydrological analyses—are debated: namely, (i) the uncertainties due to the presence of heavy tailed random variables, and (ii) the hydrological meaning/interpretation of the results of statistical tests. Furthermore, the suitability of the procedures proposed to fulfill the goals of the study (viz., detecting and interpreting non-stationarity) is discussed. Overall, the main recommendation is that statistical (multivariate) investigations may represent a necessary step, though they may not be sufficient to assess hydrological (environmental) hazards.

1. Introduction

In the current hydrological practice, water engineering works, including dams, levees, detention basins, and sewers, are designed under the hypothesis of stationarity of the random variables at play. This assumption entails the time invariance of the probability distribution of the variables (strong stationarity), which implies that the family of distribution and its parameters are fixed and constant, or the time invariance of the statistical moments (weak stationarity)—see [1,2].

In hydrological literature, the stationarity assumption has mainly been investigated under a univariate framework [3]. For instance, this is the case of streamflow, traditionally considered as the design variable. In some cases, stationarity can be considered as a valuable working assumption, and as a valid approximation of the reality [4,5]. In turn, possible departures from stationarity are judged to slightly affect the results [3].

However, in addition to possible environmental factors, in [6], it is claimed that many human activities, like urbanization, deforestation, change of agricultural practice, anthropogenic emissions, and river engineering constructions, may contribute to introduce significant changes in the components of the hydrological cycle, including streamflow. In particular, it is argued that these activities might have compromised the assumption of stationarity in hydrology.

Should the stationarity assumption not be valid, then the probability distribution of the variables at play would be time-varying, with a possible presence of trends or abrupt shifts, entailing changes in hazard evaluation and assessment. In such a case, a revision of the design criteria may be necessary, in order to avoid underestimation or overestimation of design variables, with a consequent inadequacy of the designed works on the one hand, or an increase of costs of the structure on the other hand [7,8]. Thus, testing the stationarity may represent a fundamental step.

Several techniques are used to test stationarity in hydrologic time series: e.g., trend analysis, spectral analysis, multi-resolution methods, etc.—see [3] for a review. In this work, we outline a further possible approach of distributional nature, involving the multivariate probability distribution of the variables at play, by using statistical procedures recently introduced in literature—see the discussion after Equation (1) below: indeed, according to [9], the full distributional (probabilistic) features of hydrologic time series are rarely explored.

Furthermore, we focus the attention on the possible consequences of violating the stationarity assumptions from a hazard assessment perspective. Specifically, Change-Point identification and hazard management are carried out by exploiting the Theory of Copulas, and the notion of Failure Probability, respectively. Thanks to the separability of the dependence structure and the marginal distributions (via Sklar’s Theorem—see Equation (1) later), the copula approach may simplify the assessment of the impact of potential violations of the stationarity assumption on the joint hazard, as well as on the marginal profiles. As an illustration, maximum annual flood peak and volume are used as the variables of interest for the assessment of a dam spillway, and related problems typical of hydrological analyses are discussed (e.g., how heavy tailed variables may affect the analyses, and the interpretation of the statistical tests).

2. Materials

The hydrological data investigated in the following are collected at the Ceppo Morelli dam (Northern Italy), and are the same ones considered in [10,11,12,13], to which the reader is referred. Maximum annual flood peaks Q and volumes V are identified and selected for 49 years, from 1937 to 1994 (some years are missing—see Figure 1). In turn, the case study used to illustrate the general-purpose methodology outlined in this paper—which can be adopted for investigating any environmental hazard—concerns flood hazards involving the traditional variables peak and volume.

Interestingly enough, this database represents an exceptional case study as compared to traditional hydrological ones present in literature: in fact, in hydrological practice, usually the maximum annual flood peak is first selected as a design variable, and then the volume associated with the same flood event is considered, even if it does not represent the annual maximum volume—see, e.g., [14,15]. Instead, in our case, 48 out of 49 of the occurrence dates of the Qs and the Vs are the same (viz., they took place during the same flood episode). In turn, on the one hand, the use of a Block Maxima approach for the Extreme Value analysis is well justified (here, the blocks are the years), and, on the other hand, the “hydrological consistency” of the database is preserved.

3. Methods

A convenient way to deal with multivariate phenomena, where the variables at play are generally non-independent, is to use Copulas [16,17,18]. Since the publication of [19], a number of papers in hydrology, as well as in other environmental areas, have shown the theoretical and practical advantages of using a copula approach, and support its usage. For an overview concerning different ways of quantifying the hazard of compound events, see, among others, [20,21,22,23,24]. In particular, concerning selection/estimation/test statistical procedures for copulas, the interested reader may refer to [25,26,27,28,29,30], and references therein. Note that valuable software for working with copulas, developed for the R package [31], is freely available online [27,32]. The results presented later are obtained using the methodologies outlined in the cited works, to which the reader is referred. In particular, in  the following, the same notation used in [11,12,13] is adopted.

Let $\mathbf{X}=(X_1,\dots ,X_d)$ be the random vector describing the phenomenon under investigation, with univariate marginals $F_i$s, i.e., $F_i\left(x_i\right)=\mathbf{P}(X_i\le x_i)$. According to Sklar’s Theorem [33], the joint distribution function $\mathbf{F}$ of $\mathbf{X}$ can be written as

$$\mathbf{F}\left(\mathbf{x}\right)=\mathbf{P}(X_1\le x_1,\dots ,X_d\le x_d)=\mathbf{C}(F_1\left(x_1\right),\dots ,F_d\left(x_d\right)),$$

(1)

where $\mathbf{C}$ is the Copula (viz., the Dependence Structure) of the variables $X_i$s, assumed to be continuous to ensure the unicity of $\mathbf{C}$ (as is usually the case in hydrological applications).

As recently discussed in [34], the representation provided by Sklar’s Theorem provides a valuable theoretical tool in order to assess the presence of possible changes of the distributions of hydrological/climate variables (called Change-Points in Statistics)—see also [35,36,37,38] for previous different approaches. In particular, the non-stationarity of $\mathbf{F}$ may be due to

• changes of any of the marginal distributions $F_i$s, or
• changes of the copula $\mathbf{C}$, or
• both of the previous cases.

In turn, non-parametric Change-Point statistical tests recently outlined in [39,40,41,42] (as well as the corresponding software [43]) can be used to check whether the (multivariate) Null distributional assumption “${\mathcal{H}}_0$: $\mathbf{F}$ has no Change-Points” (viz., $\mathbf{F}$ does not change with time) should be rejected or not, by investigating whether either the marginals, or the copula, or both, show Change-Points—see, later, Section 4.

Furthermore, let $\mathbf{I}$ denote the unit interval $[0,1]$, and let ${\mathcal{L}}_t$ be the level set of $\mathbf{F}$ at $t\in \mathbf{I}$: practically, ${\mathcal{L}}_t$ is the set of points in the d-dimensional Euclidean space ${\mathbf{R}}^d$ such that $\mathbf{F}\left(\mathbf{x}\right)=t$—sometimes ${\mathcal{L}}_t$ is also referred to as a “critical layer” of level t [11]. As will be made clear below, ${\mathcal{L}}_t$ plays the role as of a (critical) multivariate threshold, with dimension $d-1$.

Due to the assumption that the $F_i$s are strictly increasing—as is the case of the distributions traditionally used in hydrological practice, the  Probability Integral Transform (hereinafter, PIT) viz., the relations

$$\mathbf{U}=(U_1,\dots ,U_d)=(F_1\left(X_1\right),\dots ,F_d\left(X_d\right))={\mathcal{T}}_{\mathbf{F}}\left(\mathbf{X}\right),$$

(2)

and

$$\mathbf{X}=(X_1,\dots ,X_d)=(F_1^{-1}\left(U_1\right),\dots ,F_d^{-1}\left(U_d\right))={\mathcal{T}}_{\mathbf{F}}^{-1}\left(\mathbf{U}\right),$$

(3)

are one-to-one, where the $U_i=F_i\left(X_i\right)$s are Uniform random variables on $\mathbf{I}$. These formulas map the vector $\mathbf{U}$ living in ${\mathbf{I}}^d$ onto the vector $\mathbf{X}$ living in the d-dimensional Euclidean space ${\mathbf{R}}^d$ (and vice-versa)—see [16,17,44]. Since copulas are invariant for strictly increasing transformations of the variables at play ([16] Theorem 2.4.3), $\mathbf{U}$ and $\mathbf{X}$ share the same copula. Note that the PIT uniquely maps the probabilities of the events in ${\mathbf{I}}^d$ (as induced by the copula $\mathbf{C}$) onto ${\mathbf{R}}^d$ (as induced by $\mathbf{F}=\mathbf{C}(F_1,\dots ,F_d)$), and vice-versa. The role played by the univariate marginals is only to geometrically re-map such probabilities from ${\mathbf{I}}^d$ onto suitable regions in the Euclidean space ${\mathbf{R}}^d$ (and vice-versa), without affecting them. By the same token, also ${\mathcal{L}}_t$ is uniquely mapped from ${\mathbf{R}}^d$ onto ${\mathbf{I}}^d$ (and vice-versa), thus becoming a level set of $\mathbf{C}$.

3.1. Hazard Scenarios

The notion of Hazard Scenario introduced in [13] is fundamental and is briefly recalled below.

Definition 1.

Let $\mathbf{X}$ model the phenomenon of interest. A Hazard Scenario (hereinafter,HS) of level $\alpha \in (0,1)$ is any Upper Set $\mathcal{S}={\mathcal{S}}_{\alpha }\subseteq {\mathbf{R}}^d$ such that the following relation holds:

$$\mathbf{P}(\mathbf{X}\in \mathcal{S})=\alpha .$$

(4)

A Hazard Scenario is simply a set containing occurrences $\mathbf{x}$s that may damage a structure. By the very definition of Upper Set, if $\mathbf{x}\in \mathcal{S}$ and $\mathbf{y}\ge \mathbf{x}$ component-wise, then also $\mathbf{y}$ could be considered as dangerous, since it “exceeds” $\mathbf{x}$.

Given a realization $\mathbf{x}\in {\mathbf{R}}^d$, there exist several ways [13] to associate $\mathbf{x}$ with a suitable HS ${\mathcal{S}}_{\mathbf{x}}\subseteq {\mathbf{R}}^d$. Note that, via Equations (2) and (3), there exists a one-to-one correspondence between ${\mathcal{S}}_{\mathbf{x}}$ and a specific region ${\mathcal{S}}_{\mathbf{u}}\subseteq {\mathbf{I}}^d$. In turn, the knowledge of the copula at play may suffice to calculate the level of ${\mathcal{S}}_{\mathbf{u}}$, and hence of ${\mathcal{S}}_{\mathbf{x}}$ (since they share the same probability).

In general (see, e.g., [13,23]), the choice of the HS should depend upon two criteria: (i) the type of dangerous events, and (ii) their probabilities of occurrence. The  approach of interest here is the “OR” one [13], as defined below: in fact, it is sufficient that EITHER Q, OR V, OR both, be large in order to affect/damage the dam (spillway) of interest.

Definition 2.

Given $\mathbf{x}\in {\mathbf{R}}^d$, the associated d-dimensional “OR” Hazard Scenario ${\mathcal{S}}_{\mathbf{x}}^{\vee }$ is given by the region

$${\mathcal{S}}_{\mathbf{x}}^{\vee }=\bigcup _{i=1}^d(\mathbf{R}\times \cdots \times (x_i,+\infty )\times \cdots \times \mathbf{R}),$$

(5)

or, equivalently,

$${\mathcal{S}}_{\mathbf{x}}^{\vee }=\{\mathbf{z}\in {\mathbf{R}}_{+}^d:{\Psi }_{\mathbf{x}}\left(\mathbf{z}\right)=max\left\{\frac{z_1}{x_1},\dots ,\frac{z_d}{x_d}\right\}\ge 1\},$$

(6)

where the function Ψ (introduced in [45]) models and rules the occurrence of dangerous events (viz., when $\Psi \ge 1$). The corresponding level α is given by

$${\alpha }_{\mathbf{x}}^{\vee }=\mathbf{P}\left(\mathbf{X}\in {\mathcal{S}}_{\mathbf{x}}^{\vee }\right)=1-\mathbf{C}(F_1\left(x_1\right),\dots ,F_d\left(x_d\right)).$$

(7)

For the realization of the event $\{\mathbf{X}\in {\mathcal{S}}_{\mathbf{x}}^{\vee }\},$ it is sufficient that one (or more) of the variables $X_i$s exceed the corresponding critical threshold $x_i$ (which, usually, is specified by the Regulation, or is inferred by the problem at hand).

3.2. The Failure Probability Approach

A consistent way to assess the hydrological hazard is to compute the corresponding Failure Probability (hereinafter, FP)—for a recent discussion, see [23,24,46]. In this section, the Failure Probability is calculated for the “OR” Hazard Scenario of interest here.

Let ${\mathbf{X}}_1,\dots ,{\mathbf{X}}_T$ be the d-variate vectors describing the phenomenon of interest at times $1,\dots ,T$, where $T>0$ is an arbitrary design life time for a given structure: without loss of generality, here T is measured in years. As an example, think of a series of floods $(Q_t,V_t)$s, where $Q_t$ and $V_t$ represent, respectively, the annual maximum flood peak and volume observed in the t-th year.

Now, let $\tilde{\mathbf{x}}$ denote a bivariate critical design threshold (e.g., the one prescribed by the Regulation). Then, each occurrence ${\mathbf{x}}_i$s can be associated with a precise “OR” HS, given by the region in the $(Q,V)$ plane where either the peak, or the volume, or both exceed the corresponding component of $\tilde{\mathbf{x}}$. In turn, a set ${\mathcal{S}}_1,\dots ,{\mathcal{S}}_T$ of “OR” HSs is generated, and each ${\mathcal{S}}_i$ has level ${\alpha }_i$ given by Equation (7). From a practical point of view, the  floods occurring in the complementary regions ${\mathcal{S}}_1^c,\dots ,{\mathcal{S}}_T^c$ of the HS’s ${\mathcal{S}}_i$s could be labelled as “safe”.

According to ([47] Chapter 12) and ([48] Chapter 9), the  Failure Probability $p_T$ can be computed as

$$p_T=1-\mathbf{P}\left({\mathbf{X}}_1\in {\mathcal{S}}_1^c,\dots ,{\mathbf{X}}_T\in {\mathcal{S}}_T^c\right).$$

(8)

In turn, considering the case of independent and identically distributed ${\mathbf{X}}_i$s, according to Equations (7) and (8), the FP corresponding to the “OR” HS (the one of interest here), for events sharing a common multivariate critical threshold $\tilde{\mathbf{x}}$, is given by:

$$p_T^{\vee }=1-{\left(\mathbf{C}(F_1\left({\tilde{x}}_1\right),\dots ,F_d\left({\tilde{x}}_d\right))\right)}^T,$$

(9)

or, equivalently,

$$p_T^{\vee }=1-\mathbf{P}\left({\Psi }_{\tilde{\mathbf{x}}}\left(\mathbf{X}\right)\le 1\right)=1-\mathbf{P}\left(max\left\{\frac{X_1}{{\tilde{x}}_1},\dots ,\frac{X_d}{{\tilde{x}}_d}\right\}\le 1\right).$$

(10)

In order to provide valuable information for the estimate of, e.g., suitable multivariate design quantiles, further work may be required. In fact, all the infinite realizations lying on the critical layer ${\mathcal{L}}_t$ are associated with the same value of the Failure Probability, since $p_T^{\vee }$ is constant over the level set ${\mathcal{L}}_t$ [13]. This may leave undetermined the assignment of a specific design occurrence, once T and $p_T$ have been chosen. In turn, valuable design values, associated with given Life Times and Failure Probabilities, could be calculated via suitable strategies, e.g., mimicking those outlined in [11,30,49,50,51]. A practical example is given below in Section 4.

4. Results and Discussion

A thorough statistical analysis indicates that Generalized Extreme Value (hereinafter, GEV) marginal distributions adequately fit the observations of both Q and V (see

Table 1. Maximum Likelihood estimates of the parameters of the GEV distributions fitting the variables Q (in m${}^3$/s) and V (in 10${}^6$ m${}^3$), either considering all the data, or only those collected, respectively, before and after the Change-Point year (1971)—see text. Also shown are estimated standard errors and approximate Monte Carlo Goodness-of-Fit test p-Values (based on Kolmogorov–Smirnov statistics).

Variable	Shape	Scale	Position	p-Value
	All data
Q	0.37	36.21	59.36	77%
s.e.	0.11	5.04	5.71
V	0.61	1.52	1.72	91%
s.e.	0.13	0.25	0.24
	Before Change-Point
Q	0.02	24.46	50.05	87%
s.e.	0.11	3.79	5.36
V	0.12	0.95	1.39	99%
s.e.	0.16	0.16	0.21
	After Change-Point
Q	0.71	42.96	71.74	98%
s.e.	0.30	11.75	10.92
V	1.07	2.03	2.20	92%
s.e.	0.34	0.67	0.50

): this is coherent with the Extreme Value approach adopted, since the variables are annual maxima. The Monte Carlo p-Values of the Kolmogorov–Smirnov Goodness-of-Fit tests are larger than 10%, supporting the validity of the assumed model from a statistical point of view.

Concerning the bivariate dependence structure of the $(Q,V)$s, the two variables turn out to be non-independent: in fact, the estimates of the Kendall’s ${\tau }_{QV}\approx 0.66$ and the Spearman’s ${\rho }_{QV}\approx 0.80$ are both statistically significantly positive at a 5% level. Here, a survival-Clayton 2-copula ${\mathbf{C}}_{QV}$ is selected among dozens of competing bivariate models, and provides a valuable fit (as indicated by the results of the Goodness-of-Fit (GoF) tests shown in

Table 2. Maximum Likelihood estimates of the survival-Clayton 2-copula parameter $\theta$ fitting the pairs $(Q,V)$s, either considering all the data, or only those collected, respectively, before and after the Change-Point year (1971)—see text. Also shown are estimated standard errors and approximate p-Values (via a Multiplier Method) of the Cramér–von Mises Goodness-of-Fit test for Copulas based on the copula empirical process [32].

	All Data	Before Change-Point	After Change-Point
$\theta$	4.33	1.53	11.69
s.e.	1.37	0.72	5.31
p-Value	9%	47%	44%

). Once the marginals and the copula have been fixed, the joint distribution ${\mathbf{F}}_{QV}$ can be computed via Sklar’s Theorem: Figure 2 shows selected isolines of ${\mathbf{F}}_{QV}$.

For the sake of illustration, critical bivariate design thresholds $\tilde{\mathbf{x}}$s, to be used in Equations (5), (7), and (9) for the calculation of the quantities of interest, are fixed as follows—clearly, other criteria could be adopted, and other values could be chosen. Here, ${\tilde{x}}_1$ and ${\tilde{x}}_2$ are assigned the empirical 90%, 95%, and 99% quantiles of, respectively, Q and V (see

Table 3. The design pairs $\tilde{\mathbf{x}}$s—see text: the units are m${}^3$/s for Q, and 10${}^6·$m${}^3$ for V. Also shown are the corresponding estimates of the “OR” HS levels ${\alpha }^{\vee }$s.

	Quantile (%)	Q	V	${\mathit{\alpha }}^{\vee }$
${\tilde{\mathbf{x}}}_1$	90%	197	14	0.0880
${\tilde{\mathbf{x}}}_2$	95%	282	17	0.0453
${\tilde{\mathbf{x}}}_3$	99%	439	31	0.0172

). Figure 3 shows the three design pairs $\tilde{\mathbf{x}}$s, as well as the corresponding critical layers and the “OR” Hazard Scenarios of interest. The corresponding approximate levels ${\alpha }_i^{\vee }=1-{\mathbf{C}}_{QV}(F_Q\left({\tilde{x}}_{i,1}\right),F_V\left({\tilde{x}}_{i,2}\right))$, with $i=1,2,3$, can be calculated via Equation (7), and are reported in Table 3. As expected, (component-wise) “larger” $\tilde{\mathbf{x}}$s yield smaller scenario levels ${\alpha }^{\vee }$s.

As an example, in this work, a life time T, varying from one to 50 years, is chosen. In the first instance, here it is assumed that the occurrences are independent and identically distributed (i.i.d.): while physical and statistical reasons suggest accepting the independence of annual maxima, by contrast the stationarity of their joint distribution might be questionable (e.g., due to a changing climate or human activities), and will be discussed later.

The results under the i.i.d. assumption are shown in Figure 4(top): here, Monte Carlo Confidence Bands (at a 90% level) for the Failure Probabilities $p_T^{\vee }$s associated with the design pairs $\tilde{\mathbf{x}}$s are plotted. In  all cases, the width of the bands is very large: the explanation is as follows. The  shape parameters of the fitted GEV marginals are close to (or larger than) 1/2—see Table 1: interestingly enough, the parameters have also been estimated via the L-moments method, yielding comparable results (not shown). In turn, the existence of the second order moments (i.e., the variances) may be questionable, and large fluctuations during the simulations have to be expected, which may adversely affect the Monte Carlo procedure and yield very wide bands. Unfortunately, situations like these are common in hydrological practice, but little can be done to reduce the uncertainties shown in the plots: it is a problem intrinsic to the very structure of the available data.

Concerning the computation of suitable multivariate design quantiles, here a Most Likely strategy is used [11,30]. Practically, the level set ${\mathcal{L}}_t$ associated with selected values of T and $p_T$ is considered, and the occurrence ${\delta }^{*}=(Q^{*},V^{*}),$ which maximizes the density ${\mathbf{f}}_{QV}$ over ${\mathcal{L}}_t$, is taken as a design pair. Figure 5(top) shows the design pairs corresponding to a lifetime $T=50$ years and Failure Probabilities of order, respectively, 10%, 5%, and 1%, as well as Monte Carlo Confidence Intervals at a 90% level: as already discussed above, large unavoidable uncertainties are present.

As a second instance, it may be of interest to study how the previous results might be affected by a change of the joint distribution ${\mathbf{F}}_{QV}$ of the phenomenon under investigation (for a similar hydrological case study see, e.g., [9,52]). Specifically, a distributional change of the behavior of the pair $(Q,V)$ may be due to:

• a change of the univariate distribution $F_Q$ (respectively, $F_V$), or
• a change of the copula ${\mathbf{C}}_{QV}$ associated with $(Q,V)$, or
• both the previous instances.

Here, possible changes are detected according to different tests recently introduced in literature—see [39,40,41,42], and  implemented in the R package npcp [43].

In the present case, apparently, the overall stochastic behavior shown in Figure 1 might have changed around the year 1971: see below, and Table 1 and Table 2. More precisely, according to the approximate p-Value of the test (about 4%), the Null hypothesis that the joint distribution ${\mathbf{F}}_{QV}$ has not changed could not be rejected at the standard 1% level, but should be rejected at the 5% and 10% ones. Actually, $F_Q$ might also have changed around 1971 (the p-Value is about 3%), as well as ${\mathbf{C}}_{QV}$ around 1970 (the p-Value is about 2.5%), and $F_V$ around 1975 (the p-Value is about 2.8%). In turn, given the intrinsic uncertainties of the test [39], the critical year associated with the Change-Point is reasonably and practically assumed to be 1971 (also considering that, in hydrology, Q is traditionally indicated as the regulation variable). Apparently, this is consistent with the claims of [53], stating that “… the 1970s are known as a period of major climate and environmental changes observed in several proxies and several fields at the global scale, as attested to in many recent studies [54,55,56,57,58,59].”

Whether or not the dynamics of the phenomenon has really changed (concerning this latter issue, instructive is the recent work by [4]), it is however interesting to study the possible effects on the strategies of hazard assessment, e.g., in terms of the Failure Probabilities discussed in this work. For this purpose, both the marginals and the copula of the subsets of observations before and after 1971 were computed: the results are reported in Table 1 and Table 2, and illustrated in Figure 2 and Figure 6. In turn, both the GEV univariate laws and the survival-Clayton 2-copula still provide valuable models for these data sets, but the corresponding parameters have changed. More specifically, the observations can be fitted by distributions and copulas showing two distinct regimes (before and after 1971), and modeled by the same univariate and bivariate functions with different parameters. Note that no further Change-Points are detected in each of the two sub-periods before and after 1971.

The behavior of the FPs, assuming that the joint probability distribution has changed in 1971, is plotted in Figure 4(bottom): in general, the failure probabilities assuming stationarity are smaller than the ones after 1971, while the uncertainties are again large. Apparently, there is a remodeling of the dynamics of the temporal series of floods: actually, both univariate and multivariate parameters change. In particular, an increase of the statistical dependence between flood peak and flood volume may be likely. These results may indicate an intensification of the flood events in the period 1971–1994, both in the dependence structure and in the marginals, and consequently an increase of the corresponding failure probabilities. Similar conclusions can be drawn concerning the computation of suitable design pairs, as shown in Figure 5(bottom): again, for the same reasons mentioned above, large uncertainties are present.

5. Conclusions

The possible presence of a distributional Change-Point in a hydrologic bivariate dataset of maximum annual flood peaks and volumes has been investigated, in order to evaluate the consequences of possible violations of the stationarity assumption from a hazard assessment perspective. In particular, the assessment of the environmental hazard has been carried via the theory of Copulas and the Failure Probabilities. The impact is evaluated adopting a risk-manager perspective.

The investigation carried out in this work makes it evident that it may be important to test the stationarity assumption, viz. to check whether (or not) the stochastic behavior of the variables at play might have changed with time, in terms of a change of the families of marginal probability distributions, and/or a change of the associated copula family, and/or a change of the corresponding parameters’ values.

Apparently, considering the case study investigated here, the statistical analyses indicate that a distributional Change-Point might be present around 1971: actually, the estimates of the parameters of both the marginals and the copula are different in a time period before and after 1971 (whereas the families of the distributions might not have changed). In turn, a due analysis (in terms of Failure Probabilities) of the effects of the possible presence of a Change-Points on the hazard assessment has been carried out.

As a conclusion, under the assumption that the hydrological regime has really changed, the corresponding increase of the FPs yields a strengthening of the estimated threatening: clearly, neglecting such a dynamics might provide an unrealistic hazard assessment. However, such a diagnosis is essentially of a statistical nature: in our opinion, further physical considerations are needed in order to correctly guide the decisions of the Water Managers, which only a proper and correct engineering practice may provide. In our opinion, further analyses are needed in order to understand whether a similar behavior is present also in other neighboring sites, and to check whether analogous conclusions could be drawn also for, e.g., precipitation, droughts, etc. (see, e.g., the scenarios’ procedures outlined in [60,61], as well as [62] for recent multivariate characterizations). Actually, combined together, such pieces of information could provide valuable evidence concerning possible climate changes, and supply indications about plausible future hydrologic scenarios. Finally, the presence of heavy tailed marginals (a common situation in hydrological time series) further complicates the analyses by introducing additional uncertainties, as discussed in the paper.

As a general recommendation, it is important to stress that the statistical outcomes should always be “validated/supported/confirmed” via additional practical investigations, such as a regional analysis of other relevant hydrological variables like discharge, precipitation, etc. In terms of flood hazard management, the results presented here suggest that, for checking the adequacy and the functionality of existing water works, a statistical (multivariate) survey represents a necessary step, though it may not be sufficient.