An overall optimization and solution framework for urban historical and future DRIF under climate change Open Access

As one of the most sensitive and vulnerable areas affected by climate change, extreme weather and disaster events occur frequently in China. Geological disasters such as urban waterlogging, farmland waterlogging, mountain torrents and mudslides caused by heavy rainfall have caused serious influence to people's production and life. With regard to the future global climate change, the fifth assessment report of the IPCC (Intergovernmental Panel on Climate Change) indicated that evolutions in precipitation and melting of snow and ice in many regions are changing the hydrological system, affecting the quantity and quality of water resources (medium confidence) (Pachauri et al. 2014). From 1880 to 2012, the global average temperature increased by 0.85 °C (Zhang et al. 2021). Under the influence of climate change, glaciers close to the global scale are shrinking continuously (high confidence), affecting the downstream runoff and water resources (medium confidence). Climate change is causing warming and thawing of permafrost at high latitudes and altitudes (high confidence). The extreme weather and climate events caused by the impact of human activities include the decrease of low-temperature extreme events, the increase of high-temperature extreme events, the increase of extremely high sea levels and the increase of the number of heavy precipitation events in some regions. The number of heavy precipitation events may show an increasing trend in more land areas than in areas with decreased precipitation. All the emission scenarios evaluated predict that the surface temperature will increase in the 21st century. The frequency of heat wave is likely to be higher and longer, and the intensity and frequency of extreme precipitation in many areas will increase. The oceans will continue to heat up and acidify, and the global average sea level will also continue to rise. The results of the above analysis show that the probability and frequency of extreme weather will have a significant increasing trend for a long period in the future, which brings great challenges to countries around the world to deal with disastrous weather.

The design rainstorm intensity formula (DRIF) is an important basis for reflecting the regularity of rainfall, guiding the design of urban drainage and waterlogging prevention projects and the construction of related facilities. The accuracy of its calculation has an important impact on the project cost, drainage capacity, people's property and life safety (Mirhosseini et al. 2013; Fadhel et al. 2017). The rainstorm intensity of short-duration heavy precipitation becomes larger under climate change. Compared with the original DRIF, the precipitation for the same return period and the duration will increase. The original DRIF can no longer meet the design requirements of future projects.

The DRIF is obtained by fitting the intensity-duration-frequency (IDF) curves. The IDF assumes that the historical precipitation is stationary, so it can be used to characterize future changes in extreme precipitation. When climatic conditions change rapidly, IDF curves generated only based on historical data will erroneously represent future climate changes (Mailhot et al. 2007; Srivastav et al. 2014). In urban hydrology, the design of drainage infrastructure has traditionally been based on statistical analysis of historical precipitation records. It is assumed that the intensity and frequency of past events are statistically representative of what may happen in the near future. However, in the context of climate change, the assumption must be revisited and design criteria for drainage infrastructure revised to account for expected changes in the intensity and frequency of heavy rainfall. Otherwise, it will have a serious impact on the future design, operation and maintenance of drainage facilities (Mailhot et al. 2007; Tfwala et al. 2017).

Considering the increase of observed heavy rainfall events, IDF curves should be updated to take into account the changing climate, especially for the design of urban infrastructure. Beijing is the capital of China, with complex terrain, high population concentration, and severe rainstorms (Lu & Cui 2022). For example, in 2012 and 2016, Beijing experienced torrential rains on July 20th and July 21st, respectively, which caused severe urban waterlogging and massive economic losses. To the best of our knowledge, there are currently no studies evaluating the impact of climate change on the IDF curve in Beijing. Few studies have calculated future DRIF based on IDF curves under climate change. At present, the research on IDF curve mainly has the following problems: (1) When the extreme value of sub-daily precipitation in one year is much higher than that in other years, it will lead to poor fitting effect of the tail of the frequency distribution curve. Overfitting to the tail will result in smaller rainstorm intensity for the fitting of head; (2) For the calculation of future IDF, many studies have modified the historical IDF curve to obtain the future IDF according to the proportion of rainstorm intensity values in different return periods of future and historical climate models. This solution method using proportional increments does not consider the correlation between future IDF and observed values.

