Revisiting the Statistical Scaling of Annual Discharge Maxima at Daily Resolution with Respect to the Basin Size in the Light of Rainfall Climatology

Over the years, several studies have been carried out to investigate how the statistics of annual discharge maxima vary with the size of basins, with diverse findings regarding the observed type of scaling (i.e., simple scaling vs. multiscaling), especially in cases where the data originated from regions with significantly different hydroclimatic characteristics. In this context, an important question arises on how one can effectively conclude on an approximate type of statistical scaling of annual discharge maxima with respect to the basin size. The present study aims at addressing this question, using daily discharges from 805 catchments located in different parts of the United Kingdom, with at least 30 years of recordings. To do so, we isolate the effects of the catchment area and the local rainfall climatology, and examine how the statistics of the standardized discharge maxima vary with the basin scale. The obtained results show that: (a) the local rainfall climatology is an important contributor to the observed statistics of peak annual discharges, and (b) when the effects of the local rainfall climatology are properly isolated, the scaling of the standardized annual discharge maxima with the area of the catchment closely follows that commonly met in actual rainfields, deviating significantly from the simple scaling rule. The aforementioned findings explain to a large extent the diverse results obtained by previous studies in the absence of rainfall information, shedding light on the approximate type of scaling of annual discharge maxima with the basin size.


•
The statistical scaling of annual discharge maxima with the drainage area A of catchments is revisited in the light of rainfall climatology.

•
Neglecting the dependence of annual discharge maxima on rainfall characteristics may lead to diverse findings regarding the observed type of scaling.

