Interactive comment on “ Spatially dependent Intensity-Duration-Frequency curves to support the design of civil infrastructure systems

In general, the paper is well written. However, I have some concerns regarding the real contribution (novelty), connection with the literature and in particular with copula studies, as well as comparison with other models. Main comments: 1. Some important papers related to the topic are missing and more importantly the comparison with them not only in terms of results but also in terms of advantages and drawbacks (e.g. Bardossy and Pegram, 2009, Durocher et al. 2016 and Requena et al. 2018). 2. Regarding the issues motivating the study: the first one seems to be already fixed by Le et al. 2018b (as indicated on page 5), and the second issue is not clear (seems to be written as a statement not as an issue). 3. The topic can also be closely related to regional frequency analysis or estimation at ungauged basins. The authors did not make

1 Line 241: "An example of an asymptotically independent model is the inverted max-stable process (Wadsworth and Tawn, 2012). A general description of all continuous inverted max-stable processes that have standard exponential margins on a spatial domain is where are points of a unit Poisson process on (0, ∞) and the ( ) are independent replicas of a continuous, nonnegative stochastic process ( ) in the spatial domain , with { ( )} = 1 for all ∈ . It is convenient to work with a simple inverted max-stable process with unit Fréchet margins, because the marginal distribution can easily be transformed back to the GPD scale. To transform the process Ω( ) to unit Fréchet margins, the following transformation is used: Figure R1. The exceedance relationship of a conditional distribution across multiple durations based on the example in Figure 8 in the manuscript. The blue line is the relationship between 10-year unconditional return levels (at the location of the blue star in Figure 8) and durations, and the red line is the relationship between one in 10 chance conditional return levels (at the location of the blue star in Figure 8) and durations, given a 20-year event for 36 hr extremes happens at location of the red star (in Figure 8) for the centroid of the Kalang River catchment.

Major comment #2:
I'm puzzled about the GPD fits. If I understood correctly, GPD are fitted to 9 and 36 hr rainfall exceedances. If moving windows are considered, then there is a very strong auto-correlation for both the 9 and 36 hr rainfall values. Have you taken this into account in the fits? A declustering method should be applied. This may be the reason why the fits for 36 hr extremes are usually poorer than for 9 hr extremes (see Figs S5 and S6).
Response: Thank you. We did not consider moving windows; instead, we used restricted time periods for 36 hr rainfall (e.g. 01/01 00:00 to 02/01 12:00; 02/01 12:00 to 04/01 00:00; …). The use of a restricted estimates avoids the need for declustering to undo the effect of a moving window. We used a conversion factor of 1.13 to account for the difference between sliding (unrestricted) d hr rainfall maxima and restricted d hr maxima. This value is based on guidance from Australian Rainfall and Runoff (where Table  2.3.4. from Green et al. (2016) gives the 24-hr factor as 1.15 and the 48-hr factor as 1.11).

Response:
The green segments are to indicate the interpolation of the individual element failure probability to a system failure probability (discussion line 530). We have added this detail to the figure caption so the description is self contained. 9 Minor comment #31: l. 529: 1% annual exceedance prob → 1% AEP Response: Thank you. We have fixed this.

Minor comment #32:
l. 573: 1.74 → I guess this number depends on the considered levels Response: Yes, this number depends on the pair of locations that we analyse the conditional probability as well as the considered levels, so we have added a clarification of the considered levels in the revised version of the manuscript. 10 Minor comment #33:

l. 611: inverted max-stable → inverted max-stable process
Response: Thank you, we will fix it when revising the manuscript.  Table 1, I would have expected to have points at A=91, 294, 341, 771, 1020, which is not the case. Fig. S1 provides relationships between areal reduction factors (ARFs) and area (in km 2 ) for different return periods for the case study catchments. These relationships are calculated through the simulation of inverted Brown-Resnick process over equally sized grid points. To get the ARFs for each of subcatchments in the case study (corresponding to area A=91, 294, 341, 771, 1020), we need to interpolate these relationships. We will improve the explanation in the revised version of the manuscript. 9 Line 527: "The green lines help to interpolate the individual element failure probability from a given system failure probability of 1%. Both horizontal axis and vertical axis are constructed at a double log scale for viewing purposes." 10 Line 567: "for the two catchments having the strongest dependence (Fig. 7). The corresponding conditional flows were then estimated using a hydrological model WBNM and shown to be strongly related to the ratio of conditional and unconditional rainfall extremes (Fig. 9)."