The objectives of this study are as follows: (1) Based on the minute-level precipitation data, the influence of this extreme value on the shape and fitting accuracy of the historical frequency distribution curve is analyzed: (2) The calculation process of future IDF curve based on the Coupled Model Intercomparison Project Phase 6 (CMIP6) is given; (3) The variation of the equivalent return period before and after climate change is analyzed; (4) The historical and future DRIFs of Beijing under climate change are calculated. This research can provide a reference for the planning and design of water resources projects in Beijing or other areas in the future.

The minute-level precipitation for the observation comes from the Guanxiangtai Weather Station in Beijing, the range of which is from 1951 to 2012. CMIP6 was adopted as projected precipitation for future climate change analysis. Its range includes historical periods (1951–2012) and future periods (2023–2100). In this study 12 models were selected for accuracy evaluation, which are listed in Table 1. The results of this evaluation show that the accuracy of CanESM5 is the best. The PBIAS and RMSE of CanESM5 are the smallest for annual precipitation, which are −7.33 and 239.66, respectively. The trend of the change is more consistent with the observed value. The CC of CanESM5 and the observed monthly precipitation is 0.59. CanESM5 is chosen to describe changes in precipitation for future period. The process of evaluation and selection for CMIP6 is not described in detail in this study. For the method of this evaluation, refer to Dong & Dong (2021); and Guo et al. (2021). The latest CMIP6 combines RCPs and shared socioeconomic pathways (SSPS) to produce a more reasonable future scenario (Peng & Li 2021; Su et al. 2021). Four combinations of SSP-RCP scenarios are adopted in this study,which are SSP1-2.6 (sustainability), SSP2-4.5 (middle of the road), SSP3-7.0 (regional rivalry), SSP5-8.5 (fossil fuel development) (Su et al. 2021). ‘Nominal Resolution’ is 100 and 250 km. ‘Variant Label’ is r1i1p1f1. They describe possible future worlds and represent different combinations of mitigation and adaptation challenges. SSP1-2.6 has significant changes for landuse,and cover ratio of its global forest will increase. The landuse and aerosol pathways of SSP2-4.5 are not extreme relative to other SSPs. It combines intermediate societal vulnerability with an intermediate forcing level. SSP3-7.0 is a scenario with a large number of land use changes (especially decreased global forest cover) and high NTCF emissions (especially SO2). SSP3-7.0 is a scenario with both substantial land use change. SSP5-8.5 is the only SSP scenario with emissions high enough to produce a radiative forcing of 8.5 Wm⁻² in 2100.

Outliers are extreme values that stand out from the distribution of the data in the graph or table. Outliers in the series will affect the statistical results of samples, such as mean, standard deviation, coefficient of variation, skewness coefficient and kurtosis coefficient, thus affecting the parameters of the distribution function. Seven detection methods of outliers were adopted in this study, which are Z-Score (ZS), 3-sigma (3σ), Modified Z-Score (MZS), median absolute deviation (MADe),Box-Plot Method (BPM), Grupps-Beck Test (GBT), Stedinger Test (ST) (Grubbs 1969; Crosby 1994; Kannan et al. 2015; Asikoglu 2017; Lu et al. 2018).

• (1)
  The Z-score assumes that sample x follows normal distribution, which is ⁠. The Z-score is calculated as follows:

  

  (1)

  where is the sample X; is mean of sample X; variance of sample X. When Z_i > γ, x_i is the outlier. For the γ, the thresholds used in Asikoglu (2017) and Kannan et al. (2015) are 2.5 and 3, respectively. The γ can be adjusted according to the distribution range of abnormal values of sample data, which is 2.5 in this study. The Z-score is greatly affected by outliers and is not suitable for the detection of outliers in small samples.
• (2)
  The 3-sigma is also known as the 68-95-99.7 rule. It is similar to the calculation principle of the Z-Score. When the γ of the Z-Score is equal to 3, it is the 3-sigma. It assumes that the sample follows a normal distribution, and when the positive or negative moments of the sample data are outside σ, 2σ or 3σ, the data is considered to be an outlier (Figure 3). The probabilities of x falling within different confidence intervals are as follows:

  

  (2)

  

  (3)

  

  (4)
• (3)
  The modified Z-score adopts median and median absolute deviation instead of mean and standard deviation, which solves the problem of poor detection ability of Z-score for outliers in small sample data (Kannan et al. 2015). When > 3.5, there are outliers in the sample data. The formula is as follows:

  

  (5)

  

  (6)

  where is the sample X; is the number of data.
