Impacts of Spatial Heterogeneity and Temporal Non-Stationarity on Intensity-Duration-Frequency Estimates—A Case Study in a Mountainous California-Nevada Watershed

: Changes in extreme precipitation events may require revisions of civil engineering standards to prevent water infrastructures from performing below the designated guidelines. Climate change may invalidate the intensity-duration-frequency (IDF) computation that is based on the assumption of data stationarity. E ﬀ orts in evaluating non-stationarity in the annual maxima series are inadequate, mostly due to the lack of long data records and convenient methods for detecting trends in the higher moments. In this study, using downscaled high resolution climate simulations of the historical and future periods under di ﬀ erent carbon emission scenarios, we tested two solutions to obtain reliable IDFs under non-stationarity: (1) identify quasi-stationary time windows from the time series of interest to compute the IDF curves using data for the corresponding time windows; (2) introduce a parameter representing the trend in the means of the extreme value distributions. Focusing on a mountainous site, the Walker Watershed, the spatial heterogeneity and variability of IDFs or extremes are evaluated, particularly in terms of the terrain and elevation impacts. We compared observations-based IDFs that use the stationarity assumption with the two approaches that consider non-stationarity. The IDFs directly estimated based on the traditional stationarity assumption may underestimate the 100-year 24-h events by 10% to 60% towards the end of the century at most grids, resulting in signiﬁcant under-designing of the engineering infrastructure at the study site. Strong spatial heterogeneity and variability in the IDF estimates suggest a preference for using high resolution simulation data for the reliable estimation of exceedance probability over data from sparsely distributed weather stations. Discrepancies among the three IDFs analyses due to non-stationarity are comparable to the spatial variability of the IDFs, underscoring a need to use an ensemble of non-stationary approaches to achieve unbiased and comprehensive IDF estimates. IDF analyses were performed using the bias-corrected RESM outputs for each grid and compared with IDFs obtained using data from weather stations at six selected grid cells.


Introduction
Intensity-duration-frequency (IDF) curves have been widely used in water resources engineering to support the planning and design of hydrological infrastructure and facilities in flood risk and water management [1,2]. As rainfall-runoff modeling is typically used for flood estimation [3], IDF curves at the watershed scale are critical inputs to watershed models for designing storm water conveyance and flood control structures [4,5].

Method
This study used model outputs from the Regional Earth System Model (RESM) [25,26] developed at the Pacific Northwest National Laboratory (PNNL) based on the Weather Research and Forecasting (WRF) model [27] coupled with the Community Land Model (CLM) [28] through the flux coupler (CPL7) of the Community Earth System Model (CESM) [29] that facilitates an exchange of fluxes in a conservative manner. RESM was applied to a North American domain at 20 km grid resolution, with lateral boundary conditions and sea surface temperature and sea ice data provided by CESM. The CESM simulations are part of the Coupled Model Intercomparison Project Phase 5 (CMIP5) archive [30]. Driven by CESM boundary conditions, the RESM historical simulation covers the period of 1975-2004. Two future projections were simulated to cover the period of 2005-2100 following two representative concentration pathways (RCP), RCP4.5 and RCP8.5, corresponding to an emission mitigation scenario and a business-as-usual scenario, respectively. The bias correction spatial disaggregation (BCSD) method of Wood et al. (2004) was applied to the RESM historical and future simulations to provide bias-corrected time series of precipitation at a 0.125° latitude-longitude grid. Precipitation time series were then extracted from the bias-corrected RESM outputs to the 13 × 12 = 156 grid cells, with a latitude range of [37.9375, 39.3125] and a longitude range of [−119.8125, −118.3125], covering the Walker River Basin (see Figure 1). Among all the numerical grids, six were found to have adjacent weather stations, representing a variety of terrain and potentially different climate conditions. IDF analyses were performed using the bias-corrected RESM outputs for each grid and compared with IDFs obtained using data from weather stations at six selected grid cells.

