Flood frequency analyses are usually based on the assumption of stationarity, which might be unrealistic if changes in climate, land uses or urbanisation impact the study catchment. Moreover, most non-stationarity studies only focus on peak flows, ignoring other flood characteristics. In this study, the potential effect of increasing urbanisation on the bivariate relationship of peak flows and volumes is investigated in a case study in the northwest of England, consisting of an increasingly urbanised catchment and a nearby hydrologically and climatologically similar unchanged rural (control) catchment. The study is performed via Kendall's tau and copulas. Temporal trends are studied visually and by formal tests, considering variables individually and jointly. Bivariate joint return period curves associated with consecutive time periods are compared to understand the joint implications of such bivariate trends. Although no significant bivariate trends were detected, hydrologically relevant trends were found in the urbanised catchment.

## INTRODUCTION

Accurate design flood estimates with specified return periods are necessary for designing and operating hydraulic structures, such as culverts and dams. Traditionally, flood frequency analyses are based on an assumption of stationarity, i.e., assuming that the flood generating processes remain unchanged over time (e.g., Stedinger *et al.* 1993; Goel *et al.* 1998; Yue *et al.* 1999; Shiau *et al.* 2006). It has long been recognised by hydrologists that stationarity is, at best, a simplified working assumption when changes in urbanisation, land uses or climate are involved in the problem under analysis (e.g., Benkhaled *et al.* 2014), as such impacts can affect the behaviour of hydrological variables, e.g., leading to changes in flood characteristics.

In cases where these effects are considered important, a non-stationarity approach should be applied. Mathematical implementations of non-stationarity into flood frequency models, such as using a non-stationary generalised extreme value (GEV) or Pearson type 3 (PE3) distributions are relatively straightforward. However, assessing what is the correct model structure that best describes the impact of changing drivers on the characteristics of the flood series as well as giving realistic predictions of future impacts is far more difficult (Stedinger & Griffis 2011).

The effects of urbanisation on the characteristics of flood runoff have been the subject of several scientific investigations, and it is generally recognised that urbanisation will result in increased runoff volumes (decreased infiltration rates) and reduced catchment lag-times (e.g., Rose & Peters 2001; Shuster *et al.* 2005; Kjeldsen 2009), thus potentially being a significant cause of non-stationarity in flood series impacting both peak flow values and runoff volumes (e.g., Sheng & Wilson 2009).

Much attention has been given to studying trends in peak flow values as a function of time (Petrow & Merz 2009; Wilson *et al.* 2010; Mediero *et al.* 2014), with some exceptions such as the approach presented by Bender *et al.* (2014) for analysing bivariate non-stationary in an uncommonly long peak-volume flood data record (191 years). Indeed, as it is expected that increased urbanisation will lead to changes in other flood characteristics other than the peak flow, in particular the flood volume, a study of the potential changes in multiple variables is necessary to better understand the changes that could affect flood risk under increasing urbanisation. The present study aims to introduce and discuss a simple and general framework to investigate changes of multivariate flood characteristics under increasing urbanisation. This is performed by studying the univariate and bivariate properties of peak (*Q*) and volume (*V*) in two paired case catchments located in the northwest of England. Changes in the univariate peak flow values for these two catchments have already been assessed in Prosdocimi *et al.* (2015), who found that the increase of urbanisation in Catchment 70005 was connected to increasing trends, especially for summer flows. The present work wishes to complement such a study by proving the conceptual framework needed to investigate the bivariate behaviour of flow peaks and flood volume. The analysis of this case study is especially relevant due to the available information about urbanisation levels and flow records of high quality, something not easily found. Nevertheless, the available hydrological record is relatively short and the results presented in this work should be considered as preliminary and taken with caution. The two catchments are hydrologically and climatologically similar, except for the increasing urbanisation levels which affect one of them but not the other. The differences in the changes in the flood characteristics in the two catchments can be imputed to the changes in the land cover of the urbanised catchment. Changes in the association of the bivariate distribution of (*Q,V*) are investigated by both the Kendall's tau (), a rank-based dependence measure, and copulas for bivariate design flood analysis. Copulas (e.g., Joe 1997), which have found several applications in multivariate hydrological analysis (e.g., De Michele *et al.* 2005; Renard & Lang 2007; Klein *et al.* 2010; Ganguli & Reddy 2013; Requena *et al.* 2016), allow obtaining the multivariate joint distribution of multiple random variables by characterising the relation of dependence among them, incorporating the corresponding univariate marginal distributions that can belong to different families.

## CASE STUDY AND DATA EXTRACTION

The rural catchment drained by the River Conder is a nearby hydrologically and climatologically similar catchment, where flow data are recorded near Galgate (gauging station numbered 72014). This is a predominantly rural catchment, which has experienced little change in the study period. Hereafter the catchments are referred to with their gauging station number: the urbanised catchment corresponds to Catchment 70005 and the rural one to Catchment 72014.

*et al.*2015), and the catchments used in this study were the best match given the paucity of long, high quality flow records. In the present study, catchments have similar geographical conditions, being nearby and entailing a similar gauging station elevation. As well, flood events are of the synoptic type for both catchments (Mediero

*et al.*2015). They also entail a similar annual precipitation (SAAR in Table 1), as well as similar oscillations for different quantiles of the catchment average daily rainfall series (Figure 3). Hence, a similar climate is considered.

. | Catchment . | |||||
---|---|---|---|---|---|---|