• (4)
  MADe is a stable detection method that is less affected by extreme values. The study adopts the 3MADe method, as follows:

  

  (7)

  

  (8)

  

  (9)

  where is the outlier; is the sample X.
• (5)
  The boxplot obtains the interquartile range IQR by calculating the difference between the upper quartile Q3 and the lower quartile Q1 of the sample, and which is suitable for dealing with symmetric and skewed data. In the study, IQR15 is adopted as the judgment criterion for outliers (Saleem et al. 2021):

  

  (10)

  

  (11)

  

  (12)

  

  (13)

For extreme values, this study focuses on outliers that have a greater effect on the higher return period of the frequency curve. When the series of the sample is judged to have an extreme value, the extra-large value will be processed in three ways: (1) Extracting the second maximum value of the daily precipitation in the ith year(SMV); (2) Extract the average value of the daily precipitation in the ith year that is greater than the critical value(AVG); (3) Extracting the maximum value of the daily precipitation in the ith year which is less than the critical value (MVL). First, the maximum daily precipitation for each year is extracted to form a statistical sample. When it is judged that there is an extreme value in this sample, the data of the year in which the extreme value is located will be processed. The data of this year are sorted in descending order, and the second largest value is extracted and used to replace the ith extreme value in the sample. This processing method is the principle of SMV. The calculation methods of AVG and MVL are similar to SMV. They extract respectively the average value of daily data greater than the critical value and the maximum value of daily data less than the critical value in the year where the extreme value is located.

Step 1: The climate model was corrected using a quantile delta mapping method of daily precipitation based on frequency (DFQDM). The method can correct the frequency and value of wet days of precipitation. This 95th percentile was used to divide monthly precipitation into normal and extreme precipitation. A mixed Gamma distribution was used to fit the cumulative distribution frequency(CDF) of the two parts of the daily precipitation. A quantile delta map was used to correct for precipitation.

Step 2: is the maximum value of the extracted observations in the jth duration of the ith year. The number of durations of the precipitation in this study contains 15, which are 5, 10, 15, 20, 30, 45, 60, 90, 120, 150, 180, 240, 360, 720 and 1440. This duration takes into account both short and long-duration precipitation. The simulation of urban drainage system is more concerned with the simulation of short-duration precipitation process.

Step 3: and are the maximum daily precipitation (1440 min) of the ith year in the historical and future periods of the extracted model, respectively;

GEV, Gumbel, Pearson3, Exponential are adopted to compare and obtain the best distribution function of empirical frequency points in this study. In the process of parameter optimization of the frequency distribution curve and DRIF, this study adopts several widely used solution methods. They are respectively L-moment (Kotz & Nadarajah 2000), Multistart search least squares algorithm(MSLSA) (Ugray et al. 2007), Non-dominated sorting genetic algorithm II(NSGAII) (Deb et al. 2002), Multi-objective particle swarm optimization (MOPSO) (Coello & Lechuga 2002), Multi-objective evolutionary algorithm based on decomposition(MOEA/D) (Zhang & Li 2007), Nondominated sorting and local search (NSLS) (Chen et al. 2014). The formula for the error of the fitting accuracy for DRIF is as follows:

• (1)
  Mean absolute root mean square error (MARMSE):

  

  (29)
• (2)
  Mean relative root mean square error (MRRMSE):

  

  (30)

  where is rainstorm intensity fitted to theoretical frequency curve (⁠⁠) or DRIF (⁠⁠), mm/min; is the rainstorm intensity fitted by empirical frequency curve (⁠⁠) or theoretical frequency curve (⁠⁠), mm/min; n is the number of sample data.

• The selection of the distribution function is based on the performance of MARMSE and MRRMSE for the DRIF. In China, the regulations promulgated by the Ministry of Housing and Urban-Rural Development (MOHURD) require that when the return period is two to 20 years, MARMSE should not be greater than 0.05 mm/min in areas with general rainstorm intensity, and MRRMSE should not be greater than 5% in areas with large rainstorm intensity. The accuracy of this calculation was evaluated using two- to 100-years when the DRIF was fitted. The interval of this return period is taken as 0.5 years. The accuracy requirements for MRRMSE and MARMSE are for DRIF rather than IDF in this study. These two indicators can be used as an evaluation index for the accuracy of the fitting of the IDF curve.

