Application of weather research and forecasting (WRF) model in evaluating depth–duration–frequency estimates under a future climate scenario

The impact of climate change on precipitation patterns is an area of interest for scientists and those designing and managing stormwater conveyance systems. The design of urban drainage systems is currently based on the statistical analysis of historical precipitation data (Chow 1953; Mailhot & Duchesne 2010; Forsee & Ahmad 2011) provided by the National Oceanic and Atmospheric Administration (NOAA). The NOAA National Weather Service (NWS) Hydrometeorological Design Studies Center used statistical analysis of historical data to produce point precipitation frequency estimates (PPFEs) available on the precipitation frequency data server (PFDS). The PPFEs published by the NOAA are viewed as a reliable resource for the design of stormwater management systems, including drainage networks and culverts (Forsee & Ahmad 2011; Merkel et al. 2014).

The applied methodology in generating the current PPFEs in NOAA Atlas 14 Volume 1 is based on stationary climate assumption (Bonnin et al. 2004). The regional observations were used to calculate the shape and scale parameters at gauge locations. The location parameter was based on at-station data. The annual maximum series (AMS) at each gauge station requires a minimum of 20 years period of record (POR) (Bonnin et al. 2011). The stationary climate assumption was validated in NOAA Atlas 14 through the evaluation of the linear trends in mean values and variance and shifts in mean values of 1-day annual maximum time series. This validation showed negligible effects of climate change on the AMS and the parameters that were used in the development of PPFE (Bonnin et al. 2011). The stationary method in NOAA Atlas 14 is based on the AMS extracted at the gauge locations. The AMS data used to develop the estimates for the State of Arizona extend up to December of 2000. However, it is reported in the literature that the intensity, frequency, and distribution of precipitation events are likely to be impacted as a result of heat-trapping greenhouse gas influx in the atmosphere (Guo 2006; Rosenberg et al. 2010; Forsee & Ahmad 2011; Moore et al. 2016).

A recent debate has focused on the circumstances under which an assumption of the non-stationarity should be considered in the statistical modeling of precipitation extremes (Milly et al. 2008; Montanari & Koutsoyiannis 2014; Serinaldi & Kilsby 2015). This consideration does not mean the stationarity which is also an attribute of a model is ‘dead’ (Bayazit 2015). The assumption of non-stationarity is more appropriate for the current study because it is targeting climate variability/change. However, the stationarity or non-stationarity assumption is decided based on how well they represent statistical estimations drawn from a time series (Koutsoyiannis & Montanari 2015; Dimitriadis & Koutsoyiannis 2018).

To examine the credibility of stationary climate assumption in designing water infrastructure under climate change, a variety of modeling and statistical methods have been applied in the literature. For example, Cheng & AghaKouchak (2014) showed that a stationary climate assumption would underestimate precipitation extremes by 60% in White Sands National Monument Station in New Mexico. In another study, Mailhot et al. (2007) performed statistical analysis for precipitation in southern Quebec using the Canadian Regional Climate Model (CRCM) precipitation outputs during 2041–2070. According to their results, the precipitation intensity would decrease in the future by 50% for the 2- and 6-h duration events and by 70% for 12- and 24-h duration events. Later, Mirhosseini et al. (2013) evaluated the impacts of climate change in the State of Alabama using precipitation outputs from six different climate models. According to their results, the rainfall intensity values were expected to decrease for short rainfall durations. On the other hand, for >4-h duration events, they performed probability-based analysis due to the large uncertainty associated with the projected rainfall. In a recent study by Shrestha et al. (2017), future (2046–2065) precipitation projections were generated using Long Ashton Research Station Weather Generator (a stochastic weather generator) and the rainfall disaggregation tool in Bangkok, Thailand. Their findings represented an increasing trend in future rainfall intensities for all the return periods.