Descriptors . | 70005 (urban) . | 72014 (rural) . | ||||

AREA [km^{2}] | 54.5 | 28.99 | ||||

BFIHOST | 0.473 | 0.443 | ||||

FARL | 0.964 | 0.975 | ||||

FPEXT [%] | 0.14 | 0.08 | ||||

SAAR [mm] | 1021 | 1183 | ||||

Statistics | Season | Season | ||||

Winter | Summer | Annual | Winter | Summer | Annual | |

Mean Q | 22.18 | 17.31 | 23.62 | 16.06 | 11.42 | 17.43 |

Median Q | 22.78 | 15.85 | 22.95 | 13.85 | 9.64 | 16.45 |

25th quantile Q | 18.03 | 9.73 | 19.78 | 11.78 | 8.02 | 12.95 |

75th quantile Q | 24.75 | 22.18 | 25.97 | 21.98 | 15.45 | 22.65 |

Mean V | 0.57 | 0.32 | 0.56 | 0.28 | 0.14 | 0.28 |

Median V | 0.50 | 0.22 | 0.50 | 0.26 | 0.13 | 0.24 |

25th quantile V | 0.40 | 0.15 | 0.38 | 0.16 | 0.08 | 0.17 |

75th quantile V | 0.63 | 0.47 | 0.66 | 0.38 | 0.20 | 0.34 |

. | Catchment . | |||||
---|---|---|---|---|---|---|

Descriptors . | 70005 (urban) . | 72014 (rural) . | ||||

AREA [km^{2}] | 54.5 | 28.99 | ||||

BFIHOST | 0.473 | 0.443 | ||||

FARL | 0.964 | 0.975 | ||||

FPEXT [%] | 0.14 | 0.08 | ||||

SAAR [mm] | 1021 | 1183 | ||||

Statistics | Season | Season | ||||

Winter | Summer | Annual | Winter | Summer | Annual | |

Mean Q | 22.18 | 17.31 | 23.62 | 16.06 | 11.42 | 17.43 |

Median Q | 22.78 | 15.85 | 22.95 | 13.85 | 9.64 | 16.45 |

25th quantile Q | 18.03 | 9.73 | 19.78 | 11.78 | 8.02 | 12.95 |

75th quantile Q | 24.75 | 22.18 | 25.97 | 21.98 | 15.45 | 22.65 |

Mean V | 0.57 | 0.32 | 0.56 | 0.28 | 0.14 | 0.28 |

Median V | 0.50 | 0.22 | 0.50 | 0.26 | 0.13 | 0.24 |

25th quantile V | 0.40 | 0.15 | 0.38 | 0.16 | 0.08 | 0.17 |

75th quantile V | 0.63 | 0.47 | 0.66 | 0.38 | 0.20 | 0.34 |

The water year in the UK runs from October to September: events occurring between October and March (included) are classified as winter events, while the period April to September constitutes the summer months. The water year 1988–1989 was removed from the study period, as no summer events were available for the gauging station 72014 in this year.

The period 1976–2008 was divided into two equally sized time windows (e.g., Shastri *et al.* 2015), representative of periods of low and high urbanisation levels for Catchment 70005, respectively. Note that by using two equally sized time periods the uncertainty in the estimates for the two time windows can be assumed to be of a similar scale. Only two time windows were considered because of the relatively short common data, although if longer data records were available, the procedure could be applied considering a greater number of time windows. The first time window (named as W1) runs from 1976 to 1992, while the second period (W2) runs from 1993 to 2008. Therefore, the data series (ordered in time), with and the total number of water years, is also divided into and , with . Here, represent the first pairs of and the last pairs. To simplify formulas, hereafter the pairs are presented as with . Depending on the time period considered, makes reference to , or , with *n* the corresponding data length in each case.

*V*, adopted in this study is defined as the part of the event hydrograph above a threshold set at 40% of the flood peak

*Q*, i.e.,

*V*is considered as the volume associated with the upper 60% of the flood event. By considering only flow above a relatively high threshold, the event volumes were not unduly influenced by post-peak small amounts of rainfall causing the flow to increase part way down the recession curve. The 40% threshold used in this study was found to be sufficiently high to remove the nuisance effects caused by secondary rainfall inputs while maintaining the generality of the results. Also, note that volumes extracted using different thresholds were found to be correlated. In this regard, for instance, Karmakar & Simonovic (2007) reported a highly significant correlation between peak and volume regardless of the discharge threshold level. Note that a more detailed analysis of each individual event would not result in more informative results, but rather lead to a less transparent analysis. As an example, the identification of the pair associated with a given water year

*i*is shown in Figure 4.

For both catchments, the annual, summer and winter maxima of the instantaneous peak flow value are identified and the corresponding volumes are extracted. Seasonal maxima are also investigated to better understand whether the changes seen in the annual series are driven by changes in a specific type of events. The autocorrelation for the annual and seasonal series has been plotted and tested, indicating that the standard iid assumption is verified (not shown). Types of floods were not analysed because floods in the UK mainly belong to the synoptic type (Mediero *et al.* 2015). In summary, two catchments are studied, 70005 (urban) and 72014 (rural), two event characteristics are considered (*Q* and *V*), three maximum series of events are extracted (winter, summer and annual maximum flow events), and three time periods are considered: 1976–1992 (W1), 1993–2008 (W2) and 1976–2008 (whole series).