Response to the Reviewer #3
In general, the paper is well written. However, I have some concerns regarding the real contribution (novelty), connection with the literature and in particular with copula studies, as well as comparison with other models. Main comments: Response: Thank you for your comments. We respond in detail below (your comments in italic font and our responses in normal font).

Major comment #1:
1. Some important papers related to the topic are missing and more importantly the comparison with them not only in terms of results but also in terms of advantages and drawbacks (e.g. Bardossy andPegram, 2009, Durocher et al. 2016 and Requena et al. 2018).
Response: Thank you for the suggestion. We have added discussion on these paper to the revised manuscript. 11

Major comment #2:
2. Regarding the issues motivating the study: the first one seems to be already fixed by  (as indicated on page 5), and the second issue is not clear (seems to be written as a statement not as an issue).
Response: Thank you for pointing this out. The second issue relates to the spatial properties of asymptotic dependence (explored in . While these two issues have been separately addressed in previous papers, the contribution is to show how to combine the methods to solve a realistic design problem. 3. The topic can also be closely related to regional frequency analysis or estimation at ungauged basins. The authors did not make this connection or show the difference. In the first case (similarity or connection), a huge literature exists and should be considered.
Response: Thanks for your comment. We have discussed differences to regional frequency analysis and methods of estimation in the revised manuscript. 12 12 Line 72: "Regional frequency analysis is one type of method to estimate IDF curves, where the precision of at-site estimates is improved by pooling data from sites in the surrounding region (Hosking and Wallis, 1997). These methods can be combined with spatial interpolation methods to estimate parameters for any ungauged location of interest (Carreau et al., 2013). To determine an effective mean depth of rainfall over a catchment with the same exceedance probability as at a gauge location, the pointwise estimate of extreme rainfall is multiplied by an areal reduction factor (ARF) . However, such methods do not account for information on the spatial dependence of extreme rainfall-whether for single storm duration, or for the more complex case of different durations across a region (Bernard, 1932;Koutsoyiannis et al., 1998). The lack of dependence prevents these approaches from being applied to estimate conditional or joint flood risk at multiple points in a catchment or across several catchments, as would be required for a civil infrastructure system." 4. The paper focused on a case study (a given set of data). However, the effect of some factors on the performance of the model as not discussed and not studied: for instance, and not limited to, the dimensionality (number of sites) and the size of the subgroups.

Response:
Thanks for your comment. This is beyond the scope of the current study.

Major comment #5:
5. An important missing element from the paper is the notion of copulas which is the most important when dealing with dependence. There is a huge literature in both hydrology and statistics (even in spatial dependence). I'm surprised to not see it in the paper.

Response:
We have added literature on copulas into the revised manuscript. 13

Major comment #6:
6. In section 4: why the GPD is used directly without model selection procedure? Why it is the same for all sites? The GPD is usually asymptotically justified which is not enough (and less justified in hydrology because of the sample sizes) and does not depend on the data at hand. It should be considered as a distribution among others (like GEV for block extremes).
Response: Thank you for this comment. We used the GPD because, in contrast to block maxima, it allows us to consider concurrent rainfall extremes and therefore enables the study of dependence. The intention in this paper is not to work through repetitive fitting of different distributions, but to demonstrate a plausible method based on joint rainfall extremes for the design of linear infrastructure. The same distribution is used at each site with variation at each site carried by the parameters. The marginal model adopted is not perfect, but it is plausible, and sufficient for the intent of showing the application of rainfall dependence to design.

Major comment #7:
7. Lines 245-248: please provide other alternative models and justify the choice of your model.
Response: Thank you. We have added justification of the choice of the Brown-Resnick model in the revised manuscript. For example,  show it has better performance than the extremal-t model. 14 Le, P. D., Davison, A. C., Engelke, S., Leonard, M., and Westra, S.: Dependence properties of spatial rainfall extremes and areal reduction factors, Journal of Hydrology, Submitted, 2018a.

Major comment #8:
8. The assumption, on page 11 line 215, is it reasonable? Is it verified in your case study?
Response: Thank you very much. The assumption of AEP neutrality in rainfall-runoff design is a standard assumption when using IDF curves. While the assumption is in widespread use, it is not without limitation as this issue was explored in to the following two papers. Response: The hydrological model (i.e. WBNM) is used to transform the conditional rainfall to conditional flow. A label has been added in the revised version of the manuscript to show this (on the arrow between the see the squares for Section 4.5 and Section 4.6 in the top-right of Figure 4).

Minor comment #1:
1. Fig 4: Why in the independent model, no fitting is required? What it means?
Response: Thank you for pointing this out. The term "the independent model" here is not clear. We have changed it to "the case of independence" and have clarified that we mean the case where rainfall extremes occur independently in space. Response: Thank you. We have more explicitly indicated that the comment on fitting relates to Figure 8 ( Figure 6 in the updated version). We have also emphasized that the main feature of the model shown in these figures is the relationship at h=0, for the case of dependence between two different durations at the same location. 17

Minor comment #8:
8. Line 538 : I'm not sure about this statement. It is not true in many situations.
Response: Thank you for your comment. We have restricted our commentary to conventional hydrological design that is based on IDF curves, which is more defensible than the original comment which was too general. By construction IDF curves are focused are point-wise estimators of extremes, thus a given design is focused on independent application of univariate statistics. Keywords: areal reduction factor, asymptotic independence, conditional probability, duration 10 dependence, extreme rainfall, flood probability, inverted max-stable process, joint probability, 11 spatially dependent Intensity-Duration-Frequency, 12 Abstract 13 Conventional flood risk methods typically focus on estimation at a single location, which is 14 inadequate for civil infrastructure systems such as road or railway infrastructure. This is because 15 rainfall extremes are spatially dependent, so that to understand overall system risk it is necessary to 16 assess the interconnected elements of the system jointly. For example, when designing evacuation 17 routes it is necessary to understand the risk of one part of the system failing given that another region 18 is flooded or exceeds the level at which evacuation becomes necessary. Similarly, failure of any single 19 part of a road section (e.g., a flooded river crossing) may lead to the wider system's failure (i.e. the 20 entire road becomes inoperable). This study demonstrates a spatially dependent Intensity-Duration-21 Frequency curve framework that can be used to estimate flood risk across multiple catchments, 22 accounting for dependence both in space and across different critical storm durations. The framework 23 is demonstrated via a case study of a highway upgrade, comprising five bridge crossings where the 24 upstream contributing catchments each have different times of concentration. The results show that 25 conditional and unconditional design flows can differ by a factor of two, highlighting the importance 26 of taking an integrated approach. There is also a reduction in the failure probability of the overall 27 system compared with the case of no spatial dependence between storms. The results demonstrate the 28 2 potential uses of spatially dependent Intensity-Duration-Frequency curves and suggest the need for 29 more conservative design estimates to take into account conditional risks. 30

Introduction 31
Methods for quantifying the flood risk of civil infrastructure systems such as road and rail networks 32 require considerably more information compared to traditional methods that focus on flood risk at a 33 point. For example, the design of evacuation routes requires the quantification of the risk that one part 34 of the system will fail at the same time that another region is flooded or exceeds the level at which conditional probabilities (e.g. the probability that one river is flooded given that the other river level 48 exceeds a specified threshold) and joint probabilities (e.g. the probability that one or both rivers are 49 flooded). However, continuous streamflow data are often not available at the locations most relevant 50 to the civil infrastructure system in question, or the catchment conditions have changed to a degree 51 that reflects historical streamflow records as unrepresentative of likely future risk. Thus, direct 52 application of streamflow data for flood risk quantification in civil infrastructure systems does not 53 represent a viable approach for the majority of situations., which can then be used to estimate both 54 conditional probabilities (e.g. the probability that one river is flooded given that the other river level 55 exceeds a specified threshold) and joint probabilities (e.g. the probability that one or both rivers are  To deal with these difficulties, two alternativeovercome common limitations of streamflow data, 63 rainfall-based approaches are commonly used. The firstOne method uses continuous rainfall data 64 (either historical or generated) to compute continuous streamflow data using a rainfall-runoff model 65 information on the spatial dependence of extreme rainfall-whether for single storm duration across a 92 region, or for the more complex case of different durations across a region (Bernard, 1932; 93 Koutsoyiannis et al., 1998). This prevents these approaches from being applied to estimate conditional 94 or joint flood risk at multiple points in a catchment or across several catchments as would be required 95 for a civil infrastructure system. 96 Although tailored multivariate approaches can be applied to estimate conditional and joint 97 probabilities of extreme rainfall for specific situations (e.g., Kao  given that it is not only necessary to preserve dependence of rainfall across space, but also to account 101 for dependence across storm burst durations, as different parts of the system may be vulnerable to 102 different critical duration storm events. To this end, arguably the most promising recent research 103 direction has been the application of max-stable process theory that is able to represent storm-level 104 dependence (de Haan, 1984;Schlather, 2002). This has been applied on a spatial domain by Padoan et 105 al. (2010), who calculated conditional probabilities for a spatial domain located in United States. 106 However, to ensure that this general approach can be applied for practical flood estimation problems, 107 two further problems need to be overcome: 108 1. The approach needs to not only account for spatial dependence for rainfall 'events' of a single 109 duration (e.g. the field of annual maximum daily rainfall data), but must also account for reduces with an increasing return period (Wadsworth and Tawn, 2012). This implies that 119 inverted max-stable models, which are asymptotically independent, are likely to be preferable 120 as an approach for representing spatially dependent IDF information. An added benefit of 121 correctly representing asymptotic dependence is that information on areal reduction factors 122 can be obtained directly from the model, rather than estimating ARF information 123 independently from the computation of the IDF curves. 124 This study addresses both these issues by demonstrating the application of the inverted max-stable 125 process to estimate joint and conditional probabilities of flood-producing rainfall in the form of 126 spatially dependent IDF curves. This approach adapts the methods developed by  to 127 inverted max-stable models, and then uses the derived spatially-dependent IDF curves combined with 128 the extracted information on AFRs as the basis for transforming the rainfall into flood flows.. 129 Regional frequency analysis is one type of method to estimate IDF curves, where the precision of at-130 site estimates is improved by pooling data from sites in the surrounding region (Hosking and Wallis, 131 1997). These methods can be combined with spatial interpolation methods to estimate parameters for 132 any ungauged location of interest (Carreau et al., 2013). To determine an effective mean depth of 133 rainfall over a catchment with the same exceedance probability as at a gauge location, the pointwise 134 estimate of extreme rainfall is multiplied by an areal reduction factor (ARF) . 135 However, such methods do not account for information on the spatial dependence of extreme 136 rainfall-whether for single storm duration, or for the more complex case of different durations across 137 a region (Bernard, 1932;Koutsoyiannis et al., 1998). The lack of dependence prevents these 7 approaches from being applied to estimate conditional or joint flood risk at multiple points in a 139 catchment or across several catchments, as would be required for a civil infrastructure system. 140 Although multivariate approaches can be tailored to estimate conditional and joint probabilities of 141 extreme rainfall for specific situations (e.g., Kao  and Singh (2007)), the development of a unified methodology that integrates with existing IDF-based 143 flood estimation approaches remains elusive. This is particularly challenging given that it is not only 144 necessary to preserve dependence of rainfall across space, but also to account for dependence across 145 storm burst durations, as different parts of the system may be vulnerable to different critical duration 146 storm events. To this end, max-stable process theory has been demonstrated to represent storm-level 147 dependence (de Haan, 1984;Schlather, 2002) and used to calculate conditional probabilities for a 148 spatial domain (Padoan et al., 2010). Copulas including the extremal-t copula (Demarta and McNeil, 149 2005), and the Husler-Reiss copula (Hüsler and Reiss, 1989) have also been used to model rainfall 150

dependence. 151
This study applies a max-stable approach with an emphasis on practical flood estimation problems: 152 1. The approach needs to account for, not only the spatial dependence of rainfall 'events' of a 153 single duration, but also the dependence across multiple durations. This was addressed by Le 154 et al. (2018b), who linked the max-stable model of Brown and Resnick (1977) with the 155 duration-dependent model of Koutsoyiannis et al. (1998), to create a model that could be used 156 to reflect dependencies between nearby catchments of different sizes. 157 2. Given that often the interest is in rare flood events, the model needs to capture appropriate 158 asymptotic properties of spatial dependence as the events become increasingly extreme. 159 Recent evidence is emerging that rainfall has an asymptotically independent characteristic (Le The case study is designed to address two related questions: (i) "What flood flow needs to be used to 169 design a bridge that will fail only once on average every times (e.g., = 10 for a 10-year event) 170 that a neighbouring catchment is flooded?"; and (ii) "What is the probability that the overall system 171 fails given that each bridge is designed to a specific exceedance probability event (e.g., the 1% annual 172 exceedance probability event)?" The method for resolving these questions represents a new paradigm 173 in which to estimate flood risk for engineering design, by focusing attention on the risk of the entire 174 system, rather than the risk of individual system elements in isolation. 175 In the remainder of the paper, Section 2 emphasises the need for spatially dependent IDF curves in 176 flood risk design, followed by Section 3 which outlines the case study and data used. Section 4 177 explains the methodology of the framework, including a method for analysing the spatial dependence 178 of extreme rainfall across different durations. It also includes an algorithm with which to use that 179 information in estimating the conditional and joint probabilities of floods. The results, and a 180 discussion on the behaviour of flood due to the spatial and duration dependence of rainfall extremes, 181 are provided in Section 5. Conclusions and recommendations follow in Section 6. 182

The need for spatially dependent IDF curves in flood risk estimation 183
The main limitation of conventional methods of flood risk estimation is that they isolate bursts of 184 rainfall and break the dependence structure of extreme rainfall. Figure 1 demonstrates a traditional  185 process of estimating at-site extreme rainfall for two locations (gauge 1, gauge 2) and three durations 186 (1, 3, and 5 hr) (Stedinger et al., 1993) (Stedinger et al., 1993). The process first involves extracting 187 the extreme burst of rainfall for each site, duration and year from the continuous rainfall data, and 188 then fitting a probability distribution (such as the Generalised Extreme Value (GEV) distribution) to 189 the extracted data. Figure 1 demonstrates that, through the process of converting the continuous 190 rainfall data to a series of discrete rainfall 'bursts', this process breaks both the dependence with 9 respect to duration and space. Firstly, the duration dependence is broken by extracting each duration 192 separately, whereas for the hypothetical storm in Fig. 1 it is clear that the annual maxima from some 193 of the extreme bursts come from the same storm. Secondly, the spatial dependence is broken because 194 each site is analysed independently. Again, for the hypothetical storm of Fig. 1 it can be seen that the 195 5 hr storm has occurred at the same time across the two catchments, and this information is lost in the 196 subsequent probability distribution curves. Lastly, there is cross-dependence in space and duration. 197

203
Having obtained the IDF curves for individual locations in Fig. 1, the next step is commonly to 204 convert this to spatial IDF maps by interpolating results between gauged locations. Figure 2 shows 205 hypothetical IDF curves from individual sites, with a separate spatial contour map usually provided 206 for each storm burst duration. In a conventional application the respective maps are used to estimate 207 the magnitude of extreme rainfall over catchments for a specified time of concentration. The IDF 208 curves are combined with an areal reduction factor (ARF) to determine the volume of rainfall over a 209 region (since rainfall is not simultaneously extreme at all locations over the region). However, 210 because the spatial dependence was broken in the analysis of IDF curves, the ARF come from a 211 separate analysis and are an attempt to correct for the broken spatial relationship within a catchment 212

222
The process in Fig. 1 breaks out the dependence of the observed rainfall, which makes the 223 conventional approach unable to analyse the dependence of flooding at two or more separate 224 locations. Instead, this paper advocates for spatially dependent IDF curves which are developed by 225 retaining the dependence of observed rainfall in the estimation of extremal rainfall. By applying 226 spatially dependent IDF curves to a rainfall-runoff model, the dependence of flooding between 227 separate locations can be achieved. 228

Case study and data 229
The region chosen for the case study is in the mid north coast region of New South Wales, Australia. The case study has five main catchments that are numbered in sequence in , it is nonetheless likely that some level of spatial dependence would exist and need to be 244 integrated into the risk calculations. This is particularly of relevance given extremal rainfall in this 245 region is strongly associated with 'east coast low' systems off the eastern coastline, whereby extreme 246 hourly rainfall bursts are often embedded in heavy multi-day rainfall events.   The black circles in Fig. 3 represent the sub-daily rain stations used for this study. There were 7 sub-

Methodology 261
This section provides the method used to estimate the conditional and joint probabilities of flood for 262 civil infrastructure systems based on rainfall extremes, which is explained according to the steps 263 shown in Fig. 4. First, the generalized Pareto distribution (GPD) is used as marginal distribution to fit 264 to observed rainfall for all duration at each locations (Section 4.1). After that, an inverted max-stable 265 process is introduced and then fitted to rainfall extremes of identical or different durations (Sections is followed by the simulation to calculate areal reduction factor (ARF) in Section 4.5. An event-based 268 rainfall-runoff model is employed in Section 4.6 to transform conditional rainfall to conditional flows. 269 With an assumption ofthat there is a one-to-one correspondence between rainfall intensity and flow 270 rate, the joint flood probability for the case study is equal to the joint probability of rainfall. An

Dependence model for spatial rainfall 294
Consider rainfall as a stationary stochastic process associated with a location in a region of 295 interest. Models (the notation for spatial extremes often use the conventionthe stochastic process is 296 simplified from ( ) to ). Without loss of transforming marginal values to generality it can be 297 assumed that the margins of have a unit Fréchet distribution. An important property of dependence 298 in the extremes is whether or not two variables are likely/unlikely to co-occur as the extremes become 299 rarer, as this can significantly influence the estimate of frequency for flood events of large magnitude. This is referred to as asymptotic dependence/independence, respectively. For the case of asymptotic 301 independence, the dependence structure becomes weaker as the extremal threshold increases, which is 302 formally defined as lim →∞ { 1 > | 2 > } = 0 for all 1 ≠ 2 . The spatial extent of a rainfall event 303 with asymptotically independent extremes will diminish as its rarity increases. 304 An example of an asymptotically independent model is the inverted max-stable process (Wadsworth 305 and Tawn

Fitting the dependence model 326
One simple way to calibrate dependence models is to fit them to data by matching a suitable statistic. 327 The dependence structure of the inverted max-stable process is represented by the pairwise residual 328 tail dependence coefficient (Ledford and Tawn, 1996) (Ledford and Tawn, 1996). 329 For a generic continuous process associated with a specific location the empirical pairwise 330 residual tail dependence coefficient for each pair of locations ( 1 , 2 ) is 331 The value of ∈ (0,1] indicates the level of extremal dependence between 1 and 2 (Coles et al., 333 1999) (Coles et al., 1999), with lower values indicating lower dependence. An example of how to 334 calculate the residual tail dependence coefficient is provided in Appendix A for a sample dataset. 335 To estimate the dependence structure of an inverted max-stable model, the theoretical residual tail 336 dependence coefficient function is usually fitted to its empirical counterpart. Here the residual tail 337 dependence coefficient function is assumed to only depend on the Euclidean distance between two 338 locations ℎ = ‖ 1 − 2 ‖. The theoretical residual tail dependence coefficient function for the inverted 339 Brown-Resnick model is given as: 340 where Φ is the standard normal cumulative distribution function, ℎ is the distance between two 342 locations, and (ℎ) belongs to the class of variograms (ℎ) = ‖ℎ‖ ⁄ for > 0 and ∈ (0,2). The 343 models are then fitted to the empirical residual tail dependence coefficients by modifying parameters 344 and until the sum of squared errors is minimized. 345 In the case that extreme rainfall at locations 1 and 2 are of identical duration (i.e. both 36 hr), then 346 the inverted max-stable process is fitted to the observations by minimizing the sum of the squared 347 errors of the residual tail dependence coefficients. This information can be directly applied to the case 348 where two catchments have a similar time of concentration owing to their similar shape and size.

Estimate of conditional and joint probabilities of rainfall extremes 376
The conditional probability { 2 > 2 | 1 > 1 } is obtained from the bivariate inverted max-stable 377 process cumulative distribution function (CDF) in unit Fréchet margins (Thibaud et al., 2013), which 378 is given as: 379 where 1 is the return period corresponding to the return level 1 . 390 This section introduces general concepts for evaluating a conditional probability and a joint 391 probability for a bivariate case. A detailed method is then presented for estimating the conditional 392 probability and the joint probability for the realistic case of rainfall extremes. 393

401
To explain how the joint and conditional probabilities are calculated, their definitions are provided in 402 Table 2 with reference to the regions of Fig. 5. Rather than consider the specific case of a theoretical 403 model of extremal rain (e.g. inverted max stable), Table 2 presents these concepts more simply using 404 only two variables and with generic probability estimates. Equations for both dependence and 405 independence are provided in Table 2. 406 Table 2. Definition of joint and conditional probabilities and how to calculate them for the case of bivariate independent and 407 dependent variables.

Case
Definition Calculation 1. Conditional prob. dependent The detailed formulae are of the same nature as those in Table 2, and are used in this study to estimate 413 conditional maps for return periods once the model has been fitted to all durations. 414 Case 2: Using the definition of { 2 > 2 | 1 > 1 } = { 1 > 1 , 2 > 2 } { 1 > 1 } ⁄ for the 415 independent case results in the exceedance probability for 2 , which is ( ) + ( ) (since intuitively 416 1 has no effect on exceedances of 2 ). 417 Case 3: For the case of dependent variables the joint exceedance is defined by ( ). For the case of 418 only two locations, the probability that there is at least one location that has an extreme event 419 exceeding a given threshold is calculated as { 1 > 1 or 2 > 2 } = { 1 > 1 } + { 2 > 2 } − 420 { 1 > 1 , 2 > 2 }. Here, { 1 > 1 , 2 > 2 } can be easily obtained from the bivariate CDF for 421 inverted max-stable process in Eq. (B.1). However, for the case of multiple locations (five different 422 locations for this paper), it is difficult to derive the formula for this probability because there are 423 dependences between extreme events at all locations. So this probability is empirically calculated 424 from a large number of simulations of the dependent model (see the description of the simulation 425 procedure for an inverted max-stable process in Section 4.5). It is also noted that the case study 426 contains five catchments, which have approximate times of concentration of either 36 hr or 9 hrs. 427 Case 4: Joint probability for The joint probability for independent variables is broken down as the 428 product of the marginals. The exceedance probability for 1 is ( ) + ( ) and the exceedance 429 probability for 2 is ( ) + ( ), and by definition their independent product will result in the joint 430 probability. In order to compare with a situation of no spatial dependence of rainfall extremes, theThe 431 probability that there is at least one location that has an extreme event exceeding a given threshold for 432 the case that all of events are independent can be calculated based on the addition rule for the union of 433 probabilities, as: 434

Evaluation of model for space-duration rainfall process 513
A GPD with an appropriate threshold was fitted to the observed rainfall data for 36 hr and 9 hr 514 durations, and the Brown-Resnick inverted max-stable process model was calibrated to determine the 515 spatial dependence. 516 Analysis of the rainfall records led to the selection of a threshold of 0.98 for all records as reasonable 517 across the spatial domain and the GPD was fitted to data above the selected threshold. Figure 75  518 shows QQ plots of the marginal estimates for a representative station for two durations 36 and 9 hr. 519 Overall the quality of fitted distributions is good and plots for all other stations can be found in the 520 supplementary material (Fig. S5S6 and S7).

524
The inverted max-stable process across different durations was calibrated to determine dependence 525 parameters. The theoretical pairwise residual tail dependence coefficient function between two 526 locations ( 1 and 2 ) was calculated based on Eq. (35) and Eq. (46), and the observed pairwise 527 residual tail dependence coefficient was calculated using Eq. (2). The model has a reasonable fit to 528 the observed data given the small number of dependence parameters.(4). Figure 86 shows the pairwise 529 residual tail dependence coefficients for the Brown-Resnick inverted max-stable process versus 530 distance. The black points are the observed pairwise residual tail dependence coefficients, while the 531 red lines are the fitted pairwise residual tail dependence coefficient functions. A coefficient equal to 1 532 indicates complete spatial dependence, and a value of 0.5 indicates complete spatial independence. 533 The top-left panel shows the dependence between 36 hr extremes across space, with the distance h = 0 534 corresponding to "complete dependence". It also shows the dependence decreasing with increasing 535 distance. Figure 6 indicates that the model has a reasonable fit to the observed data given the small 536 number of dependence parameters. Although the theoretical coefficient (red line) does not perfectly at 537 long distances, the main interest is in short distances, especially at ℎ = 0 for the case of dependence 538 between two different durations at the same location. 539

27
The remaining panels of Fig. 86 show the dependence of 36 vs. 9 hr extremes, 36 vs. 6 hr extremes, 540 and 36 vs. 3 hr extremes, with the latter two duration combinations not being used directly in the 541 study but nonetheless showing the model performance across several durations. As expected, the 542 dependence levels are weaker compared with 36 vs. 36 hr extremes at the same distance, especially at 543 thezero distance of 0. This is expected, as the dependence at the same site between annual 544 maximaexceedances at different durations will be lower than between annual maximaexceedances at 545 the same duration. This is because the annual maximaexceedances of different durations may arise 546 from different storm events (Zheng et al., 2015).  The recommended approach for estimating conditional rainfall extremes is demonstrated by 555 considering a hypothetical evacuation route across location 2 , given a flood occurs at location 1 , 556 evaluated using Eq. (B.39). This approach is applied to a case study of the Pacific Highway upgrade 557 project that contains five main river crossings (from Fig. 3). For evacuation purposes, we need to 558 know "what is the probability that a bridge fails only once on average every times (e.g., = 10 559 for a 10-yearan one in 10 chance conditional event) that its neighbouring bridge is flooded?" This 560 section provides the conditional estimates for two pairs of neighbouring bridges in the case study that 561 have the shortest Euclidean distances, i.e. pairs ( 1 , 2 ) and ( 2 , 3 ). The comparisons of 562 unconditional and conditional maps are given in Fig. 97 and Fig. 108, and the corresponding 563 unconditional and conditional flows are given in Fig. 11.9. In order to obtain the maps in Fig. 7 and 564 interest is the centroid of the Deep Creek catchment (the blue star in Fig. 108) and the conditional 579 point is the centroid of the Kalang River catchment (the red star in Fig. 108). The pointwise 10-year 580 unconditional and one in 10 chance conditional return levels at the location of the blue star are 134 581 mm and 194 mm, respectively. The relative difference between the conditional and unconditional 582 return levels is only 1.45 times, compared with 1.74 times for the case in Fig. 97. This is because the 583 pair of locations in Fig. 108 has a longer distance than those in Fig. 97, so that the dependence level is 584 weaker. Moreover, the location pair in Fig. 108 was analysed for different durations (between 36 and 585 9 hr extremes), which has weaker dependence than the case of the equivalent durations in Fig. 97

611
The left panel of Fig. 119 indicates that the peak conditional flow at the river crossing in the Bellinger 612 catchment is almost 2.0 times higher than that for unconditional flow. The time taken to reach to the 613 peaks is the same for both cases. This is because this river crossing is affected by a large region with a 614 long time of concentration (36 hr); the impact of rainfall losses on the hydrograph is insignificant. 615 This difference is a direct result of the conditional relationship being more stringent than the 616 unconditional relationship. Given that there is an existing extreme event nearby, it is more likely for 617 an extreme event to occur at another location of interest in the region. If a bridge design were to take 618 into account this extra criterion for the purposes of evacuation planning it would require the design to 619 be at a higher level. 620 Shown in the right panel in Fig. 119, the peak of the conditional flow at the river crossing in the Deep 621 Creek catchment occurred earlier, and is around 1.7 times higher than that for the unconditional flow. 622 This is due to the fact that the river crossing in Deep Creek covers a small region with a short time of 623 concentration (9 hr) and the impact of rainfall losses on the hydrograph is significant. 624 Although Fig. 119 shows a difference in terms of the time taken to reach the peak flows, the two 625 design hydrographs are separate and this is not a physical timing difference. The relevant feature of 626 the conditional design hydrograph is the peak, and timing information is not a part of the method. 627 The difference between the maximum discharge of conditional and unconditional flows at the river 628 crossing in the Bellinger catchment is shown in Fig. 12 for the case of a 20-year event occurring in the 629 Kalang River catchment nearby. The relationship with AEP indicates that the difference between the 630 maximum discharge of conditional and unconditional flows decreases when AEP increases, and that 631 the difference approaches zero when the AEP increases to above 50% (i.e. a 2-year return period). 632 633 Formatted: Font: 9 pt  The recommended approach for estimating the overall failure probability of a system is demonstrated 639 by considering a hypothetical traffic system with multiple river crossings at locations 1 , … , . If 640 there is a one-to-one correspondence between extreme rainfall intensity over a catchment and flood 641 magnitude, the overall failure probability will be approximately equal to the probability that there is at 642 least one river crossing whose contributing catchment has rainfall extremes exceeding the design 643 level, which can be estimated using a large number of simulations from the spatial rainfall model. 644 This approach is applied to the Pacific Highway upgrade project containing five river crossings. A set 645 of 10,000 year simulated rainfall (Section 4.5) is generated from the fitted model (Section 5.1) to 646 calculate the overall failure probability of the highway section. This process is repeated 100 times to 647 estimate the average failure probability, under the assumption that all river crossings are designed to 648 the same individual failure probability. 649 Figure 1310 is a plot of the overall failure probability of the highway andas a function of the failure 650 probability of each individual river crossing (black). Similar relationships for the cases of complete 651 dependence (blue) and complete independence (red) are also provided for comparison. For the case of 652 complete dependence, when the whole region is extreme at the same time, the overall failure 653 probability of the highway is equal to the individual river crossing failure probability and it represents 654 the best case. (the lowest overall failure probability). The worst case is complete independence where 655 extremes do not happen together unless by random chance; this means the failure probability of the 656 highway is much higher than that for individual river crossings. Taking into account the real 657 dependence, there are some extremes that align and it seems from the Fig. 1310 that this is a relatively 658 weak effect. As an example from Fig. 1310, to design the highway with a failure probability of 1% 659 To illustrate how Eq. (24) in the manuscript is calculated, consider a set of = 10 observed values at 739 the two locations: 1 = (5,9,1,2,10,3,8,6,4,7); 1 and 2 = (10,1,7,6,4,3,9,2,8,5). (see Table A1).  Table A1). 744 Table A1. Observed data 1 and 2 and corresponding empirical cumulative probabilities 1 and 2 .

Appendix B. Equations for bivariate conditional and joint probabilities for inverted max-stable 751
In the context of this study, the conditional probability { 2 > 2 | 1 > 1 } is obtained from the 752 bivariate inverted max-stable process cumulative distribution function (CDF) in unit Fréchet margins 753