In the state of Arizona as well, the spatial and temporal variability of the precipitation events has challenged the ability of scientists to forecast, model the state's atmosphere, and evaluate the impact of climate change on water infrastructures (Adams & Comrie 1997; Serrat-Capdevila et al. 2007; Neiman et al. 2013; Hughes et al. 2014; Luong et al. 2017; Mohebbi et al. 2018). The previous investigations have primarily utilized global climate models (GCMs) to produce high spatial and temporal resolution precipitation data. Generating precipitation data by the GCM is convenient, especially for data-scarce regions, such as the study area in this research, but it comes with a catch. GCMs have a questionable performance in modeling precipitation in contrast to modeling global temperature rise. They show false alarm simulations for near-zero precipitation, and they underestimate smaller temporal aggregates (Anagnostopoulos et al. 2010). Their outputs are typically in low resolution (hundreds of kilometers) and need to go under serious bias correction.

Some of the issues above are improved on a regional scale through a variety of downscaling methods. The most common downscaling methods include the stochastic weather generator, dynamic downscaling, and statistical downscaling (Prudhomme et al. 2002; Gutmann et al. 2012; Wang et al. 2013; Nourani et al. 2018; Vesely et al. 2019). According to the previous studies, the dynamic downscaling method produces more homogeneous spatial patterns than the statistical downscaling method and is, therefore, more appropriate for regional-scale limited-area models (Prudhomme et al. 2002). For example, Zobel et al. (2018) stated that dynamic downscaling has a better capability to simulate small-scale surface forcings such as topography and vegetation type, even with its inherent limitations documented in the literature. Other researchers such as Schmidli et al. (2007) reported that the statistical downscaling precipitation strongly underestimated the magnitude of interannual variations. Tolika et al. (2007) also confirmed that the use of statistical downscaling would underestimate the precipitation both in spring and winter. Moreover, Tang et al. (2016) reported complex and nonlinear relationships between precipitation and circulation as one of the major drawbacks of statistical downscaling in generating regional precipitation.

According to the reviewed literature mostly an underestimation but also an overestimation of the variability of precipitation is observed in some regions, while also the spatiotemporal precipitation trends may either increase or decrease (e.g., as observed through the discrepancies between the modeled depth–duration–frequency (DDF) curves and NOAA Atlas 14). While the choice of dynamic downscaling seems to be justified for a regional model, there are still certain concerns reported in the literature. The most crucial one is their ability in preserving the stochastic characteristics of fine-resolution precipitation (i.e., marginal and dependence structure). This means they may not correctly preserve the autocorrelation structure of rainfall (Lombardo et al. 2012) which is shown to exhibit fractal and long-range dependence behaviors at small (sub-daily) and large (over annual climatic) scales, respectively. This translates into impacting the output variability of the precipitation and its stochastic characteristics (Koutsoyiannis & Dimitriadis 2021). Precipitation's dependence (or simply autocorrelation) structure exhibits the so-called fractal behavior at small sub-daily scales and long-range dependence (LRD; or else known as Hurst–Kolmogorov behavior or scaling hypothesis or simply the Hurst phenomenon; Hurst 1951; Mandelbrot & Wallis 1968) at large climatic scales. Both behaviors (i.e., fractal and LRD), may highly increase the variability of precipitation (Koutsoyiannis & Dimitriadis 2021), and when combined to the heavy-tail of precipitation, then this variability further increases (Volpi et al. 2015), and on the hidden persistence in maximum rainfall records by Iliopoulou & Koutsoyiannis, 2019).

In this work, we applied a dynamic downscaling using the WRF model to produce the precipitation data for the medium climate change emission scenario (Nakicenovic et al. 2000) to extend our previous work (Mohebbi et al. 2021). The model developed by Mohebbi et al. (2021) was validated for a long-term historical data (1950–2017) minimizing the issues encountered in using GCM models for precipitation simulation. The goal of this work is to examine the climate variability impacts on future precipitation depth and frequency in the State of Arizona as a case study for an arid to semiarid climate. Assuming a non-stationary future precipitation model, we have generated PPFE and DDF curves for four different 6-year periods (2010–2015, 2030–2035, 2050–2055, and 2070–2075) at six locations in the study area to see how they deviate from the historical values published in NOAA Atlas 14.