## METHODOLOGY AND RESULTS

The investigation of potential changes in the *Q-V* relationship proceeds as follows: at first, trends in the univariate series for *Q* and *V* separately and in the association between the two variables are studied. Then, a non-parametric procedure is used to assess the statistical significance of the observed trends. Finally, changes in the bivariate return period curves computed via copulas are investigated. Results for the two catchments under study are shown directly after the description of the methodological steps.

### Analysis of univariate flood trends in Q and V series

The first step of the methodology consists in the analysis of univariate temporal trends in the flood series, using visual inspection and the widely used Mann–Kendall test (e.g., Villarini *et al.* 2009; Coch & Mediero 2015), a non-parametric test based on Kendall's *τ*. The Mann–Kendall test (Kendall 1975) is used to assess the null hypothesis of no association between two variables, and the presence of significant temporal trends can be assessed by taking time as one of the variables. The statistical significance is assessed with a two-sided test at a 95% confidence level.

*Q*and

*V*annual and seasonal series for Catchment 70005 and Catchment 72014 (see Table 1 for descriptive statistics). A smaller variation in time of

*Q*and

*V*can be seen in Catchment 70005 in comparison to the characteristics of the flood events recorded in Catchment 72014. To ease the visual identification of trends, least squares fits are superimposed to each plot. Overall, the regression slopes for

*Q*appear to be steeper than those for

*V*. Indeed,

*p*-values of the Mann–Kendall test indicate that only the trends of the annual and summer peak flow series in Catchment 70005 are statistically significant at a 5% significance level (Table 2).

Catchment . | Season . | Variable . | τ
. | p-value
. |
---|---|---|---|---|

70005 | Winter | Q | 0.230 | 0.067 |

V | 0.097 | 0.446 | ||

Summer | Q | 0.300 | 0.016 | |

V | 0.222 | 0.077 | ||

Annual | Q | 0.276 | 0.027 | |

V | 0.101 | 0.427 | ||

72014 | Winter | Q | 0.190 | 0.131 |

V | 0.214 | 0.089 | ||

Summer | Q | 0.173 | 0.168 | |

V | 0.133 | 0.292 | ||

Annual | Q | 0.157 | 0.212 | |

V | 0.099 | 0.436 |

Catchment . | Season . | Variable . | τ
. | p-value
. |
---|---|---|---|---|

70005 | Winter | Q | 0.230 | 0.067 |

V | 0.097 | 0.446 | ||

Summer | Q | 0.300 | 0.016 | |

V | 0.222 | 0.077 | ||

Annual | Q | 0.276 | 0.027 | |

V | 0.101 | 0.427 | ||

72014 | Winter | Q | 0.190 | 0.131 |

V | 0.214 | 0.089 | ||

Summer | Q | 0.173 | 0.168 | |

V | 0.133 | 0.292 | ||

Annual | Q | 0.157 | 0.212 | |

V | 0.099 | 0.436 |

Values in bold indicate statistically significant trends.

### Analysis of bivariate flood trends via Kendall's τ

Bivariate trends are assessed through an exploratory analysis of the relationship between *Q* and *V* using both graphical tools and a non-parametric measure on changes in the dependence between the variables, as estimated by Kendall's *τ*. Finally, hydrograph shape and its connection with Kendall's *τ* is also analysed. Note that the whole methodology followed in this section is first presented and later results are displayed.

First, a scatter plot of the pairs, where is the rank of and is the rank of (with ) is drawn for winter, summer and annual maximum flow events for both catchments under study. These plots provide a visual assessment of the dependence between variables. Note that scatter plots are linked to the Kendall's *τ* value, as the latter is a rank-based measure: more scattered ranks lead to smaller Kendall's *τ* values (i.e., lower association), whereas less scattered ranks will entail the opposite.

Next, changes in hydrograph shapes over time are investigated. The standardised mean flood hydrographs of the *Q-V* pairs for a given time window (whole series, W1 or W2) are plotted together by locating their time of peak on the same vertical (Miller *et al.* 2014), with the aim of visually comparing their average shape. Note that the mean flood hydrographs are standardised by the largest observed peak value () in the sample, considering the entire data length. The corresponding point-wise confidence intervals () are also obtained, where and *s* are the mean and standard deviation, respectively, of the standardised flow values at each time step. This analysis is carried out considering winter, summer and annual flood hydrographs.

*Q*and

*V*) and Kendall's

*τ*, can be understood by studying the relation between the corresponding hydrograph shapes and the pairs. A didactic example extracted from the annual data in Catchment 72014 is shown in Figure 6, in which specific pairs with are selected from the scatter plot (Figure 6(a)). Overall, three distinct clusters of events are evident. The first set is composed of pairs positioned close to the main diagonal (e.g., and ). The second set consists of pairs located below the main diagonal (e.g., and ). Finally, pairs located above the main diagonal (e.g., and ) constitute the third set. Each set is associated with a particular standardised hydrograph shape, as can be seen in Figure 6(b), generated following the procedure explained in the previous paragraph. Shapes related to points in the second set are steeper than the pairs of the first set. The opposite applies to the points of the third set, where shapes are less flashy. Therefore, the closer the pairs are to the main diagonal (i.e., the larger Kendall's

*τ*is, so the larger is the correlation), the more balanced and similar are the shape of the events. On the other hand, points located far away from the diagonal, which would lead to smaller values of Kendall's