It can be seen from Table 2 that the solution method and frequency distribution of DRIF to meet the requirements of the fitting accuracy are the Exponential optimized by NSGA-II, the Gumbel optimized by MOEA/D, and the Gumbel distribution optimized by NSLS. The optimized accuracy of NSLS is better than that of NSGA-II, NSGA-III and MOEA/D, and the Gumbel has the best performance. The MARMSE and MRRMSE of NSLS were 0.045 and 3.475%, respectively. The MRRMSE is the lowest among all distributions. The coefficient of determination (R²) is close to 1. In this 100-year return period and different durations, the intensity of the rainstorm varies from 0.169 to 4.124 mm/min. The MARMSE and MRRMSE of the L-moment do not meet the accuracy requirements. L-moment was established by Hosking in 1990 and is widely used in fields such as hydrology and civil engineering (Wan Zin et al. 2009). Compared with other methods of moments, the main advantage of the L-moment method is that it is less affected by the variability of the samples (Bílková 2012). The L moment is more stable and the solution speed is faster, and it can provide safer results in the case of small samples. When the IDF table needs to be obtained quickly, the L-moment method is an optimal choice. MOPSO has the advantages of easy implementation and fast search. The main disadvantage of this algorithm is that the non-dominated solutions of uniform distribution in the solution space are poorer and the diversity of solutions is insufficient (Zhang et al. 2016). NSGA-II and MOEA/D are the most widely used multi-objective optimization algorithms. The disadvantage of NSGA-II is that it performs poorly on high-dimensional problems when the number of objective functions exceeds three (Zhao & Li 2014; Huang et al. 2019). The probability of crossover and mutation of NSGA-II is not perfect, and its process of calculation is very time-consuming. MOEA/D has problems such as reduced population evolution efficiency and poor evolution quality in the calculation process. It may be due to the above shortcomings that the performance of these algorithms is not as well as that of NSLS. NSLS is able to find better distributions of solutions and better convergence to true pareto optimal frontiers (Chen et al. 2014). The Gumbel optimized for the NSLS will be adopted to compute the historical DRIF in this study.

Table 2

View Large

After the outliers are processed, the rainstorm intensity is slightly reduced in different return periods compared with before processing. In Table 4, the color of dark red represents a larger proportion of reduction, and the color of dark green represents a smaller proportion of reduction. The proportion of this reduction increases as the return period increases. The scale of this reduction varies from −0.1 to −3.1%, with a slight decrease in value. The Z-score method did not change the value of the sample points significantly, and the ratio was the largest at 60-min.

Table 4

View Large

The calculation results of the DRIF for the historical period are as follows:

(1) No handling of outliers	(2) Handling of outliers

(1) No handling of outliers	(2) Handling of outliers

View Large

The calculation results of DRIF from 2023 to 2100 are as follows:

SSP1-2.6:	SSP2-4.5:
SSP3-7.0:	SSP5-8.5:

SSP1-2.6:	SSP2-4.5:
SSP3-7.0:	SSP5-8.5:

View Large

Optimization results of parameters of DRIF from 2023 to 2100 are shown in Table 6.