This research was conducted based on the Intergovernmental Panel on Climate Change (IPCC) medium emission scenario, which is evaluated as the most-likely emission scenario if climate policies are not enforced (Wayne 2013). A high-resolution WRF model was developed for the state of Arizona. This model is a dynamically downscaled global climate projection system, and it is one of the most commonly used modeling systems in climate simulations (Harding & Snyder 2014; Mahoney et al. 2013). Using WRF-simulated precipitation data, the running total cumulative depth for each duration was calculated. Later, the maximum yearly precipitation depths of different durations were determined to generate the AMS. The Weibull distribution was then applied to calculate precipitation DDF. Similar to other distributions, the Weibull distribution is derived from the principle of maximum entropy but with two constraints (Mead & Papanicolaou 1984) and is related to a type III generalized extreme value (GEV) distribution. Weibull distribution has been recognized as one of the methods of precipitation frequency data estimation (Singh 1987; Heo et al. 2001; Bhunya et al. 2007). Another reason for using the Weibull distribution was its simplicity (only two constraints) compared to the other more sophisticated models available in the literature (Koutsoyiannis 2020). These two constraints, which are named shape and scale parameters, were first determined by fitting the AMS data to the Weibull probability density function (PDF) distribution (Clarke 2002; Serrano 2010; Peng & Yan 2014). Then, simulated PPFEs were generated by adopting Weibull's cumulative distribution function (CDF).

The University Corporation for Atmospheric Research (UCAR) developed the WRF model as a Unix-like open-source option for climate forecasting research (Michalakes et al. 2001). The dynamic numerical solver known as Advanced Research WRF (ARW) core was utilized in this study. The development of a WRF Preprocessing System (WPS) requires specifications on the model location, projection, spatial/temporal resolutions, and the time period. For this study, WPS models were set at a spatial resolution of 12 km×12 km and at a 1-min time resolution. The Mercator map projection was used for the models with a reference latitude and longitude of 34.5636 and −111.8543, respectively, which set the center of the model to the approximate central point in Camp Verde, Arizona. A resolution of 75 nodes in both vertical and horizontal directions was used to simulate an area of about 789,000 km² which is almost twice the size of the State of Arizona. The model was computed over a geographical domain larger than Arizona to ensure continuity at the edges of the state (Michalakes et al. 2001). The WRF model uses the physics suite option ‘CONUS’ as the modeling baseline condition that matches the continental United States. CONUS applies the same physics parameters as the National Weather Service (NWS) uses to create their United States national forecasts. The physics options and the main items of the CONUS parameterization are summarized in Table 2 (Wang et al. 2016). The numerical solver in WRF's ARW core is coded based on the third-order Runge Kutta method to resolve temporal processes (Wang et al. 2016).

The WRF model used in this study has undergone rigorous optimization and calibration in the team's previous work (Mohebbi et al. 2021). The PPFEs were determined based on 6 years of WRF-simulated data for every 20 years starting from 2030. The 6-year time duration was selected to showcase the climate change effects on precipitation events in 20-year increments as the purpose here is mainly to assess the variations of precipitation depths and frequency with respect to varying climatic parameters rather than offering a new set of PPFE values or DDF curves. Longer than 6-year period simulations would considerably increase the simulation runtime because of the limitations inherent in the WRF dynamic solver and the low priority of the simulation in the utilized high-performance computing system (6-year runtime was 3 weeks for each period). In addition, the storage required for 6 years of the high-resolution simulation was 4.7 terabytes, a little shy of the five terabytes available storage.