*τ*, are generally characterised by a larger variability in the hydrograph shapes.

*Q*and

*V*in the case study is investigated in Figure 7. The scatter plots show the interplay of the scaled ranks of the two univariate variables, and they are generated by first using the complete record, followed by a plot considering each of the two time windows W1 and W2 (Figure 7(a)). Overall, it was found that ranks related to W1 are more scattered than to W2 for Catchment 70005. The value of Kendall's

*τ*is also derived for each time window and drawn in a new plot to facilitate the visual identification of the possible trend (Figure 7(b)). 95% confidence intervals are also displayed in Figure 7(b) (Schneider

*et al.*2015). The Kendall's

*τ*value is smaller for W1 (indicating a more weak correlation) than for W2 (indicating a stronger correlation) in Catchment 70005, as expected from the results of the scatter plots. Therefore, an increase over time was observed in Kendall's

*τ*for Catchment 70005. The opposite was found for Catchment 72014.

*τ*analysis presented in the previous paragraph, showing that the difference between confidence interval boundaries decreases from W1 to W2 for Catchment 70005 (i.e., the difference between events decreases), while the opposite holds for Catchment 72014 (Figure 8). No noticeable differences were identified between mean hydrograph shapes, in neither the time windows nor catchments.

### Trend significance assessment: a permutation procedure

*Q*, the volume

*V*or the Kendall's

*τ*coefficient are statistically different in the two time windows W1 and W2. This is slightly different from testing whether time is related to the variable of interest, as in the Mann–Kendall test, and allows the same procedure to be used to assess both univariate and bivariate temporal flood trends. Also, permutation tests are non-parametric (as well as the Mann–Kendall test) so that no formal distributional assumption for the data is needed in order for the test to be valid, and, in some situations, they provide exact inference (see, e.g., Ernst (2004) and Good (2005) for an introduction on permutation methods). The testing procedure consists of the following steps: (i) choose a test statistic which gives a good representation of the scientific question at hand; (ii) compute the test statistic for the observed data ; (iii) permute without replacement the observed sample for times; (iv) for each permuted sample compute the test statistic , with ; (v) estimate the empirical distribution () of the test statistic using all the permuted samples and compute the (two-tailed)

*p*-value as: For

*p*-values greater than 0.05, the null hypothesis is accepted at a 95% confidence level.

In this study, permutation methods are used to test whether the location of the distribution of both *Q* and *V* in the two different time windows is different. The difference between the sample median ( and ) in the two time windows is chosen as the test statistic to represent the null hypothesis of equal locations. The observed difference between the sample medians is compared with the distribution of the difference in the medians for the permuted samples and the associated *p*-value is computed as in Equation (1). Similarly, to further support the visual evidence of the previous sub-section, the null hypothesis of constant association between *Q* and *V* is tested by using the difference between Kendall's *τ* in the two windows as a test statistic.

*p*-values from Equation (1) by applying the permutation procedure to the case study are shown in Table 3. Among the positive observed differences between and in W1 and W2 for both catchments considering winter, summer and annual maximum flow events, only the difference associated with summer maximum flow events of Catchment 70005 can be considered significant. Note that this trend was also found to be significant when the Mann-Kendall test was applied considering the whole data length (first sub-section of ‘Methodology and results’ section). For the Kendall's *τ,* neither the increase for Catchment 70005 nor the decrease for Catchment 72014 were found to be statistically significant for any season.

. | . | Differences in . | ||
---|---|---|---|---|

Catchment . | Season . | m
. _{Q} | m
. _{V} | τ
. |

70005 | Winter | 0.264 | 0.644 | 0.123 |

Summer | 0.002 | 0.317 | 0.961 | |

Annual | 0.150 | 0.362 | 0.132 | |

72014 | Winter | 0.213 | 0.262 | 0.393 |

Summer | 0.332 | 0.256 | 0.355 | |

Annual | 0.262 | 0.450 | 0.465 |

. | . | Differences in . | ||
---|---|---|---|---|

Catchment . | Season . | m
. _{Q} | m
. _{V} | τ
. |

70005 | Winter | 0.264 | 0.644 | 0.123 |

Summer | 0.002 | 0.317 | 0.961 | |

Annual | 0.150 | 0.362 | 0.132 | |

72014 | Winter | 0.213 | 0.262 | 0.393 |

Summer | 0.332 | 0.256 | 0.355 | |

Annual | 0.262 | 0.450 | 0.465 |

Value in bold indicates statistically significant trends.

### Analysis of bivariate flood trends by comparison of return period curves

In the bivariate (*Q-V*) space, an infinite set of events given by their *Q-V* pairs are located under the same return period curve, which can be estimated by copulas (Salvadori & De Michele 2004; Requena *et al.* 2013). In this regard, the final step of the assessment of the impact of urbanisation on flood properties entails the analysis of trends in the bivariate *Q*-*V* space by the comparison between the return period curves in windows W1 and W2, for a set of given return period values. As an initial step, the analysis involves the selection of the joint distribution of (*Q*,*V*) that best characterises the statistical behaviour of the variables, composed of the marginal distributions of *Q* and *V* and a copula.

Although the information gained by studying the seasonal series can be useful for identifying potential changes in the flood seasonality, the analysis of the bivariate return period curves will focus on the annual series, since these are the ones that would be used to estimate design floods. The methods could also be applied to summer or winter maximum flow series.

### Selection of a joint distribution: Margins and copula