Method
This study used model outputs from the Regional Earth System Model (RESM) [25,26] developed at the Pacific Northwest National Laboratory (PNNL) based on the Weather Research and Forecasting (WRF) model [27] coupled with the Community Land Model (CLM) [28] through the flux coupler (CPL7) of the Community Earth System Model (CESM) [29] that facilitates an exchange of fluxes in a conservative manner. RESM was applied to a North American domain at 20 km grid resolution, with lateral boundary conditions and sea surface temperature and sea ice data provided by CESM. The CESM simulations are part of the Coupled Model Intercomparison Project Phase 5 (CMIP5) archive [30]. Driven by CESM boundary conditions, the RESM historical simulation covers the period of 1975-2004. Two future projections were simulated to cover the period of 2005-2100 following two representative concentration pathways (RCP), RCP4.5 and RCP8.5, corresponding to an emission mitigation scenario and a business-as-usual scenario, respectively. The bias correction spatial disaggregation (BCSD) method of Wood et al. (2004) was applied to the RESM historical and future simulations to provide bias-corrected time series of precipitation at a 0.125 • latitude-longitude grid. Precipitation time series were then extracted from the bias-corrected RESM outputs to the 13 × 12 = 156 grid cells, with a latitude range of [37.9375, 39.3125] and a longitude range of [−119.8125, −118.3125], covering the Walker River Basin (see Figure 1). Among all the numerical grids, six were found to have adjacent weather stations, representing a variety of terrain and potentially different climate conditions. IDF analyses were performed using the bias-corrected RESM outputs for each grid and compared with IDFs obtained using data from weather stations at six selected grid cells.

Trend Analysis
Trend analysis can be conducted to evaluate the non-stationarity in the historical and future precipitation extremes. In this study, trend analyses including the Mann-Kendall (MK) test and Sen's slope analysis were conducted for the annual maximum series (AMS) of precipitation at each grid for different durations (6-and 24-h, which are the most commonly used event durations in extreme precipitation studies) and climate scenarios. The Mann-Kendall (MK) test [31] is a nonparametric statistical test for detecting trends in a time series. It can identify trends in the time series without assuming linearity [32]. The MK method has been widely-used for assessing the significance of trends in hydroclimatic time series data such as rainfall, temperature and stream flow [33][34][35][36][37][38]. Sen [39] developed a nonparametric procedure to estimate the slope of trends using linear model through pairs of sample points, which has become the classical method used to quantify trends in time series in hydrometeorology [40][41][42].