The WRF model input data for both past and future were adopted from the National Center for Atmospheric Research (NCAR 1994). The historical data were chosen from 2010 to 2015 to compare NOAA's reported PPFEs to the ones based on the simulated data for current climate conditions. The future prediction model was conducted for the medium emission scenario for 2030–2035, 2050–2055, and 2070–2075 based on the representative concentration pathways (RCPs) of 6.0, which corresponds to 850 parts per million of carbon dioxide equivalent (Griggs & Noguer 2002).

The forcing data used to simulate the RCP 6.0 model include the fifth phase of Coupled Model Intercomparison Project (CIMP5). However, using GCMs such as CMIP5 for regionally downscaled models can result in high bias. That is why regional studies typically have uncertainty analysis to account for the bias in their output. This is also one of the disadvantages of using a dynamical downscaling which was pointed out before. Having this in mind, the forcing data used in this study were corrected for bias by NCAR (Bruyère et al. 2014). The new data were published under the name of NCAR's Community Earth System Model (CESM) Global Bias-Corrected CMIP5 Output to support WRF research (Monaghan et al. 2014).

This study used a WRF model for the state of Arizona to produce a high-resolution gridded precipitation data set. The yearly maximum rainfall depth for different precipitation durations was extracted from the simulated data to calculate the AMS. The AMS data were then fit to a Weibull distribution to find the scale and shape parameters that best explain the probability distribution of extreme values. The AMS was generated for each time interval at the following durations: 1, 2, 3, 6, 12, and 24 h; 2 days, 3 days, 4 days, 7 days, 10 days, 20 days, 30 days, 45 days, and 60 days. These series were used to determine the Weibull parameters.

The PPFEs for each of the selected stations were also presented in the form of DDF curves. Each annual exceedance probability was plotted with precipitation depth on the -axis versus duration (⁠-axis) set as a log scale. The simulated DDF curves at each station for each time interval were plotted along with the NOAA DDF curves for comparison purposes.

Comparisons of observed data and simulations are also provided during the monsoon season (Figure 2) as the North America Monsoon (NAM) is most prevalent in the state of Arizona. The ability to project monthly precipitations during the monsoon season provides a better understanding of extreme precipitation events. The NAM is the result of organized convection that occurs at lower elevations of northwestern Mexico and the southwestern United States in the late summer months (Smith & Gall 1989; McCollum et al. 1995; Wall et al. 2012; Mazon et al. 2016). About 30–50% of the annual precipitation in the southwestern United States and an even more significant portion in northwestern Mexico occurs during the NAM season (Mazon et al. 2016). The NAM is responsible for weather hazards such as hail, lightning, dust storms, microbursts, and flash flooding (Mazon et al. 2016). The WRF-simulated precipitations for the monsoon season resulted in an R² of 0.74, NSE of 0.35, and RMSE of 37.85 mm. The discrepancies between the modeled and observed precipitation data during the monsoon season can be related to the underlying model physics in WRF, as well as high convention and unprecedented changes associated with monsoon season (Ragi et al. 2020). A couple of studies have attempted to address the issues and performance of dynamically downscaled models in forecasting the timing and intensity of the summer monsoon precipitations (Singh et al. 2015; Ragi et al. 2020). Since sub-hourly precipitation predictions associated with the monsoon season are subject to high uncertainty, the data comparison is not discussed for these events in this study.

The modeled PPFEs and the corresponding values published in NOAA Atlas 14 are demonstrated for each of the selected stations in Tables S1-S30 in the supplemental materials. The stations used in this study are representative of the climate variation experienced in the state of Arizona, with the driest and wettest regions represented by Yuma Valley and Payson, respectively. The scale parameter, a, was in the range of 1–155 mm, and the shape parameter, b, was between 1 and 5 by fitting the AMS data to Weibull distribution in different stations. As an example, the PPFEs for the Phoenix-Encanto are displayed in Table 4 (modeled and NOAA Atlas 14 values) from 2070 to 2075. Phenix-Encanto station is considered because of the proximity to the city of Phoenix, the most populous city in the state. The first column of Table 4 lists the durations of 1, 6, 24 h; 3 days, 10 days, and 60 days; the first row is an annual exceedance probability or a probability of specific precipitation recurrence interval.