*Q*and

*V*, , can be expressed as: where

*q*and

*v*are given values of the variables

*Q*and

*V*, and are the marginal cumulative distributions functions of

*Q*and

*V*, respectively; and

*C*is the copula function, i.e., a joint cumulative distribution function with uniform margins. Thus, Equation (2) can be expressed as , where and . The selection of both the marginal distributions that best represent the individual variables, and the copula function that best characterises the dependence between

*Q*and

*V*is then needed to achieve a complete description of .

Several distributions used in hydrology, such as the GEV, generalised logistic (GLO), generalised normal (GNO), generalised Pareto (GPA) and PE3 were considered as potential marginal distributions of *Q* and *V*. The distribution that best characterises the observed data was selected using the L-moment ratio diagram, in which the relations between the theoretical coefficients of L-skewness () and L-kurtosis () for different three-parameter distributions are shown through curves (Hosking & Wallis 1997). The sample estimates of the coefficients of L-skewness () and L-kurtosis () are also plotted in the same diagram, choosing the distribution curve closest to the sample values. The choice of the marginal distribution for each variable was performed using the series covering the whole data length, and the selected distribution applied to both time periods. This larger sample ensures that estimates of and are less biased. Location , scale and shape parameters of the three-parameter marginal distributions ( and ) are estimated by the method of L-moments. Since estimates of the shape parameter have a high uncertainty when using a short record (Stedinger & Lu 1995), the estimate obtained using the complete data series was used in each of the two time windows to reduce the uncertainty originating from the third-order statistic estimates.

*Q*and

*V*, a representative set of potential copulas was tested: the Clayton, Frank and Gumbel copula belonging to the Archimedean family; the Galambos (and also the Gumbel) copula belonging to the extreme-value family and the Plackett copula representing other families. The goodness-of-fit test used for identifying possible copula candidates (see Genest

*et al.*2009) is based on the Cramér–von Mises statistic (): where is the empirical copula and is the estimated copula with a parameter (obtained by the inversion of Kendall's

*τ*method). The

*p*-value needed to formally test whether the copula is suitable is estimated by a validated bootstrap procedure (Genest & Rémillard 2008). This procedure (in a similar way to the aforementioned permutation procedure) derives an empirical distribution of the test statistic, , using 10,000 simulations. The copula is acceptable if the

*p*-value is greater than 0.05.

Additional information is needed to choose the best copula among the ones that have passed the goodness-of-fit test. By assessing the upper tail dependence, a further analysis to check the ability of each copula to characterise the extreme values of the studied variables is carried out (e.g., Poulin *et al.* 2007). For this purpose, the non-parametric upper tail dependence coefficient of the observed data, , is compared with the theoretical upper tail dependence coefficient, , of each copula (formulas in Frahm *et al.* (2005)). The copula is considered more appropriate the smaller is the difference between and values. Note that, however, the reliability of is limited, especially for a short data length. Moreover, the use of the Akaike information criterion (AIC) as model selection criterion (Akaike 1974) can be helpful for ranking the candidate copulas. Based on the latter, the best copula would be that with the smallest AIC value.