Quasi-Stationary and Non-Stationary IDFs
To obtain reliable IDFs, two solutions can be implemented: (1) identify quasi-stationary time windows from the time series of interest [43,44] and compute the IDF curves using data for the corresponding time windows; (2) introduce a parameter representing the trend in the means of the extreme value distributions [45][46][47][48]. Here, we adopt these approaches and derive both quasi-stationary and non-stationary IDFs for the study site. Extreme value theory has provided a framework to analyze the climate extremes and IDFs [49][50][51], which usually involves the fitting of generalized extreme value (GEV) distributions. The GEV distribution family consists of Gumbel, Fréchet and Weibull distributions. The standard cumulative distribution of GEV can be expressed as: The GEV distribution is flexible for modeling different behaviors of extremes with three distribution parameters: the location parameter µ, the scale parameter σ and shape parameter ξ. In order to represent a dynamic distribution, the location and scale parameters can be assumed to be linear functions of time to account for non-stationarity, with the shape parameter kept constant [52,53]. Then, µ and σ can be defined as where t is the time in years, x t denotes the intensity of the T-year return period event, K T is the frequency factor, µ a and σ a are the mean and the standard deviation of the set of AMS, respectively. For the simulated 6-and 24-h AMS under both climate scenarios at each modeling grid, GEV fitting was exanimated and the shape parameter ξ was generally found to be close to 0, which corresponds to the extreme value type I (EV1) distribution [54], which is also called the Gumbel distribution. Gumbel distribution is recommended by a number of studies as a good choice for developing IDFs [8,15,[54][55][56]. After carefully examining the quasi-stationary window, EV1/Gumbel is applied for developing IDFs for a 30-year time window. The corresponding lognormal regression equations are developed and used to estimate the events with 2-, 5-, 10-, 25-, 50-, and 100-year returning periods. The equation for fitting the Gumbel distribution to precipitation AMS for different return periods T is similar to Equation (2) [57][58][59], with the Chow's Lognormal frequency factor K T developed by [60]: In a stationary process, it is assumed that the probabilistic attributes of extreme events are independent of time. In a non-stationary process, the parameters of the underlying distribution function are time dependent and the properties of the distribution would vary with time. The framework we adopted to obtain non-stationary IDFs is called Non-stationary Extreme Value Analysis (NEVA) developed by [61]. The NEVA framework is based on the GEV and successfully used in nonstationary IDF analysis [47,48,62]. A Bayesian inversion framework was adopted in NEVA to infer the uncertainty of the GEV distribution parameters estimation under non-stationary conditions to yield more realistic estimations [63]. NEVA generates a large number of realizations from the parameter joint posterior distribution using the Differential Evolution Markov Chain (DE-MC) [64][65][66], which utilizes the generic algorithm DE for global optimization over the parameter space with the Markov-Chain Monte Carlo (MCMC) sampling approach.
Considering the non-stationarity in the data, two types of IDFs (i.e., quasi-stationary time windows based vs. trending parameter introduced non-stationary) were developed at all the numerical grids. The first IDF approach involves identification of quasi-stationary time windows [67]. The Mann-Kendall (MK) test failure rate was used to identify quasi-stationary time windows during the study period, and the 30-year-long time windows were determined to be the most stationary with lowest failure rate compared to other window sizes (e.g., 10-, 15-, 20-, 40-, and 50-year). After determining the quasi-stationary window, the IDFs were calculated for 30-year time periods including the historical period  and three future periods: 2011-2040, 2041-2070, and 2071-2100, respectively. The second non-stationary IDF approach is to apply the NEVA framework to the AMS for the entire time period (1975-2100), with 5000 DE-MC posterior samples in the Bayesian calibration of IDF parameters. Both the quasi-stationary and non-stationary IDFs were evaluated at each grid location, for durations of 6 and 24 h, for the 2-, 5-, 10-, 25-, 50-, and 100-year return periods, and for both RCP4.5 and RCP8.5 climate scenarios.

Geostatistics: (Variogram and Kriging)
Geostatistics has been used to evaluate and model spatial random functions, particularly when spatial heterogeneity and correlation patterns exist. It contains a suite of statistical methods including adaptations of classical regression techniques to describe the spatial continuity of natural phenomena. These methods aim to solve for characterization of spatial attributes by employing mainly random models in a similar fashion to the way time series analysis characterizes temporal data [68,69]. A package for geostatistical analysis named geoR [70] is used in our study to fit the variogram model for IDFs to reveal the spatial pattern of extreme events in the whole watershed.

Results
Bias-corrected RESM simulations of precipitation time series were used for IDF analyses. Trend analyses were done for the 126-year record (1975-2100) at each grid for the 6-and 24-h precipitation events under the RCP4.5 and RCP8.5 climate scenarios, as shown in Figure 2. The size of each dot represents the value of Sen's Slope at a confidence level set to be 0.95. The overall trending is all positive (i.e., increasing precipitation intensity) for different events under both climate scenarios. Generally stronger trending is observed for 6-h events than 24-h ones. The Sen's slope analysis shows the large spatial heterogeneity in the temporal trends, and that trending is extremely strong in mountain areas (e.g., northwest, southwest and middle of the region).
As a comparison, the IDFs with the traditional stationary assumption from historical observations were obtained at six adjacent weather stations from the Hydrometeorological Design Studies Center of NOAA's National Weather Service. Differences in the IDFs are shown at the selected grids highlighted in Figure 1. A comparison of the IDFs taken directly from the weather stations where the IDFs were obtained by the quasi-stationary approach for the same time period reflects biases in simulating precipitation by RESM, which has its own biases in addition to being influenced by biases in CESM1.
Although bias correction has been applied to the RESM outputs, it focuses on removing biases for monthly mean precipitation so biases may still exist for short-duration precipitation extremes. As a comparison, the IDFs with the traditional stationary assumption from historical observations were obtained at six adjacent weather stations from the Hydrometeorological Design Studies Center of NOAA's National Weather Service. Differences in the IDFs are shown at the selected grids highlighted in Figure 1. A comparison of the IDFs taken directly from the weather stations where the IDFs were obtained by the quasi-stationary approach for the same time period reflects biases in simulating precipitation by RESM, which has its own biases in addition to being influenced by biases in CESM1. Although bias correction has been applied to the RESM outputs, it focuses on removing biases for monthly mean precipitation so biases may still exist for short-duration precipitation extremes.
The IDFs estimated by NEVA provide uncertainty bounds from the ensemble estimates based on DE-MC sampling. The medians and uncertainty bounds of the ensemble IDF estimates are shown in the boxplots in Figures 3 and 4, for events with durations of 6 and 24 h, respectively. The IDFs estimated using the non-stationarity approach is potentially the closest to the truth and used as references; the terms of "underestimation" and "overestimation" are used to describe the relative deviations. It is observed that IDFs estimated using the quasi-stationary approach for the historical period  suggests that the simulated precipitation during this time period is generally comparable with the six adjacent station observations. However, clear differences between IDFs obtained by different approaches are also noticeable at specific locations (e.g., grid 2) where the simulated precipitation is much lower than observed, especially for the high return periods. Such The IDFs estimated by NEVA provide uncertainty bounds from the ensemble estimates based on DE-MC sampling. The medians and uncertainty bounds of the ensemble IDF estimates are shown in the boxplots in Figures 3 and 4, for events with durations of 6 and 24 h, respectively. The IDFs estimated using the non-stationarity approach is potentially the closest to the truth and used as references; the terms of "underestimation" and "overestimation" are used to describe the relative deviations. It is observed that IDFs estimated using the quasi-stationary approach for the historical period  suggests that the simulated precipitation during this time period is generally comparable with the six adjacent station observations. However, clear differences between IDFs obtained by different approaches are also noticeable at specific locations (e.g., grid 2) where the simulated precipitation is much lower than observed, especially for the high return periods. Such biases may be due to the lack of resolution in the model (20 km) as well as the gridded data (~12 km) that was used to bias correct the simulation in resolving the local conditions at Grid 2.
In general, the IDFs obtained from the weather stations deviate significantly from the ones estimated by the non-stationary approach and tend to underestimate the extremes especially for the 24-h duration events under both RCP4.5 and RCP8.5. For IDFs from the quasi-stationary approach, the estimates are higher for later time periods than earlier ones, although there are exceptions as they are confounded by spatial heterogeneity. For 6-h duration events (Figure 3), the IDFs based on the quasi-stationary approach are within the uncertainty range given by the non-stationary NEVA approach; however, the IDFs from the weather stations with traditional stationarity assumption overestimate by more than 100% for most return periods at Grid 2. They also tend to overestimate for the shorter return periods at grids 5 and 6, and underestimate for grids 1, 3 and 4, especially under RCP8.5. For example, under RCP8.5, the IDFs from the weather stations overestimate by about 20% on Grid 5, and underestimate 25% at Grid 3 for 2-and 5-year return periods. For the 24-h duration IDFs in Figure 4, the differences between the traditional stationary IDFs with the quasi-stationary IDFs for the same historical period are clearly smaller than the differences for the 6-h duration IDFs. This is consistent with the bias correction using monthly mean data, which should have larger impacts on removing biases for longer duration than shorter duration events. The traditional stationary IDFs from the weather stations tend to underestimate for longer return periods at all the selected grids except for grid 2. Under the RCP4.5 scenario at grids 1, 4 and 5, the estimates by both traditional and quasi-stationary IDFs underestimate by up to ~60% for the 100-year return period. The quasi-stationary IDFs for both time periods 2041-2070 and 2071-2100 are comparable to the non-stationary ones at grids 2, 3, and 6; on the other hand, the traditional stationary IDFs overestimate by ~30% (~8 mm/h) at grid 2 and underestimate by ~30% (~6 mm/h) at grid 3 compared to the non-stationary IDFs, which yield an estimate of ~6 mm/h and 9 mm/h, respectively.
Under the RCP8.5 scenario, the traditional stationary IDFs underestimate by ~40% at grids 1, 3, Spatial heterogeneity in the IDFs is evident from our results. When the IDFs obtained from quasi-and non-stationary approaches are considered, the extreme precipitation magnitudes at grids 1 and 3 are remarkably higher than at other grids; for example, under the RCP8.5 scenario, the precipitation estimated by non-stationary approach at grid 3 for the 100-year return period is about 17 mm/h, which is two or three times as high as it is at grids 2, 4, 5 and 6. The elevation of grid 3 is the highest (~3400 m) among the entire region and corresponds to the highest precipitation intensities. Grids 4, 5 and 6 are at relatively lower elevations (~1600 m) and have lower precipitation intensities, although their IDFs using quasi-stationary and non-stationary approaches do show different patterns (e.g., different slopes and uncertainty bounds). These differences might be attributed to the weather non-stationary IDFs are generally larger with higher emission inputs at all grids, and these bounds increase with the return period.
The above comparison shows a general agreement between the two non-stationary IDF approaches, although at a few locations (e.g., grids 1) larger deviations are observed. Such deviations might be statistically insignificant given the corresponding large uncertainty bounds in the estimated IDFs. In the following section, we fully evaluate the differences using the two non-stationary IDF analyses approaches, by looking at both deviations from the mean/median estimates and the uncertainty bounds of IDFs, for all time periods of study, at different spatial locations, for the 6-h and 24-h events of various return periods, under the two climate scenarios. The differences between non-stationary and quasi-stationary IDFs over the entire watershed are illustrated in Figures 5 and 6, where the boxplots show the distributions of relative difference percentage, computed as (non-stationary-quasi-stationary)/non-stationary*100%, at each numerical grid. Overall, in the majority of grids, the quasi-stationary framework produced IDFs that might underestimate the risk of extreme events for both short and long durations since the differences are positive for most boxplots. For the RCP4.5 scenario in Figure 5, the quasi-stationary window 2071-2100 has the smallest differences compared with the non-stationary IDF because the mean values of the difference approach zero for all return periods. However, under the RCP8.5 scenario, the mean values of differences are all positive, e.g., the relative differences in non-stationary IDFs for 24-h 100-year events are about 25% on average. The extra precipitation in a 24-h storm would cause a significant increase in flood peak. Another observation of the boxplots under both climate scenarios For the 24-h duration IDFs in Figure 4, the differences between the traditional stationary IDFs with the quasi-stationary IDFs for the same historical period are clearly smaller than the differences for the 6-h duration IDFs. This is consistent with the bias correction using monthly mean data, which should have larger impacts on removing biases for longer duration than shorter duration events. The traditional stationary IDFs from the weather stations tend to underestimate for longer return periods at all the selected grids except for grid 2. Under the RCP4.5 scenario at grids 1, 4 and 5, the estimates by both traditional and quasi-stationary IDFs underestimate by up to~60% for the 100-year return period. The quasi-stationary IDFs for both time periods 2041-2070 and 2071-2100 are comparable to the non-stationary ones at grids 2, 3, and 6; on the other hand, the traditional stationary IDFs overestimate by~30% (~8 mm/h) at grid 2 and underestimate by~30% (~6 mm/h) at grid 3 compared to the non-stationary IDFs, which yield an estimate of~6 mm/h and 9 mm/h, respectively.
Under the RCP8.5 scenario, the traditional stationary IDFs underestimate by~40% at grids 1, 3, and 4 for the 100-year return period, while the quasi-stationary IDF estimates for the period 2071-2100 match the non-stationary ones at all the selected grids. The largest deviation (~15% from the non-stationary estimation median) is observed at grid 1. The uncertainty bounds of the non-stationary IDFs are generally larger with higher emission inputs at all grids, and these bounds increase with the return period. The above comparison shows a general agreement between the two non-stationary IDF approaches, although at a few locations (e.g., grids 1) larger deviations are observed. Such deviations might be statistically insignificant given the corresponding large uncertainty bounds in the estimated IDFs. In the following section, we fully evaluate the differences using the two non-stationary IDF analyses approaches, by looking at both deviations from the mean/median estimates and the uncertainty bounds of IDFs, for all time periods of study, at different spatial locations, for the 6-h and 24-h events of various return periods, under the two climate scenarios.
The differences between non-stationary and quasi-stationary IDFs over the entire watershed are illustrated in Figures 5 and 6, where the boxplots show the distributions of relative difference percentage, computed as (non-stationary-quasi-stationary)/non-stationary*100%, at each numerical grid. Overall, in the majority of grids, the quasi-stationary framework produced IDFs that might underestimate the risk of extreme events for both short and long durations since the differences are positive for most boxplots. For the RCP4.5 scenario in Figure 5, the quasi-stationary window 2071-2100 has the smallest differences compared with the non-stationary IDF because the mean values of the difference approach zero for all return periods. However, under the RCP8.5 scenario, the mean values of differences are all positive, e.g., the relative differences in non-stationary IDFs for 24-h 100-year events are about 25% on average. The extra precipitation in a 24-h storm would cause a significant increase in flood peak. Another observation of the boxplots under both climate scenarios is that the uncertainty grows with the return period. Also, compared with the longer duration events, the uncertainty of shorter duration events is larger. is that the uncertainty grows with the return period. Also, compared with the longer duration events, the uncertainty of shorter duration events is larger.  Geostatistical analysis was performed for 100-year precipitation intensities with durations of 6 and 24 h, respectively. Spatial maps were generated to illustrate the spatial distributions of such events under different climate scenarios. Figure 7 illustrates the spatial intensity distribution of 6-h 100-year events using both stationary and non-stationary IDF estimations. The spatial patterns during different quasi-stationary time periods vary with the changing climate. In general, the intensities are stronger in mountain regions located near the mid-west of the watershed for both stationary and non-stationary IDF analyses. Under the RCP4.5 scenario, the peak intensity is 20 mm/h from non-stationary IDF, which is about 8 mm/h (40%) or 4 mm/h (20%) higher compared to stationary historical and future IDFs respectively. The quasi-stationary time period 2041-2070 corresponds to a peak intensity of 20 mm/h near the southwest border of the watershed, which is 6 mm/h (30%) higher than from the non-stationary IDF. With the increasing carbon emission, the intensities increase as well. In the non-stationary IDF maps with lower carbon emission RCP4.5, two peak intensity areas can be identified. Under RCP8.5, one area near the middle of the watershed, with intensity over 18 mm/h, has twice as high a precipitation intensity as that under RCP4.5. Also, the northwest region shows a 30% increase in intensity from 16 mm/h to 20 mm/h. For the quasi-stationary time period of 2071-2100, the area has the largest peak intensity corresponding to the highest carbon emission.
Lastly, we note that for the precipitation pattern in a small watershed, the decadal variability can be rather significant. For RCP4.5, the precipitation intensity is strongest during 2041-2070, but for RCP8.5, the precipitation intensity is strongest during 2071-2100, suggesting a simple linear scaling of precipitation intensity changes with the global mean warming is not likely to work well. On the other hand, despite the decadal fluctuations, the spatial pattern of precipitation intensity is Geostatistical analysis was performed for 100-year precipitation intensities with durations of 6 and 24 h, respectively. Spatial maps were generated to illustrate the spatial distributions of such events under different climate scenarios. Figure 7 illustrates the spatial intensity distribution of 6-h 100-year events using both stationary and non-stationary IDF estimations. The spatial patterns during different quasi-stationary time periods vary with the changing climate. In general, the intensities are stronger in mountain regions located near the mid-west of the watershed for both stationary and non-stationary IDF analyses. Under the RCP4.5 scenario, the peak intensity is 20 mm/h from non-stationary IDF, which is about 8 mm/h (40%) or 4 mm/h (20%) higher compared to stationary historical and future IDFs respectively. The quasi-stationary time period 2041-2070 corresponds to a peak intensity of 20 mm/h near the southwest border of the watershed, which is 6 mm/h (30%) higher than from the non-stationary IDF. With the increasing carbon emission, the intensities increase as well. In the non-stationary IDF maps with lower carbon emission RCP4.5, two peak intensity areas can be identified. Under RCP8.5, one area near the middle of the watershed, with intensity over 18 mm/h, has twice as high a precipitation intensity as that under RCP4.5. Also, the northwest region shows a 30% increase in intensity from 16 mm/h to 20 mm/h. For the quasi-stationary time period of 2071-2100, the area has the largest peak intensity corresponding to the highest carbon emission.
Lastly, we note that for the precipitation pattern in a small watershed, the decadal variability can be rather significant. For RCP4.5, the precipitation intensity is strongest during 2041-2070, but for RCP8.5, the precipitation intensity is strongest during 2071-2100, suggesting a simple linear scaling of precipitation intensity changes with the global mean warming is not likely to work well. On the other hand, despite the decadal fluctuations, the spatial pattern of precipitation intensity is dominated by the topographic variations and remains similar throughout the historical and future periods.
The spatial intensity distribution maps of 24-h 100-year events from both stationary and non-stationary IDF estimations are shown in Figure 8. Compared to 6-h 100-year events, the spatial patterns show some differences in both quasi-stationary and non-stationary IDF analyses. The west region of the watershed from the northern to the southern sides of the mountain forms a high intensity area that matches well with the terrain. The east side of the watershed with lower elevations has relatively weak intensities. The relationships between IDF curves and terrain/elevation were further evaluated, where the elevation for each numerical grid was extracted from the 3D Elevation Program (3DEP) provided by U.S. Geological Survey (USGS) [71]. The boxplot in Figure 9 displays the precipitation distribution for 6-and 24-h 100-year extreme events obtained by both quasi-stationary and non-stationary IDFs with elevations. The precipitation intensity increases with elevation in general.  The spatial intensity distribution maps of 24-h 100-year events from both stationary and non-stationary IDF estimations are shown in Figure 8. Compared to 6-h 100-year events, the spatial patterns show some differences in both quasi-stationary and non-stationary IDF analyses. The west region of the watershed from the northern to the southern sides of the mountain forms a high intensity area that matches well with the terrain. The east side of the watershed with lower elevations has relatively weak intensities. The relationships between IDF curves and terrain/elevation were further evaluated, where the elevation for each numerical grid was extracted from the 3D Elevation Program (3DEP) provided by U.S. Geological Survey (USGS) [71]. The boxplot in Figure 9 displays the precipitation distribution for 6-and 24-h 100-year extreme events obtained by both quasi-stationary and non-stationary IDFs with elevations. The precipitation intensity increases with elevation in general.    When elevation increases from 1200 m to 3300 m, precipitation intensity may increase from 5 mm/h to 10 mm/h for the 6-h events; and it may increase from 3 mm/h to 7 mm/h for the 24-h events. Precipitation intensity is subject to larger variability under the higher carbon emission scenario. The spatial heterogeneity due to the terrain effect is comparable to but generally larger than the variability due to the impact of a changing climate. For example, for the spatial locations with an elevation around 1800 m, the historical IDF yields a 5-6 mm/h intensity, while the future When elevation increases from 1200 m to 3300 m, precipitation intensity may increase from 5 mm/h to 10 mm/h for the 6-h events; and it may increase from 3 mm/h to 7 mm/h for the 24-h events. Precipitation intensity is subject to larger variability under the higher carbon emission scenario. The spatial heterogeneity due to the terrain effect is comparable to but generally larger than the variability due to the impact of a changing climate. For example, for the spatial locations with an elevation around 1800 m, the historical IDF yields a 5-6 mm/h intensity, while the future quasi-stationary or non-stationary IDFs yield a 5-10 mm/h intensity (RCP4.5) or 6-12 mm/h intensity (RCP8.5).

Discussion and Conclusions
In this study, both quasi-stationary and non-stationary IDF analyses over the entire Walker watershed were carried out using 126-year precipitation AMS time series under two RCP future emission scenarios, with high spatiotemporal resolution (hourly, 0.125 • ) climate model outputs. Traditional IDFs obtained from six adjacent weather stations were used for comparison, suggesting that the bias-corrected precipitation generally reproduces the observations well, except at the station on the southwest corner of the watershed where large differences were noted for 6-h duration events. With uncertainties in the physics parameterizations of climate models and insufficient spatial resolution to capture important driving factors for the local climate, the use of climate simulations to develop IDFs may introduce biases and uncertainties.
Despite some limitations of the climate model outputs, the long time series of precipitation provided at relatively high resolution by the models have been useful for exploring different approaches to account for the effects of climate non-stationarity on IDF estimates. We find that the quasi-stationary IDFs for the future time window of 2071-2100 match the non-stationary IDFs quite well, but the traditional stationary assumption-based estimates deviate significantly from the two, due to the strong non-stationarity in the data. The quasi-stationary IDFs for the earlier quasi-stationary windows 1975-2004 and 2011-2040, generally have lower estimates of extreme precipitations, as expected. The terrain effect results in strong spatial heterogeneity in the extreme precipitation distribution, and such an effect is at least comparable to the non-stationarity effect on the IDF estimate of extreme precipitation.
Overall, our results show that the discrepancies among the three IDFs analyses due to non-stationarity are comparable to the spatial variabilities of the IDFs. This underscores a need to use an ensemble of non-stationary approaches to achieve unbiased and comprehensive IDF estimates. But both quasi-and non-stationary IDFs are better solutions than the traditional stationary approach as they incorporate information of the impacts of climate change on extreme precipitation in flood estimation. For a comprehensive hydrological infrastructure engineering design, both approaches may be considered to handle non-stationary data and incorporate the associated estimation uncertainties in risk analysis, in addition to using outputs from higher spatial resolution climate models (e.g., 4 km) to deal with spatial heterogeneity and/or ensemble climate models to address model uncertainty.