To better understand the discrepancies between the modeled and published PPFEs, a comparison of the 6-h duration precipitation with annual exceedance probability is provided in Table 5. The 1/(100 years) annual exceedance probability, related to the 100-year recurrence interval, is commonly used in the design of storm drainage systems, open channels, and culverts (City of Flagstaff 2009; Pima County Regional Flood Control District 2014; Maricopa County Flood Control 2018). As shown in Table 5, for the most part, the depth estimates based on the WRF-simulated data were lower than those in NOAA Atlas 14 except for the 2070–2075 time period. As it was mentioned earlier, the future spatiotemporal trend of precipitation was not similar everywhere, and correspondingly the PPFE values at each interval were not following a specific pattern, either. However, on average, an overall decreasing trend can be suggested for depth values.

The Phoenix-Encanto station, which represents the most densely populated region, had similar results to the Payson station (Figure 6). Except for 2010–2015, where much lower depth values were determined based on the simulated data compared to the NOAA-published DDF curves. The average precipitation was predicted to decrease during the years 2010–2015 (Figure 4), which could be the reason for reducing depth in DDF curves for all the durations.

Regarding the DDF curves in the driest region of Arizona, represented by Station 2 in Yuma Valley (Figure 7), the NOAA DDF curves were mostly higher for all durations. In contrast to the wettest and most populated regions, where long-duration storms resulted in higher depth values, in Yuma Valley, the future long-duration events were either equal or less intense than the historical amounts.

Generally, the discrepancies between the modeled DDF curves and NOAA Atlas 14 were mostly observed for sub-daily precipitation events. Other studies have also found that climate change mainly influences shorter-duration rainfall events while longer-duration events experience little or no change in precipitation depth (Trenberth et al. 2003; Willems & Vrac 2011; Mirhosseini et al. 2013). However, the changing pattern was highly correlated by the longitude of the station (Figure 4). Higher depth values in the future were observed for the long-duration events in Payson and Phoenix-Encanto, which are distributed near the central part of Arizona (Figures 5 and 6). Generally, precipitation is the lowest in southwest Arizona and increasing toward the northeast, with the maximum values happening near the central parts of the state (Karnieli 1990; Brown & Comrie 2002), where Payson and Phoenix-Encanto stations are located. The updated curves predicted lower depth values almost for all the durations and stations from 2010 to 2015 compared to NOAA Atlas 14 based on analyzing the temporal dynamics of DDF curves. This period can be considered as the driest one compared to other time periods (Figure 4).

Other factors could also result in the variation of simulated and published PPFEs and DDF curves, such as different applied methodology, probability distribution, and the fitted POR. In NOAA Atlas 14, the Parameter-elevation Regression on Independent Slopes Model (PRISM) was adapted to spatially interpolate mean annual maximum precipitation to a uniformly spaced grid (Bonnin et al. 2011). Since the availability of station data in the calculation of PPFE for Arizona was limited, the PPFEs were based on the application of regional analysis using L-moments (Hosking & Wallis 1997). Although the L-moments are considered to be a robust estimator for the marginal moments of a time series, they do not take into account the auto-correlation structure of the process. This issue could have been avoided if the K-moments were to be used where this cross-dependency between the marginal distribution and the auto-correlation function is taken into account (Koutsoyiannis 2020). The regional analysis is used to extend the at-stations POR in order to improve the estimation of the shape and scale parameters. The performance of the traditional interpolation technique employed by PRISM could be affected by the availability of observation data (Gutmann et al. 2012), specifically in Arizona, where the availability of observation data throughout the state is limited. However, WRF models employ dynamically downscaled global climate projections to provide gridded precipitation data with a high temporal resolution (in this case, 1 min) for the desired timeframe. Additionally, the WRF model accounts for geographic information such as topographic channeling of flow and local frontal systems that play a significant role in precipitation distribution (Gutmann et al. 2012).