At present, there are few relevant studies on the influence of extreme values on the IDF curve. A variety of judgment methods and replacement schemes for outlier are adopted in this study, and the purpose is to propose a more reasonable method for generating IDF curves. This combined method can prevent IDF curve from overestimating the rainstorm intensity in the high return period. In this study, this method is only compared and verified in a weather station, so it may have some uncertainty and still needs further verification. This study proposes a suitable method for replacement of the extreme value for the IDF. For the solution of IDF curve in the future, the previous research mainly includes four methods: (1) After the climate model is corrected, the frequency distribution of historical and future precipitation is fitted. It is assumed that the ratio of rainstorm intensity on daily and sub-daily of observation under different return periods does not change in the future. The daily IDF of the future period is corrected by the ratio to obtain the IDFs of the sub-daily (Tousi et al. 2021); (2) This second method is similar to the first method. Daily frequency curves for model and observation are calculated. The ratio of daily precipitation under different return periods of model for the historical and future period is calculated. This ratio is used to obtain the future IDF by correcting the frequency curve of the observed sub-daily. This method is a simplified processing method (Zhou et al. 2018); (3) First, the frequency distribution of historical periods is calculated. Quantile mapping is adopted to establish the non-linear relationship between the rainstorm intensity of modeled daily and the observed sub-daily in the historical period. Estimates of rainstorm intensity for different return periods and durations of the model over the historical period were calculated. Then, adopting the same method as (2), the estimated value of the modeled IDF of in historical period is modified to obtain the IDF of the future period. Srivastav et al. (2014) adopts this method to calculate the IDF for the future period. However, this study has many shortcomings. It does not adopt equidistant quantile mapping (EQM). The non-linear relationship of the IDF curve established by this study is the relationship between the modeled future daily and the observed sub-daily. This relationship cannot express the relationship between the IDF of future daily and sub-daily; (4) High-resolution precipitation for this future period is generated based on the random weather generator (Mirhosseini et al. 2013; Doi & Kim 2021). The method has great uncertainty for the generation of precipitation. The IDF tables for future periods generated by the models of different weather generators may be different.

The difference between this study and others is that the precipitation of the model is corrected to establish a non-linear relationship between the model and the observation in different sub-dailys of the historical period. The non-linear relationship of the maximum precipitation in different years for the historical period is also applicable in the future period, especially for the near and medium term of the future period. The daily CDF of this future period is calculated using the non-linear relationship established by the historical period to obtain the IDF table of the sub-daily. Therefore, the method established by this study should be more reasonable. In Figure 10, the IDF of different scenarios in this future period are all larger than the observation in the historical period. However, Tousi et al. (2021) indicated that the rainstorm intensity of SSP2 of MPI-ESM1-2-HR for the partial return period of 1440-min is smaller than the observation. Therefore, it may be unreliable to directly use the corrected climate model to calculate the IDF table for future periods. The choice of the climate model should first be evaluated for accuracy in the historical period. The model with the best performance of this accuracy can then be used to calculate the IDF table for the future period. The study finally calculates three types of DRIFs, which are respectively: (1) historical periods without outlier processing, (2) historical periods with outlier processing, (3) future periods with different SSP scenarios. The specification of MOHURD does not mention a requirement for the accuracy of the fit of the DRIF for future periods. According to the requirements of the accuracy of this observation, this future DRIF does not meet the conditions. However, the error of the future DRIF calculated by this study is very low, which can meet the design requirements of the infrastructure.

This research has established an overall framework for calculating DRIF, which includes historical and future periods. This DRIF is obtained based on the fitting of IDF table. Therefore, this research focuses on the generation of the IDF curve. In the historical period, the judgment of extreme value for the statistical sample and its influence on IDF curve are analyzed. Based on the non-linear relationship for different quantile values of between the model and the observation in the historical period, the study establishes a solution method for generation of IDF curve in the future period. The conclusions drawn from the study are as follows:

• (i)

  The optimization accuracy of the parameters of the DRIF by the NSLS algorithm is the best. In the historical period of Beijing, the variation range of the rainstorm intensity of different durations for the 100-year return period is 0.169–4.124 mm/min;

• (ii)

  For the processing of this outlier, the Z-Score was recommended to judge the outlier in this study. The best alternative to the outlier is the AVG method, which takes the average of the series in the sample that are larger than the critical value. The error of SMV in the high return period is still larger. The rainstorm intensity of MVL in the high return period is obviously smaller. After the outliers are processed, the value of the IDF table is slightly reduced. The percentage of reduction for IDF varies from −0.1 to −3.1%.

• (iii)

  There is a perfect non-linear relationship between the model and the observation at different quantiles over the historical period. The projected IDF under the four scenarios for this future period are all larger than the observation.

• (iv)

  Except for the two- and three-year of observation in SSP2-4.5, the return periods of observation in the historical period are all larger than equivalent return periods in future period.

Xingchen Ding: Methodology, Formal analysis, Resources, Writing – original draft, Data curation. Weihong Liao: Writing – review. Hao Wang: Writing – review. Hao Wang: Conceptualization. Xiaohui Lei: Supervision. Jiali Yang: Supervision, encouragement.

This study was supported by the National Natural Science Foundation of China (No. 52179027).

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.