Although in the present case study the observed changes in time considering jointly *Q* and *V* (i.e., Kendall's *τ* trends in the sub-section ‘Trend significance assessment’) were not identified as statistically significant (with *p*-values equal to 0.132 and 0.465 for Catchment 70005 and 72014, respectively, see Table 3), they could still be hydrologically relevant. Therefore, the implications of such trends for the flood hydrograph shape should be analysed. For this analysis, the bivariate joint distribution of the observed data was estimated for both time windows.

*V*. The sample L-moment ratios of the

*V*series for each catchment were plotted on an L-moment ratio diagram (Figure 9). Figure 9 also shows the 95% confidence ellipsoids based on the bivariate distribution of L-skewness and L-kurtosis, which are the basis of the goodness-of-fit measure introduced by Kjeldsen & Prosdocimi (2015). The thicker lines indicate the distributions that can be accepted according to the goodness-of-fit measure; the selected distribution corresponds to the one for which the measure is minimised for each catchment. The GLO distribution was selected to represent

*V*for Catchment 70005, while the GEV distribution was chosen for Catchment 72014. The estimated parameters of the marginal distributions are shown in Table 4. As can be seen, the urban catchment has larger location parameters than the rural catchment, while for both catchments and variables the skewness is negative. The estimated fitted marginal flood frequency curves are displayed in Figure 10 along with the observed data series. In the case of Catchment 70005, the marginal curve corresponding to W1 intersects that corresponding to W2 for both

*Q*and

*V*. That is, if small univariate return periods (

*T*) are considered, larger values of

*Q*and

*V*are expected for W2; the contrary happens for larger

*T*values. Intersection between marginal curves is not observed in Catchment 72014.

. | . | . | . | Estimated margin parameters . | ||
---|---|---|---|---|---|---|

Catchment . | Variable . | Distribution . | Period . | . | . | . |

70005 | Q | GLO | 1976–1992 (W1) | 21.679 | 3.774 | −0.162 |

1993–2008 (W2) | 23.824 | 2.576 | ||||

1976–2008 (whole series) | 22.739 | 3.222 | ||||

V | GLO | 1976–1992 (W1) | 0.493 | 0.154 | −0.265 | |

1993–2008 (W2) | 0.497 | 0.141 | ||||

1976–2008 (whole series) | 0.496 | 0.146 | ||||

72014 | Q | GLO | 1976–1992 (W1) | 15.717 | 3.349 | −0.079 |

1993–2008 (W2) | 18.187 | 3.950 | ||||

1976–2008 (whole series) | 16.949 | 3.675 | ||||

V | GEV | 1976–1992 (W1) | 0.186 | 0.103 | −0.122 | |

1993–2008 (W2) | 0.214 | 0.114 | ||||

1976–2008 (whole series) | 0.201 | 0.107 |

. | . | . | . | Estimated margin parameters . | ||
---|---|---|---|---|---|---|

Catchment . | Variable . | Distribution . | Period . | . | . | . |

70005 | Q | GLO | 1976–1992 (W1) | 21.679 | 3.774 | −0.162 |

1993–2008 (W2) | 23.824 | 2.576 | ||||

1976–2008 (whole series) | 22.739 | 3.222 | ||||

V | GLO | 1976–1992 (W1) | 0.493 | 0.154 | −0.265 | |

1993–2008 (W2) | 0.497 | 0.141 | ||||

1976–2008 (whole series) | 0.496 | 0.146 | ||||

72014 | Q | GLO | 1976–1992 (W1) | 15.717 | 3.349 | −0.079 |

1993–2008 (W2) | 18.187 | 3.950 | ||||

1976–2008 (whole series) | 16.949 | 3.675 | ||||

V | GEV | 1976–1992 (W1) | 0.186 | 0.103 | −0.122 | |

1993–2008 (W2) | 0.214 | 0.114 | ||||

1976–2008 (whole series) | 0.201 | 0.107 |

Results for the copula selection criteria are shown in Table 5. Most of the considered copulas passed the goodness-of-fit (i.e., *p*-values greater than 0.05), with the exception of the Clayton copula that is rejected for several cases. Since all cases present upper tail dependence, i.e., is greater than zero, the previously accepted copulas are cut down according to the values. In this regard, the Gumbel and Galambos copulas show values close to . Consequently, both of them could be chosen, as the results of the AIC are also very similar. Owing to the results in Table 5 and its properties, the Gumbel copula was selected for the three time periods and both catchments. This copula was also selected as the best copula in characterising the dependence between peak flow and volume in other studies (e.g., Zhang & Singh 2006; Poulin *et al.* 2007; Karmakar & Simonovic 2009; Requena *et al.* 2013; Sraj *et al.* 2015).

Catchment . | Time period . | Copula . | . | S
. _{n} | p-value
. | . | . | AIC . |
---|---|---|---|---|---|---|---|---|

70005 | 1976–1992 (W1) | Clayton | 2.286 | 0.055 | 0.126 | 0.604 | 0 | −5.719 |

Frank | 6.377 | 0.043 | 0.454 | 0 | −8.870 | |||

Gumbel | 2.143 | 0.046 | 0.310 | 0.618 | −8.036 | |||

Galambos | 1.429 | 0.046 | 0.299 | 0.616 | −8.330 | |||

Plackett | 13.869 | 0.044 | 0.403 | 0 | −8.009 | |||

1993–2008 (W2) | Clayton | 7.231 | 0.049 | 0.066 | 0.803 | 0 | −3.634 | |

Frank | 16.636 | 0.035 | 0.439 | 0 | −25.318 | |||

Gumbel | 4.615 | 0.036 | 0.341 | 0.838 | −20.228 | |||

Galambos | 3.906 | 0.036 | 0.345 | 0.837 | −19.892 | |||

Plackett | 95.667 | 0.037 | 0.354 | 0 | −22.482 | |||

1976–2008 (whole series) | Clayton | 3.043 | 0.060 | 0.004 | 0.664 | 0 | −2.391 | |

Frank | 8.022 | 0.034 | 0.219 | 0 | −27.034 | |||

Gumbel | 2.522 | 0.031 | 0.244 | 0.684 | −26.040 | |||

Galambos | 1.810 | 0.031 | 0.241 | 0.682 | −26.101 | |||

Plackett | 21.622 | 0.034 | 0.196 | 0 | −25.201 | |||

72014 | 1976–1992 (W1) | Clayton | 3.371 | 0.039 | 0.544 | 0.686 | 0 | −16.647 |

Frank | 8.717 | 0.032 | 0.868 | 0 | −13.654 | |||

Gumbel | 2.685 | 0.031 | 0.899 | 0.706 | −12.614 | |||

Galambos | 1.974 | 0.031 | 0.903 | 0.704 | −12.794 | |||

Plackett | 25.516 | 0.032 | 0.889 | 0 | −12.932 | |||

1993–2008 (W2) | Clayton | 1.705 | 0.080 | 0.009 | 0.555 | 0 | −5.572 | |

Frank | 5.057 | 0.063 | 0.069 | 0 | −5.036 | |||

Gumbel | 1.853 | 0.056 | 0.116 | 0.546 | −5.055 | |||

Galambos | 1.136 | 0.056 | 0.123 | 0.543 | −4.796 | |||

Plackett | 9.102 | 0.062 | 0.081 | 0 | −5.400 | |||

1976–2008 (whole series) | Clayton | 2.062 | 0.060 | 0.010 | 0.575 | 0 | −20.427 | |

Frank | 5.876 | 0.041 | 0.113 | 0 | −17.136 | |||

Gumbel | 2.031 | 0.033 | 0.249 | 0.593 | −15.451 | |||

Galambos | 1.316 | 0.033 | 0.244 | 0.591 | −14.958 | |||

Plackett | 11.908 | 0.038 | 0.142 | 0 | −18.086 |

Catchment . | Time period . | Copula . | . | S
. _{n} | p-value
. | . | . | AIC . |
---|---|---|---|---|---|---|---|---|

70005 | 1976–1992 (W1) | Clayton | 2.286 | 0.055 | 0.126 | 0.604 | 0 | −5.719 |

Frank | 6.377 | 0.043 | 0.454 | 0 | −8.870 | |||

Gumbel | 2.143 | 0.046 | 0.310 | 0.618 | −8.036 | |||

Galambos | 1.429 | 0.046 | 0.299 | 0.616 | −8.330 | |||

Plackett | 13.869 | 0.044 | 0.403 | 0 | −8.009 | |||

1993–2008 (W2) | Clayton | 7.231 | 0.049 | 0.066 | 0.803 | 0 | −3.634 | |

Frank | 16.636 | 0.035 | 0.439 | 0 | −25.318 | |||

Gumbel | 4.615 | 0.036 | 0.341 | 0.838 | −20.228 | |||

Galambos | 3.906 | 0.036 | 0.345 | 0.837 | −19.892 | |||

Plackett | 95.667 | 0.037 | 0.354 | 0 | −22.482 | |||

1976–2008 (whole series) | Clayton | 3.043 | 0.060 | 0.004 | 0.664 | 0 | −2.391 | |

Frank | 8.022 | 0.034 | 0.219 | 0 | −27.034 | |||

Gumbel | 2.522 | 0.031 | 0.244 | 0.684 | −26.040 | |||

Galambos | 1.810 | 0.031 | 0.241 | 0.682 | −26.101 | |||

Plackett | 21.622 | 0.034 | 0.196 | 0 | −25.201 | |||

72014 | 1976–1992 (W1) | Clayton | 3.371 | 0.039 | 0.544 | 0.686 | 0 | −16.647 |

Frank | 8.717 | 0.032 | 0.868 | 0 | −13.654 | |||

Gumbel | 2.685 | 0.031 | 0.899 | 0.706 | −12.614 | |||

Galambos | 1.974 | 0.031 | 0.903 | 0.704 | −12.794 | |||

Plackett | 25.516 | 0.032 | 0.889 | 0 | −12.932 | |||

1993–2008 (W2) | Clayton | 1.705 | 0.080 | 0.009 | 0.555 | 0 | −5.572 | |

Frank | 5.057 | 0.063 | 0.069 | 0 | −5.036 | |||

Gumbel | 1.853 | 0.056 | 0.116 | 0.546 | −5.055 | |||

Galambos | 1.136 | 0.056 | 0.123 | 0.543 | −4.796 | |||

Plackett | 9.102 | 0.062 | 0.081 | 0 | −5.400 | |||

1976–2008 (whole series) | Clayton | 2.062 | 0.060 | 0.010 | 0.575 | 0 | −20.427 | |

Frank | 5.876 | 0.041 | 0.113 | 0 | −17.136 | |||

Gumbel | 2.031 | 0.033 | 0.249 | 0.593 | −15.451 | |||

Galambos | 1.316 | 0.033 | 0.244 | 0.591 | −14.958 | |||

Plackett | 11.908 | 0.038 | 0.142 | 0 | −18.086 |

Values in bold indicate copulas that pass the test.

### Comparison between bivariate joint return period curves

*q*or

*v*are exceeded by the random variable

*Q*‘or’

*V*, respectively (Equation (4)). where is the mean inter-arrival time between two successive events, with for annual maximum series. Finally, return period curves are obtained in original units by transforming the () pairs with the same into () pairs by Equation (5), using the previously selected marginal distributions: where is the inverse marginal distribution of

*Q*, . The same holds for . The bivariate return period curves associated with the separate time windows W1 and W2 and the whole data length are calculated.

*p*(i.e., points fulfilling , where

*p*is the simultaneous non-exceedance probability of the two variables) for both study catchments are shown in Figure 11. The plot also shows 100,000 () pairs randomly generated from the fitted copula (grey points), showing that smaller Kendall's

*τ*values lead to more scattered data. Results for W1, W2 and the entire record are plotted and compared for both catchments. In Catchment 70005, curves related to W2 (with larger Kendall's

*τ*) are located below those corresponding to W1 (with smaller Kendall's

*τ*), while converging in the extremes (i.e., the asymptotes). The opposite occurs for Catchment 72014. As expected, curves related to the whole data length are located between W1 and W2 curves.

*Q*than for

*V*, reflecting the findings for the univariate distributions (Figure 10). This means that flood events tend to have a lower peak value, while at the same time, the flood hydrographs tend to be less flashy. For instance, the vertex of the 100-year return period curve for the urban catchment undergoes a 9.4% decrease in peak and a 7.6% decrease in volume from W1 to W2. Also, in accordance with the results presented in Figure 10, such a trend differs for small return periods, as margins of different time windows cross at for

*Q*values, while at for

*V*values. However, the opposite is observed for Catchment 72014, as return period curves move upward to the right. Such a shift is larger for

*Q*, meaning that flood events would become larger, while flood hydrographs are steeper. This is also in accordance with the results obtained in the univariate case (Figure 10), where the increase of

*Q*is greater than that of

*V*. For instance, the vertex of the 100-year return period curve for the rural catchment undergoes a 20.3% increase in peak and a 14.1% increase in volume from W1 to W2. Note that for both catchments, a higher return period generally results in a larger shift of the curve. This could be caused by the increasing uncertainty related to increasing return periods, in particular with the small sample sizes available in this study.

In the case that an estimate of a design flood is required, the effect of changes in urbanisation extent on return period curves can be assessed using the results associated with W2. Alternatively, if stationarity is assumed, the most reliable estimate will be obtained by using the return period curves associated with the whole data length. Also, differences between return period curves in the two time windows (of Catchment 70005) give some insight into how the curves move in time because of the urban development. Consequently, this behaviour could be extrapolated in the future by using predictions of the urbanisation increase level expected in a given catchment.

## DISCUSSION

Significant univariate trends in the observed peak value *Q* for summer maximum flow events in the urban catchment (70005) have been identified by both the Mann–Kendall test and the permutation procedure. The Mann–Kendall test also found a significant trend in the observed annual maximum flow events in this catchment. No significant univariate trends were found, neither in the observed volume *V* nor in any of the variables extracted from the rural catchment (72014). The results of trends in the seasonal series help understanding from which types of events the trends in the annual series are most likely to be driven. The results suggest that for the urbanised catchment the potential changes in the annual series would mostly be driven by change for the summer events.

The visual analysis of bivariate (*Q*,*V*) series found an increase in Kendall's *τ* over time for Catchment 70005, leading to increasingly more regular hydrograph shapes; whereas the opposite was found for Catchment 72014, resulting in a larger variability of flood hydrograph shapes. The visual analysis of the hydrograph shape variability in time, using confidence intervals, confirmed this result.

The analysis of bivariate trends in the characteristics of flood events in the urban catchment found no statistically significant trends based on Kendall's *τ*. However, opposite results to those found in the rural catchment were obtained. The implications of such trends when considering bivariate return period curves were also opposite, suggesting that the trend in the urban catchment may be caused by changes in the flood generating processes that may not be statistically detected. Thus, it is likely that, if indeed a change is occurring in some of the flood characteristics, a larger sample size would be needed for standard statistical tests to detect it (see also Prosdocimi *et al.* (2014), for a discussion on sample size problems in the analysis of hydrological series). The selection and fit of both margins and copula would be more powerful and accurate if longer data series were available; and large uncertainties in the estimate of both margins and copula could have affected the results.

In addition, two larger flood events observed in the first time window (W1) in Catchment 70005 could lead to larger quantiles for high return periods than in the case of W2. Therefore, as W1 and W2 are short, both univariate and bivariate estimates for high return periods are highly dependent on the magnitude of the observed flood events and, consequently, on the magnitude of rainfall events that drove such flood events in each period. Finally, it is also interesting to highlight that the effect of the bivariate trend found (for high return periods) in the urban catchment is different from what would normally be expected in an urbanised catchment, as in general, urbanisation should lead to steeper hydrograph shapes. A possible explanation could also be the influence of the sewer system or local flood mitigation measures when high return period floods are considered, although looking into the causes is beyond the scope of this paper. Also, results may point towards a more complex interaction between urbanisation and flood characteristics than commonly assumed.

Note that in future studies, the results of the present preliminary analysis could be compared with those obtained by applying the proposed methodology to a much longer data length (when it is available), as well as to flood events identified by applying a detailed hydrograph separation method, for which an exhaustive analysis of rainfall and streamflow should be performed.

## CONCLUSIONS

A simple and general framework to investigate the effect of changes in a catchment land cover on the univariate and bivariate behaviour of some flood characteristics is introduced. The case study is composed of two nearby hydrologically and climatologically similar catchments in the northwest of England, where the most important difference is the increasing urbanisation extent in the urban catchment; hence any difference observed in time can mostly be attributed to urbanisation. In general, no statistical evidence of temporal change was identified in the univariate series, apart from an increasing trend in summer peak flows in the urban catchment. It should be mentioned that the permutation test used for trend significance assessment on the differences between the location of the distribution of a given variable might be applicable to other hydrological analyses.

The potential bivariate trend due to increasing urbanisation in the urban catchment was found to lead to smaller flood peaks and less flashy flood hydrographs. However, these results could be conditioned to the short available records and the use of larger data sets could be advisable for its confirmation. In addition, further research in the identification and modelling of the process control on storm runoff in urban catchments could help in understanding this finding.

The methodology presented in this work could be applied to any pair of catchments that can be considered hydrologically and climatologically similar except for one major characteristic, which has changed in one of the two catchments. Finally, the proposed methodology can help practitioners to describe trends in flood characteristics, in order to improve estimates of the design floods by a non-stationarity approach.

## ACKNOWLEDGEMENTS

This work has been supported by the COST Office grant ES0901 *European procedures for flood frequency estimation* (FloodFreq), via the Short Term Scientific Mission (STSM) program. The authors are also grateful for the financial contribution made by the *Carlos González Cruz* Foundation, the project ‘MODEX-Physically-based modelling of extreme hydrologic response under a probabilistic approach. Application to Dam Safety Analysis' (CGL2011-22868), funded by the Spanish Ministry of Science and Innovation (now the Ministry of Economy and Competitiveness) and the project CGL2014-52570 funded by the Spanish Ministry of Economy and Competitiveness. The authors thank James Miller for the urbanisation extent data shown in Figure 2.