Moreover, the GEV probability distribution was applied in NOAA Atlas 14 as the most appropriate for precipitation frequency calculations based on the goodness-of-fit to the precipitation data and the climatologic and geographic consistency considerations, while Weibull distribution was adapted in this study for frequency calculations. Additionally, the assumption of a 6-year time duration to showcase the climate change effects on precipitation events in 20-year increments versus the 20-year POR can affect the depth and frequency estimates resulting from simulated data.

The impact of changes in PPFEs for urban drainage systems is unclear without further analysis of different emission scenarios explained in IPCC. Although the spatial dependence of precipitation events could increase the uncertainty associated with future estimates of precipitation frequency and depth, rainfall events of different durations appeared to be affected by climate change regardless of location or time interval. Therefore, the use of historical PPFEs could result in applying overestimated precipitation frequency in the design of water infrastructure, which could lead to excessive sediment deposition (Fan 2004; Wilson 2011; Tang 2016; Song et al. 2018). Oversized and mildly sloping drainage systems run the risk of having a flow velocity less than the particle-settling velocity. The accumulation of sedimentation in drainage systems can result in a loss of flow-carrying capacity that may restrict flow and cause upstream surcharge, local flooding, and enhanced solids deposition (Fan 2004).

Another risk associated with significant sediment deposition in stormwater drainage systems is the potential to shock load the receiving waters. The loss of carrying capacity of drainage systems combined with an extreme precipitation event can provide a flow velocity capable of re-suspending the deposited sediment and transporting it to the receiving waters (Stein et al. 1999). The high-suspended solids concentration inflow can be detrimental to the receiving waterbody (Fan 2004).

Further, it is worth mentioning that stormwater drainage systems are inherently complex and extended systems with a design life, reaching almost 100 years (Mailhot & Duchesne 2010). The modification of urban drainage to adapt for changes in precipitation depths and distributions is going to take time, and the adaptation strategy developed will need to be periodically revisited in the context of an uncertain future (Griggs & Noguer 2002). The model used to project PPFEs and DDF curves will need to evaluate numerous likely emission scenarios with a more extended range of data (50 years or more) to decrease the uncertainties related to modeling and the underlying assumptions.

Redefining the design criteria of urban drainage infrastructure will require careful evaluation of the expected design life, the annual exceedance probability (or average recurrence interval), and duration accepted. The long-term impact of climate change on PPFE will likely require new infrastructure to be oversized for a fraction of its design life (Mailhot & Duchesne 2010). This oversizing could cause adverse effects due to excessive deposition, as discussed before. The oversizing of urban drainage systems will require routine cleaning to avoid the negative impacts associated with excessive sediment deposition.

This study intended to evaluate the stationary climate assumption in the production of PPFEs in NOAA Atlas 14. A WRF model was developed to generate precipitation data based on the medium emission scenario for the following timeframes: 2010–2015, 2030–2035, 2050–2055, and 2070–2075. Based on the comparison between the simulated PPFEs and NOAA-published values, there seems to be a large variability and uncertainty in the estimation of precipitation extremes. On average, the precipitation depth values were decreasing for shorter-duration precipitations. However, the wettest areas were less sensitive to climate change compared to the driest part of the state, which experienced considerably lower depths for all storm durations. Various time range was defined in the model to showcase the climate change impacts on precipitation events in 20-year increments. Years 2010–2015 could be considered as a dry period where the discrepancies between the historical and modeled PPFEs were the highest. The period of 2070–2075 was the wettest, with the DDF curves presenting a similar trend and values to the historical ones in the NOAA. This conclusion refers to the fact that the annual precipitation seems to be decreasing in most parts of the west and southwest of the United States, even though there are significant regional differences in precipitation patterns.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.