•
The multiplicative fluctuations of discharge maxima relative to their corresponding means closely follow the scaling properties of actual rainfields. • For spatial scales below approximately 100 km 2 , where rainfall deviates from multiflactal scale invariance, a break of scaling is also observed for the multiplicative fluctuations.
In general terms, RFA includes the transfer of statistical information of hydrologic quantities between catchments in a homogeneous region (see e.g., the detailed review in Hosking and Wallis [13]), and exhibits a long history with starting point the index flood method initially developed by Dalrymple [1] in 1960, who made a first attempt to standardize the maximum annual discharges measured within a statistically homogeneous geographical region by the corresponding local means, and describe them using a single probability distribution model see e.g., [4,6,7,10,11,14,15]. Soon afterwards, USGS suggested a quantile regression technique where flood quantiles are regressed against catchment characteristics to develop regional flood prediction equations see [2,16], thus creating a shift towards the direction of modeling regional variability of river discharges based on rainfall climatology and catchment physical variables see e.g., [6,[17][18][19][20][21][22][23][24][25][26][27].
More than half a century later, river flow prediction in ungauged basins (PUB) still remains an open challenge in hydrology, as it becomes apparent also from the PUB initiative set by the International Association of Hydrological Sciences (IAHS) for the decade 2003-2012 (see https: //iahs.info/pub/index.php), with primary aim to reduce uncertainty in hydrological predictions.
In the context of PUB initiative, an attractive approach to model the statistical characteristics of maximum annual discharges at ungauged locations of a river network, is to use scale invariance arguments to link them to the statistical properties of maximum annual discharges measured at gauged locations within the same statistically homogeneous geographical region see e.g., [3,4,[8][9][10]14,17,[19][20][21][28][29][30][31][32]. More precisely, define Q (i) max (A i ) to be the annual maximum discharge at gauged location i of a river network that drains area A i , and Q ( j) max A j to be the annual maximum discharge at an ungauged location j upstream of i that drains area A j < A i . In the most general case of statistical scale invariance, referred to as stochastic self-similarity, multiscaling, or multifractality see e.g., [3,4,[7][8][9][10][11]14,17,19,20,28,30,32,33] one can obtain the statistical properties of Q ( j) max A j as a function of the statistical properties of Q (i) max (A i ) through (see e.g., the review in Veneziano and Langousis [32] and references therein): where G is a random function that depends on the ratio r = A i /A j and is stochastically independent from Q (i) max (A i ), and d = denotes equality in all finite dimensional distributions. Note that Equation (1) includes self-similarity (also referred to as simple scaling) as a sub-case when G is deterministic. It follows from Equation (1) that the moments E Q max (A) q of different orders q depend on the catchment size A in a log-linear way see e.g., [3,4,32]: where K(q) is a non-linear (linear) function of q in the multiscaling (simple-scaling) case. Note that Equation (2) is a necessary but not sufficient condition for stochastic self-similarity to hold, as it describes the marginal statistics of Q max as a function of scale A. More precisely, Equation (2) implies that: where md = denotes equality in the marginal distributions of Q max (i.e., obtained across independent basins of various areas A), whereas Equation (1) refers to all finite dimensional distributions of the maxima field (i.e., obtained across all sub-catchments of various areas A of a basin). Gupta and Waymire [3] (see also Lovejoy and Schertzer [34]) introduced the term "wide sense multiscaling" for the scaling of the marginals imposed by Equations (2) and (3), and "strict sense multiscaling" for the scaling imposed by Equation (1). In what follows, and since the present study deals with scaling of annual discharge maxima as a function of the drainage area A across geographically distinct basins, we skip the qualification term "wide sense" and refer to the scaling imposed solely by Equations (2) and (3) as multiscaling (simple scaling), in case the K(q) function in Equation (2) is a non-linear (linear) function of q.
Based on Equation (2), a convenient way to distinguish between multiscaling and simple scaling (or self-similarity) of annual discharge maxima, is by studying the form of the moment scaling function K(q), defined as: Where Q max (A) = Q max (A)/A . An alternative way that has been suggested in order to distinguish between these two different types of scaling, is to study the dependence of the coefficient of variation (CV) of Q max (A) as a function of the drainage area A see e.g., [4,6,7,10,11,14,17,19,33]. In the simple scaling case, CV is independent of A, whereas in the case of multiscaling, CV becomes a function of scale see e.g., [4,10,32].
The idea that the catchment area A can serve as the most critical factor in explaining flood variability has a long history, as briefly outlined below. The starting point was in the late 20th century, where new methods were developed to estimate design floods see e.g., [4,6,10,15,28,32,35,36]. For example, NERC [37] argued that E[Q] ∝ A θ with θ = 0.73 can explain approximately 69% of the sample variability of mean annual floods in Great Britain and Ireland. Cadavid [38] analyzed instantaneous streamflow records from different regions of the United States, pointing towards multiscaling behavior. Smith [14] studied flood records from the Appalachia, showing that Q max varies approximately log-linearly with A, with CV values displaying an ascending trend with A for basins smaller than a critical area A c 50 km 2 , and a descending trend for basins with larger areas.
In a subsequent study, Gupta et al. [4] studied the scaling properties of Q max using daily discharges from 270 drainage basins in central Appalachia, including the 104 catchments studied by Smith [14]. The study concluded that there is a break of scaling for catchment sizes larger than a critical area A c 50 km 2 , characterized by significant uncertainties due to the dispersion of the considered sample. Robinson and Sivapalan [7] used the CV concept outlined above, to analyze the same dataset as Gupta et al. [4]. They concluded that the CV values exhibited significant variability with a slight descending trend for catchments larger than a critical area A c 100 km 2 , and an ascending trend for catchments with A < A c .
A moment scaling analysis was also performed by Pandey [8] using streamflow records from 180 stations in Canada. The main conclusion of the study, obtained by studying how the moments E Q max (A) q of different orders q vary with the catchment area A, was that Q max displays simple-scaling behavior with A. Similar findings on the self-similarity of annual discharge maxima using moment scaling analysis have been reported by Bhatti [39], who studied annual discharge maxima from 2150 unregulated basins from the USGS database, Vogel and Sankarasubramanian [40] who used streamflow records from 1433 river basins across the United States, and Ishak et al. [11] who used streamflow data from 91 catchments from New South Wales (NSW) in Australia. Additionally, Dodov and Foufoula-Georgiou [10] provided further evidence on the self-similar behavior of peak annual discharges, by using 99 daily flow series from the US and showing that E Q max (A) q ∝ A θq exhibits two distinct simple-scaling regimes: θ ≈ 0.7 for A < 700 km 2 and θ ≈ 0.3 for A > 700 km 2 .
It follows from the discussion above that while some type of scaling of Q max (A) with A exists, one cannot easily conclude on its precise type. This is primarily due the type of the conducted analysis, and the large dispersion of the analyzed samples. Regarding the first issue, all aforementioned analyses were performed by studying the log-linearity of the moments of Q max (A) against A. Since Q max (A) increases almost proportionally with the drainage area A, the strong log-linear dependence of E Q max (A) q on A makes statistical interpretation of any type of non-linear behavior almost impossible from data. Regarding the second issue, note that for rainfall-triggered flood events, river discharges increase almost proportionally with both the drainage area A and the rainfall depthĪ. Hence, unless all catchments fall within the same hydrologically homogeneous geographical region (i.e., in terms of rainfall accumulations; see e.g., Gupta et al. [31]), regressing Q max solely against A cannot resolve the variability induced by the different hydroclimatologies and, more importantly, it may produce biased results in favor of a simple scaling rule (see discussion below). This explains to an extent the significant dispersion of the samples analyzed by the aforementioned studies, as well as some of the findings in favor of a simple scaling rule.
More precisely, according to the main assumption of the index flood method see e.g., [1,7], the discharge maxima at different locations within the same statistically homogeneous geographical region exhibit the same probability distribution when standardized by their mean, or some other index discharge. For example, considering that the mean discharge E Q ( j) at some location j increases almost proportionally with the annual rainfall depth, thus, being indicative of the local climatology, the index flood method can be written in the form: where Q ( j) max denotes the maximum annual discharge at location j. Assuming proportionality of discharges with rainfall intensity and some power θ (0, 1) of the drainage area (i.e., generalized rational method; see [15,26,[41][42][43][44]), if the mean annual rainfall intensity is assumed constant over the region of interest, then the assumption E Q ( j) /E Q (i) = A j /A i θ holds in good approximation, and the index flood approach in Equation (5) reduces to a simple scaling rule of annual discharge maxima with the drainage area A: Evidently, in the most general case when locations i, j exhibit different hydroclimatic characteristics, as is the case when investigating data originating from different regions, Equation (5) may lead to more complex types of scaling than the simple scaling rule in Equation (6).
It follows from the discussion above that, from a theoretical point of view, any type of moment scaling analysis of Q max with A cannot be conclusive regarding the actual type of scaling of Q max . This is because when regressing the moments of Q max against A, one a priori assumes a constant mean rainfall intensity field over the region of interest (see discussion on the derivation of Equation (6)) and, consequently, a simple scaling rule for the annual discharge maxima.
In this context, an important question arises on how one can effectively conclude on an approximate type of scaling of annual discharge maxima. Here we aim at answering this question, by investigating the effect of rainfall climatology on the spatial scaling of maximum annual floods using daily discharge data from 805 stations located in different parts of the United Kingdom. To the best of our knowledge, no such investigation has been conducted so far. This is done by: (1) studying how the mean value of the standardized discharge maxima Q max = Q max /A depends on the catchment area A and the average precipitation in 30-year climatic periods, hereafter referred to as SAAR (Standard-period Average Annual Rainfall), and (2) studying how the distribution of the ratio Q max /E[Q] , also referred to as index flood ratio (see e.g., Hosking and Wallis [13]) or amplification factor (see e.g., Langousis and Veneziano [45]), scales with the area A, along the lines of the index flood concept in Equation (5).
Section 2 provides necessary information regarding the data used. In Section 3, we start by illustrating the strong linkage between the rainfall characteristics of catchments and the statistics of annual discharge maxima, and then focus on the application of two alternative standardization approaches to study the scaling properties of annual discharge maxima with the size of basins A, while discussing the origins of the observed differences: (a) standardization by the catchment area A, through the ratio Q max = Q max /A (i.e., used in most studies focusing on the scaling of peak discharges; see above), and (b) standardization by the corresponding discharge mean, through the amplification factor γ max = Q max /E[Q] (i.e., following the original version of the index flood method in Equation (5)). A comprehensive summary of the main findings of the study, as well as future research directions are presented in Section 4.

Data
The hydrologic information used in this study originates from NRFA (National River Flow Archive; see https://nrfa.ceh.ac.uk/), an extensive open access database of flow information and catchment descriptors in the United Kingdom, hosted by the Centre for Ecology and Hydrology (CEH), holding over 60,000 station-years of data. Information used in this paper includes: (a) daily discharges from 805 catchments with at least 30 years of recordings, located in different parts of the United Kingdom, (b) catchment size information, and (c) average annual rainfall in mm for two standard 30-year periods; i.e., 1941-1970 and 1961-1990 referred to as SAAR1 and SAAR2, respectively. Catchment sizes vary from 0.9 to 9948 km 2 with a mean value of 359 km 2 . For the total of 805 stations, the average SAAR value (i.e., SAAR = (SAAR1 + SAAR2)/2), corresponding to the mean annual precipitation in mm for the 60-year period 1941-1990, varies from 558.5 to 2910.5 mm with a mean value of 1059.5 mm. Figure 1 shows the spatial distribution of the average SAAR for the 805 considered catchments across the United Kingdom. SAAR values tend to decrease when moving from the West to the East coast, with lowest values located in the southeastern region of the United Kingdom. The observed rainfall gradient is directly linked to Gulf Stream and the local topography, as moist air masses from the sea move upslope and cool causing precipitation to form see e.g., [44,46].

Data
The hydrologic information used in this study originates from NRFA (National River Flow Archive; see https://nrfa.ceh.ac.uk/), an extensive open access database of flow information and catchment descriptors in the United Kingdom, hosted by the Centre for Ecology and Hydrology (CEH), holding over 60,000 station-years of data. Information used in this paper includes: (a) daily discharges from 805 catchments with at least 30 years of recordings, located in different parts of the United Kingdom, (b) catchment size information, and (c) average annual rainfall in mm for two standard 30-year periods; i.e., 1941-1970 and 1961-1990 referred to as SAAR1 and SAAR2, respectively. Catchment sizes vary from 0.9 to 9948 km 2 with a mean value of 359 km 2 . For the total of 805 stations, the average SAAR value (i.e., SAAR = (SAAR1 + SAAR2)/2), corresponding to the mean annual precipitation in mm for the 60-year period 1941-1990, varies from 558.5 to 2910.5 mm with a mean value of 1059.5 mm. Figure 1 shows the spatial distribution of the average SAAR for the 805 considered catchments across the United Kingdom. SAAR values tend to decrease when moving from the West to the East coast, with lowest values located in the southeastern region of the United Kingdom. The observed rainfall gradient is directly linked to Gulf Stream and the local topography, as moist air masses from the sea move upslope and cool causing precipitation to form see e.g., [44,46].  by the catchment area A is important for the effects of local rainfall climatology to be revealed, as river discharges increase with both the basin size and accumulated rainfall.   The strong dependence of the statistics of ∕ on the local rainfall climatology explains to a large extent the diverse findings of previous studies see e.g., [8,10,11,28,39,40] for evidence on simple scaling, and [4,14,17] for evidence on multiscaling, when the effects of rainfall signature are not taken into account; see below.

Analysis and Results
A way to simultaneously isolate the effects of the area of the basin A, and the local rainfall characteristics when studying the scaling properties of annual discharge maxima, it is to follow the exact version of the index flood method in Equation (5)  , which removes the variability induced by the different sizes and rainfall characteristics of the considered catchments (see above), vary with the drainage area A of the basin. We do so by: (1) classifying the selected 805 catchments, according to their size, into 10 approximately equally sized groups (i.e., ~80 catchments per group; see Table 1), and (2) studying how the ensemble averages of the initial moments [( ′ ) ] and [( ) ] of different orders q inside each group vary with the average area Ā of the grouped catchments. Similar plots have been used in the literature to study the scaling properties of discharges as a function of the basin size see e.g., [4,[8][9][10][11]14,19,28,39,40]. The strong dependence of the statistics of Q max /A on the local rainfall climatology explains to a large extent the diverse findings of previous studies see e.g., [8,10,11,28,39,40] for evidence on simple scaling, and [4,14,17] for evidence on multiscaling, when the effects of rainfall signature are not taken into account; see below.
A way to simultaneously isolate the effects of the area of the basin A, and the local rainfall characteristics when studying the scaling properties of annual discharge maxima, it is to follow the exact version of the index flood method in Equation (5) (i.e., without additional assumptions on the functional dependence of E[Q] on A), standardize the annual discharge maxima Q ( j) max at each location j by the corresponding discharge mean E Q ( j) , and study the distribution of the amplification factor In what follows, we examine how the distributions of: (a) the standardized discharge maxima Q max = Q max /A (i.e., used in most studies focusing on the scaling of peak discharges; see Introduction), and (b) the amplification factor γ max = Q max /E[Q], which removes the variability induced by the different sizes and rainfall characteristics of the considered catchments (see above), vary with the drainage area A of the basin. We do so by: (1) classifying the selected 805 catchments, according to their size, into 10 approximately equally sized groups (i.e.,~80 catchments per group; see Table 1), and (2) studying how the ensemble averages of the initial moments E (Q max ) q and E (γ max ) q of different orders q inside each group vary with the average areaĀ of the grouped catchments. Similar plots have been used in the literature to study the scaling properties of discharges as a function of the basin size see e.g., [4,[8][9][10][11]14,19,28,39,40].

E[Q']
(a)   Figure 4 shows log-log plots of the empirical moments E (Q max ) q (points) of the standardized discharge maxima inside each group, as a function of the average areaĀ of the grouped catchments, for different moment orders q = 0.5, 1, 1.5, 2, 2.5, 3. The reason why we limit our analysis to moment orders q ≤ 3, is because for highly variable random fields (as is the case of discharge maxima), higher moment orders are underestimated with high probability see e.g., [32,[47][48][49][50][51][52][53]. Figure 4 shows log-log plots of the empirical moments [( ′ ) ] (points) of the standardized discharge maxima inside each group, as a function of the average area Ā of the grouped catchments, for different moment orders = 0.5, 1, 1.5, 2, 2.5, 3. The reason why we limit our analysis to moment orders  3, is because for highly variable random fields (as is the case of discharge maxima), higher moment orders are underestimated with high probability see e.g., [32,[47][48][49][50][51][52][53]. Similar to the findings of previous studies see e.g., [4,[8][9][10][11]14,19,28,40], one sees that for all moment orders considered, the initial moments of the standardized discharge maxima ′ vary log-linearly with the drainage area A (see least-squares (LS) fitted lines), with negative slopes corresponding to the empirical moment scaling function ( ) in Equation (4). Clearly, the shape of the ( ) function is non-linear (see circles in Figure 5), indicating significant deviations from the simple scaling rule. Since river discharges are the hydrological response of a catchment to the rainfall input, runoff maxima can be seen as the filtered output of a rainfall signal, and therefore any observed type of scaling is expected to be inherited (at least to some extent) from the statistical properties of the input see e.g., [32,[54][55][56]. Therefore, the observed deviations from simple scaling can be attributed to the multifractal structure of actual rainfields within finite but practically important ranges of scales: typically from below 1 hour to several days in time and from a few kilometers (say 9 to 10) to more than 100 km in space see e.g., [32,34,51,. For smaller spatial scales, a break of the log-linearity of the initial moments of areal rainfall with the scale of spatial averaging occurs see e.g., [32,50,61,62,73,85,87], with much smaller log-log slopes, indicating a power deficit of high frequency fluctuations. This break of scaling is not observable in Figure 4, as it is smeared out by the considerable variability of the rainfall climatologies of the analyzed basins; see discussion above.  Similar to the findings of previous studies see e.g., [4,[8][9][10][11]14,19,28,40], one sees that for all moment orders q considered, the initial moments of the standardized discharge maxima Q max vary log-linearly with the drainage area A (see least-squares (LS) fitted lines), with negative slopes corresponding to the empirical moment scaling function K(q) in Equation (4). Clearly, the shape of the K(q) function is non-linear (see circles in Figure 5), indicating significant deviations from the simple scaling rule. Since river discharges are the hydrological response of a catchment to the rainfall input, runoff maxima can be seen as the filtered output of a rainfall signal, and therefore any observed type of scaling is expected to be inherited (at least to some extent) from the statistical properties of the input see e.g., [32,[54][55][56]. Therefore, the observed deviations from simple scaling can be attributed to the multifractal structure of actual rainfields within finite but practically important ranges of scales: typically from below 1 hour to several days in time and from a few kilometers (say 9 to 10) to more than 100 km in space see e.g., [32,34,51,. For smaller spatial scales, a break of the log-linearity of the initial moments of areal rainfall with the scale of spatial averaging occurs see e.g., [32,50,61,62,73,85,87], with much smaller log-log slopes, indicating a power deficit of high frequency fluctuations. This break of scaling is not observable in Figure 4, as it is smeared out by the considerable variability of the rainfall climatologies of the analyzed basins; see discussion above. Along these lines, Figure 6 shows log-log plots of the empirical moments [( ) ] (points) of the amplification factors inside each group as a function of the average area Ā of the grouped catchments, for different moments = 0.5, 1, 1.5, 2, 2.5, 3. One sees that a break of the log-linearity occurs for spatial scales below approximately 100 km 2 , with much smaller slopes, which concurs with the observed break of scaling in spatial rainfall (see above). In addition, for spatial scales larger than 100 km 2 , the empirical moment scaling function ( ) in Equation (4) (see stars in Figure 5), estimated as the negative log-log slopes of the LS fitted lines in Figure 6, remains non-linear and close to that of the standardized discharge maxima ′ . Similar findings on the scaling properties of flood annual maxima recorded in Southern Italy have been obtained, also, by Iacobellis et al. [56]. The latter study reported a break of scaling around 100 km 2 and, also, showed that although significantly influenced by soil infiltration and abstraction Along these lines, Figure 6 shows log-log plots of the empirical moments E (γ max ) q (points) of the amplification factors inside each group as a function of the average areaĀ of the grouped catchments, for different moments q = 0.5, 1, 1.5, 2, 2.5, 3. One sees that a break of the log-linearity occurs for spatial scales below approximately 100 km 2 , with much smaller slopes, which concurs with the observed break of scaling in spatial rainfall (see above). In addition, for spatial scales larger than 100 km 2 , the empirical moment scaling function K(q) in Equation (4) (see stars in Figure 5), estimated as the negative log-log slopes of the LS fitted lines in Figure 6, remains non-linear and close to that of the standardized discharge maxima Q max . Along these lines, Figure 6 shows log-log plots of the empirical moments [( ) ] (points) of the amplification factors inside each group as a function of the average area Ā of the grouped catchments, for different moments = 0.5, 1, 1.5, 2, 2.5, 3. One sees that a break of the log-linearity occurs for spatial scales below approximately 100 km 2 , with much smaller slopes, which concurs with the observed break of scaling in spatial rainfall (see above). In addition, for spatial scales larger than 100 km 2 , the empirical moment scaling function ( ) in Equation (4) (see stars in Figure 5), estimated as the negative log-log slopes of the LS fitted lines in Figure 6, remains non-linear and close to that of the standardized discharge maxima ′ . Similar findings on the scaling properties of flood annual maxima recorded in Southern Italy have been obtained, also, by Iacobellis et al. [56]. The latter study reported a break of scaling around 100 km 2 and, also, showed that although significantly influenced by soil infiltration and abstraction  Similar findings on the scaling properties of flood annual maxima recorded in Southern Italy have been obtained, also, by Iacobellis et al. [56]. The latter study reported a break of scaling around 100 km 2 and, also, showed that although significantly influenced by soil infiltration and abstraction mechanisms, the dependence of the CV of flood annual maxima on the area A of basins is theoretically linked to rainfall scaling. In any case, a definite argument on the origins of the observed scaling cannot be reached in the absence of concurrent rainfall and discharge data at high temporal resolution (i.e., at least hourly).
Based on the aforementioned findings, one may conclude that by neglecting the dependence of the statistics of annual discharge maxima on the rainfall climatology, the dispersion/variability of the analyzed samples increases, with important consequences on the diversity of the corresponding findings, and the accuracy of the assessments made (see Introduction).

Conclusions
The present work aimed at investigating the marginal statistics of annual discharge maxima after proper standardization, so that the effects of the catchment area and the local rainfall climatology are effectively isolated. The analysis was conducted using daily discharge data from 805 stations located in different parts of the United Kingdom, with at least 30 years of recordings, by: (a) standardizing the annual discharge maxima Q ( j) max at each location j by the corresponding discharge mean E Q ( j) , following the exact version of the index flood method in Equation (5) (without additional assumptions on the functional dependence of E[Q] on A; see discussion on the derivation of Equation (6)), and (b) studying the scaling of the initial moments of the amplification factor γ max = Q max /E[Q] with the area A of the basin. Note that γ max corresponds to the multiplicative fluctuations of the annual discharge maxima relative to the discharge mean E[Q], with the latter been strongly linked to both the sizes and rainfall characteristics of catchments.
The obtained results show that the regional rainfall climatology significantly affects the mean value of the standardized discharge maxima Q max /A (see Figures 2 and 3b, and discussion in Section 3) and, therefore, by neglecting the dependence of the statistics of annual discharge maxima on rainfall, the dispersion/variability of the analyzed samples increases, with important consequences on the diversity of the obtained findings, and the accuracy of the assessments made. This explains to a large extent the diverse results of previous studies regarding the observed type of scaling (i.e., simple scaling vs. multiscaling) of peak annual discharges with the catchment area (see Introduction).
When the effects of the local rainfall climatology are properly isolated, through proper standardization of annual discharge maxima by the corresponding discharge means, the resulting multiplicative fluctuations closely follow the scaling properties of actual rainfields, deviating significantly from the simple scaling rule. In addition, for spatial scales below 100 km 2 where significant deviations of rainfall from multiflactal scale invariance have been observed see e.g., [32,50,61,62,73,85,87], a break of scaling is also observed for the amplification factor γ max (see Figure 6). This break is not observable when studying the initial moments of Q max (or the ratio Q max /A ) with A (see Figure 4), as it is smeared out by the considerable variability of the rainfall characteristics of basins.
While the presented analysis and findings shed light on the important role of local rainfall climatology on the statistical scaling of annual discharge maxima, one should not neglect the important influence of catchment properties and runoff routing mechanisms, as well as the space-time rainfall distribution embodied in different precipitation patterns see e.g., [88], especially at small temporal resolutions (below daily) where infiltration and abstraction effects become more influential see e.g., [54][55][56]. The aforementioned topics and, in particular, how the statistical scaling of annual discharge maxima with the size of basins depends on the characteristics of dominant rainfall generating mechanisms, and the scale of temporal averaging (i.e., resolution) of observed discharges, still constitute an open research challenge that remains to be addressed in the context of a large-scale study, when concurrent rainfall and discharge data at high temporal resolution (e.g., at least hourly) become openly available.