Calculation method of short-duration rainstorm intensity formula considering nonstationarity of rainfall series: impacts on the simulation of urban drainage system

In recent years, the impact on global urban waterlogging has been intensified, and the frequency and degree of urban flood is increasing (Douglas et al. 2008; Kundzewicz et al. 2014). The urban drainage system is the main engineering system of urban sewage and rainwater transportation and drainage. It is the infrastructure to protect the city from flood damage. There are many impervious areas on the urban surface, and the natural rivers and lakes have been transformed. On the one hand, the rainfall interception and infiltration amount are significantly reduced, and the net rainfall is increased; on the other hand, the convergence speed is obviously accelerated. Once the rainstorm exceeds the design standard of drainage system, waterlogging disaster is likely to occur in the city, causing incalculable loss and harm (Zheng et al. 2016; Anker et al. 2019; Abass et al. 2020). At present, the design of urban drainage system depends on the intensity and frequency of historical extreme heavy rainfall, which is usually summarized as intensity-duration-frequency (IDF) curves or further generalized into the form of rainstorm intensity formula. The existing infrastructure is generally designed by the static rainstorm intensity formula. It is assumed that the environment formed by rainstorm does not change and the probability density of hydrological variables does not change with time, which shows that the statistical characteristics do not change (Strupczewski et al. 2001; Milly et al. 2008). However, the fifth assessment report of the Intergovernmental Panel on Climate Change of the United Nations (IPCC AR5) pointed out that the intensity and frequency of rainstorm in most land areas of the world may increase due to the impact of global climate change and land-use change. This means that the existing hydrological steady-state law has been broken, and the rainstorm intensity formula relying on the assumption of stationarity may underestimate the extreme rainfall events, which will bring unnecessary losses if applied to the design of urban drainage systems (Khaliq et al. 2006; Willems et al. 2012; Forzieri et al. 2018). In addition, some studies have shown that with the increase of research time, the deviation of nonstationary and stationary design values will be larger and larger (Agilan & Umamahesh 2016; Lima et al. 2018; Hosseinzadehtalaei et al. 2020). Therefore, it is necessary to explore the calculation method of nonstationary rainstorm intensity formula and its influence on urban drainage system simulation so as to provide a theoretical basis for urban drainage system design.

Since 1901, the global average temperature has increased by 0.89 °C, which has seriously affected the precipitation pattern, and this trend has become more obvious in recent years (Banerjee et al. 2021). The rainfall patterns of many cities in the world have been studied extensively (Mishra et al. 2015; Rahman et al. 2017; Yurekli 2020). Although the short-duration rainstorm has an obvious change trend, it is not necessarily reflected in the long-term rainstorm (Madsen et al. 2009; Fujibe 2013; Zhou et al. 2017). It is worth mentioning that the short-duration rainstorm is one of the main causes of urban waterlogging. There are a lot of impervious surfaces on the urban surface, less rainfall interception and infiltration, and the urban surface roughness is small. The natural river channel is transformed into an artificial river channel with fast convergence speed (Chen et al. 2015; Bertilsson et al. 2019; Pour et al. 2020). It is of great significance to calculate the nonstationary short-duration rainstorm intensity formula and comprehensively evaluate the drainage system performance for the planning and reconstruction of urban drainage infrastructure in the future.

Many methods have been used to derive the return period values of nonstationary hydrological time series (Coulibaly & Baldwin 2005; Cheng & AghaKouchak 2014; Salles et al. 2019). GAMLSS is a widely used method to establish the relationship between the distribution function of time series and the explanatory variables (Salas & Obeysekera 2014; Balistrocchi & Grossi 2020). The GAMLSS model has become an important tool for hydrological extreme value time series analysis because of its good flexibility and adaptability (Villarini et al. 2009; Villarini & Serinaldi 2012; Hao et al. 2019). It has shown excellent performance in the frequency analysis of nonstationary flood series (Li & Tan 2015; Debele et al. 2017; Zhang et al. 2018), but there is little research on its application to the short-duration rainstorm intensity formula.

Recent studies on the impact of environmental change on extreme rainstorm and urban drainage systems show that the risk of damage to the existing drainage system infrastructure caused by the design rainstorm calculated under nonstationary conditions will increase, and the change of short-duration rainstorm will have an important impact on the outflow, overload and overflow of urban drainage systems (Denault et al. 2006; Huong & Pathirana 2013; Andimuthu et al. 2019; Yang et al. 2020). Unfortunately, most studies have not carried out a full and detailed impact analysis of urban drainage systems, such as the flooded area, the maximum flooded depth, the risk-bearing pipeline and node and the detailed inundation process. The two-dimensional MIKE model can remedy the defects of the one-dimensional model and simulate more detailed flood processes (Bisht et al. 2016).

The purpose of this study is to consider the change of short-duration rainstorm intensity formula under nonstationary state and provide a detailed performance evaluation of urban drainage systems in the study area. It can provide possible short-duration extreme rainstorm scenarios in the study area in the future and provide support for the reconstruction and optimization of existing drainage infrastructure. The design value of the nonstationary return period is calculated by the GAMLSS method, and the short-duration rainstorm intensity formula of different periods (1981 represents the past, 2010 represents the present and 2030 represents the future) is fitted, which is compared with the traditional stationary rainstorm intensity formula. The comprehensive evaluation of the drainage system is realized by dynamically coupling the urban water simulation model and the MIKE 21 two-dimensional surface flow model.

Zhengzhou University campus, located at 113°39″ E and 34°47″ N, is located in the continental monsoon climate zone of north temperate zone, with four distinct seasons. The annual average precipitation is between 600 and 700 mm, and the rainfall from June to September accounts for more than 50% of the annual rainfall. The local urban floods occur frequently in recent years, which may be caused by climate change and rapid urbanization (Wang et al. 2020). The original fixed rainstorm intensity formula may not be appropriate. How to consider the nonstationarity has become an important problem. The fixed rainstorm intensity formula of the 2000s is often used in local urban planning, and it is vital to develop a rainstorm intensity formula calculation method which can consider the nonstationarity.

This study uses the observed IDF data of Zhengzhou meteorological station from 1981 to 2010, including 11 annual maximum rainfall series with different rainfall durations (5, 10, 15, 20, 30, 45, 60, 90, 120, 150 and 180 min) and the actual observed rainfall data on 12 August 2017 (Figure 1). Zhengzhou meteorological station is very close to the study area, and its rainfall characteristics can represent the characteristics of rainstorm in the study area. In addition, the basic data used in the study include the urban drainage network planning map of shp format from Zhengzhou University logistics office, including 369 pipelines and a drainage outlet (Figure 2(c)). The digital elevation map (DEM) is obtained by Kriging interpolation based on discrete elevation points (Figure 2(b)). Before interpolation, we correct the elevation of some abnormal points. According to the remote sensing image map (Figure 2(a)), the land-use types can be divided into four types: building, road, greenbelt and lake (Figure 2(d)).

First, the stability test is used to verify whether the rainfall sequence is nonstationary. If the rainfall sequence is stable, the design value of recurrence period is calculated by the traditional frequency analysis method, and the formula of rainstorm intensity is fitted. If not, the GAMLSS method is used to calculate the recurrence period design value and the rainstorm intensity formula.

The stability test of rainstorm data mainly includes trend test, abrupt test and periodic test. This paper uses the Mann–Kendall trend test, sliding t-test, ordered clustering method and wavelet analysis method to verify the nonstationarity of rainfall sequence. The variable UF and the variable UK exceed the critical range at the significance level of 0.05, and there are intersection points between the UF curve and the UK curve in the critical range (Figure 3), which means that all rainfall sequences may have trend variation and mutation. At the significance level of 0.05, sliding t-test and ordered clustering were used to detect significant mutation points in all rainfall duration rainstorm sequences (Table 1). The wavelet real-part isoline shows that the maximum rainfall intensity series of 5–90 min of rainfall has many abundant and dry oscillations in the time scale of 8–16 years. The maximum rainfall intensity sequence of 120–180 min oscillates obviously between the time scale of 3–8 years. According to the wavelet error map, the possible periods are 4, 7 and 11 years (Figure 4). The trend, mutation and periodicity may exist, which indicates that the stability of rainfall series has been destroyed.

GAMLSS is a kind of (semi) parametric regression model (Rigby & Stasinopoulos 2005). It is a generalized additive model based on location, scale and shape parameters developed by Rigby and Stasinopoulos in order to overcome the limitations of correlation models such as generalized linear model and generalized additive model. GAMLSS allows us to describe the linear and nonlinear relationship between statistical parameters and covariates, which give GAMLSS obvious advantages in dealing with nonuniform frequency analysis.

In this study, the position parameters and scale parameters corresponding to the mean and variance are analyzed, and the covariant is time t. There may be large error if only the linear trend is used to simulate nonstationary series (Agilan & Umamahesh 2017). In this paper, we also use the maximum polynomial of degree 3 to study the distribution of covariates. Four commonly used two-parameter distribution functions are selected: Gumbel (GU), Gamma (GA), Lognormal (LOGNO) and Weibull (WEI). On this basis, seven models are set up: (a) and are both invariant; (b) is invariant, is a linear function of time; (c) is a linear function of time, is invariant; (d) and are both linear functions of time; (e) is invariant, is polynomial function of time; (f) is a polynomial function of time, is invariant and (g) and are both polynomial functions of time. The Akaike information criterion (AIC) is used to select the optimal model.

The changing nature of the Earth's climate and rapid urbanization lead to the change of rainfall characteristics in urban areas, and the stability of rainfall series is destroyed; the failure of fixed assumptions in urban drainage system design may threaten the current drainage facilities. Therefore, nonstationarity should be considered in the design of urban drainage systems. The GAMLSS method is used to calculate the design value of return period under nonstationary state, the Marquardt method is used to fit the short-duration rainstorm intensity formula, and the short-duration rainstorm intensity formula of three periods under nonstationary state is calculated so as to analyze the difference with the commonly used rainstorm intensity formula under stationary state.

After considering the influence of nonstationary factors on the short-term rainstorm intensity formula, the influence of nonstationary factors on urban drainage infrastructure planning and design is analyzed by urban rainfall–runoff simulation. The models of urban hydrology MIKE URBAN and two-dimensional surface flow model MIKE 21 are coupled by MIKE flood platform. In terms of rainfall–runoff model, MIKE URBAN uses the time–area curve method to calculate. In terms of the pipe network hydrodynamic model, MIKE URBAN adopts Saint Venant equations (continuity equation and momentum equation) of one-dimensional free surface flow and uses the Abbott–Ionescu six-point implicit scheme finite difference method to solve the problem. In the aspect of two-dimensional surface flow simulation, MIKE 21 uses the implicit finite difference method to solve the two-dimensional unsteady flow equations.

The DEM data are transformed into the dfs2 format that can be used, and the buildings and roads are superimposed to represent the complex terrain of urban areas. The runoff flows into the drainage system through the connected nodes. When the drainage capacity of the drainage system is exceeded, the flood will overflow to the surface and form submergence.

The model can be set to operate under various boundary conditions (such as rainfall–runoff and external input flow). In this model, only the rainfall time series (dfs0 format) is used as the boundary condition to analyze the operation of the drainage system under different short-duration design rainfall conditions.

According to the AIC, the variation of annual maximum rainfall intensity sequence with rainfall durations of 5 and 90 min can be simulated by GU, although other sequences using WEI are better. At the same time, the model of and changing linearly with time performs well in 5 and 90 min series, while the model of changing linearly with time is more suitable for other series. This fact shows that although we consider the influence of nonlinear link function, for the data of this study, the linear function performs well and can achieve good simulation effect.

In the residual normal worm chart (Figure 5) of the optimal models, the scatter points are distributed in the 95% confidence interval of the confidence level and close to the red curve. This shows that the model fitting effect is good.

According to the fitting results of the GAMLSS model, the quantile curves are summarized and predicted, and the design value curves of different return periods are obtained (Figure 6). It can be seen from the figure that although the change range is not consistent, the annual maximum rainfall intensity series change is mainly increased, and this trend may be more obvious in the future.

We set up four scenarios to analyze the development of the short-duration rainstorm intensity formula under the nonstationary assumption and the difference between the formula under the nonstationary assumption and the formula under the traditional stationary assumption. For convenience, we use to represent (a) the scenario of nonstationarity in the past 1981, (b) the current situation of nonstationarity in 2010, (c) the scenario of nonstationarity in the future 2030 and (d) the scenario of stationarity assumption. The short-duration rainstorm intensity formula suitable for the return periods of 2–50 years is obtained by means of Marquardt method (Table 2), and the goodness of fit was evaluated by mean absolute error (MAE). (In the Chinese industry standard, MAE should not be greater than 0.05 mm/min.)

The input data of rainfall used in the experiment is from an actual rainfall in 2017 (Figure 1). In addition, there is typical design rainfall calculated by the Chicago method based on the short-duration rainstorm intensity formula of the above four scenarios (the rainfall duration is 3 h and the rainfall peak coefficient is 0.5; Figure 7). The mean surface velocity of the flood simulation model is 0.3 m/s, the reduction factor is 0.9 and the initial loss is 0.0006 m. Imperviousness of different land-use types is shown in Table 3.

It is difficult to obtain the actual measured data of inundation area and time. Due to the lack of actual observation data of flood inundation, a qualitative verification method is adopted to verify the model. It is indeed found that there was a waterlogging event on 12 August 2017, and the corresponding rainfall data with a 70-min observation interval of 10 min from 14:30 to 15:40, with a total precipitation of 51 mm (close to the once-in-2-year rainfall). We simulated the event and compared it with the actual observation. The simulation results (Figure 8 and Table 3) show that the inundated area is 199,275 m², accounting for 6.86% of the total area of the study area (the study area in the model is 2,903,625 m²). Among them, the areas with serious inundation mainly include the west side of the central stadium (the maximum flood depth reaches more than 0.3 m), the intersection of Tingyun road and Heyuan West Road (the maximum flood depth is also above 0.3 m), and near Qinchunyuan supermarket (the maximum flood depth reaches more than 0.2 m). In addition, most of the maximum submergence depth in the study area is below 0.1 m. The simulated submergence position and submergence depth are consistent with the actual inundation situation in the study area. Therefore, the model can be used to simulate and evaluate the operation of the drainage system in the study area.

The operation of drainage system with return periods (2, 5, 10 and 20 years) was simulated by using the rainstorm intensity formula with stationary assumption as the input of the model. The inundation at RP 2 is similar to the observed rainstorm in 2017. With the increase of recurrence period, the location of the seriously affected area did not change significantly, but the depth and area of water accumulation increased significantly. In addition to several main pond points, the ponding mainly occurred on the road. In the case of RP 20, the main roads in the study area were flooded to varying degrees. Compared with RP 2, the inundated flood volume of RP 5, RP 10 and RP 20 increased by 80, 138 and 190%, respectively, and the maximum flood depth increased by 30, 46 and 86%, respectively. This is in accordance with the actual hydrological law and further illustrates the reliability of the model (Figure 8 and Table 4).

The nonstationary rainstorm intensity formula in different time periods is used to simulate the rainfall–runoff and compared with the stationary rainstorm intensity formula. There are obvious differences in flood range and depth in the whole study area (Figure 9 and Table 5). In terms of flood inundation area, although the serious inundation area did not change, the less serious area changed greatly and new areas were submerged. At RP 5, the difference ratio of submerged cells of (a), (b) and (c) relative to (d) is −24.5, 18.5 and 24.6%, respectively. At RP 2, this gap reaches the maximum. In terms of flood depth, the maximum flood depth of (d) is larger than that of (a), but it is obviously smaller than that of (b) and (c).

To compare the difference of submergence depth in the same place under different scenarios, 160 cells were randomly selected as the samples of each scenario, and the maximum submergence depth of 160 cells in three scenarios was extracted. Scatter plots were drawn for the maximum flood depths of 160 sampling points in each return period (Figure 10), where the X-, Y- and Z-axes represent the flood depths of the three cases (b), (c) and (d). In addition, the projection of sample points on three planes is also drawn. In all return periods, most of the points on the XZ- and YZ-planes are below the 1:1 straight line, which indicates that the submergence depth of (b) and (c) is greater than that of (d) at the same location. On the XY-plane, most of the points are located near the 1:1 straight line, but slightly higher than the straight line, which means that the submergence depth of (c) at the same location is slightly greater than but close to (b). In addition, with the increase of return period, the sample points gradually approach the 1:1 straight line, which indicates that the difference of submergence depth between different scenarios gradually decreases with the increase of return period.

To compare the overall difference of submergence depth in different scenarios, we calculated the flood depth of each unit and drew a violin diagram. From the violin diagram of flood depth in different return periods of four scenarios, it can be seen that the influence of rainstorm intensity formula in different scenarios on the unit flood depth is mainly reflected in the maximum flood depth and nuclear density, and the influence on the median is not obvious. In terms of the maximum flood depth, the maximum flood depth of (d) with a 2-year return period is greater than (a) and (b), close to (c), and (d) is the minimum in other periods. In terms of nuclear density, the violin diagram of (d) is ‘thinner’, which means more uniform distribution. In addition, under the assumption of nonstationarity, the violin diagram is more ‘fat’ over time, which means more concentrated distribution (Figure 11).

To analyze the response of detailed water accumulation process of urban surface to different scenarios, the water accumulation process lines of two main water accumulation points (intersection of Tingyun road and Heyuan West Road, and Qinchunyuan supermarket) were extracted. In the simulation process, these two water accumulation points are the two places most seriously affected. They are in the low-lying area and at the crossroads, with strong runoff generation capacity and low design standard of drainage network, so they are more vulnerable to flood damage. The flood hydrograph of (d) is much lower than that of (b) and (c) when the return period is 5 years. This means smaller flood depth and smaller flood volume. In addition, the time when the nonstationarity begins to submerge and the time when the flood depth reaches the peak are earlier than that of the stationary. This difference still exists once in 20 years, but it is not so significant. In the higher return period, the difference between different scenarios was not obvious (Figure 12).

The change of the nature of the Earth's climate and the acceleration of urbanization lead to the change of rainfall characteristics in urban areas, which destroys the stability of rainfall series. It is a popular method to overcome this difficulty to consider the nonstationary by adding covariates.

According to the results of GAMLSS model with time as covariate, the rainstorm intensity formula based on stationarity assumption does not consider the nonstationarity caused by environmental changes, which is quite different from the rainstorm intensity formula based on nonstationarity calculation. Although there is no specific rule for the parameters of short-duration rainstorm intensity formula, some clues can be found in the design rainstorm. The design rainstorm based on nonstationary design has higher total precipitation, higher peak precipitation intensity and more concentrated precipitation process. Such design rainstorms pose a greater threat to the urban drainage system, and such design rainstorms with greater threat will be more serious in the future.

The rainstorm intensity formulas of three different periods in nonstationary state are calculated and compared with those in stationary state. Although the difference of rainstorm intensity formula of four scenarios in urban drainage system simulation decreases with the increase of rainstorm scale, this difference is still very important even in high return period. When the return period is 20 years, the flood volume of (b) is 11.4% more than that of (d), and the flood volume of (c) is 18.5% more than that of (d). The flooded area of (b) is 8.8% more than that of (d), and the flooded area of (c) is 10.8% more than that of (d). The difference can even reach more than 50% in the case of small-scale rainstorm. The area that will not be submerged under the assumption of stationarity may be submerged in the case of nonstationarity. The nonstationary flood depth is higher than the stationary flood depth in the same place, and the flood depth in the future scenario is higher in the same place than in the current situation. Although the flood depth at the same location has increased, the distribution of flood depth is more concentrated and mainly concentrated below the 25% quantile. This shows that the flood depth of the newly added flooded cells is mainly lower, and no new flood hot spots can be found from the map of maximum flood depth. Among the two flood hotspots extracted, the flood in Qinchunyuan supermarket has subsided, while the flood at the intersection of Tingyun road and Heyuan West Road is more and more serious. On the whole, the flood hydrograph of nonstationary scenario is higher than that of stationary scenario. However, the difference is also different in the two places, especially in the flood hot spot of Qinchunyuan, which is mainly reflected in the process of recession. From various perspectives, nonstationary scenarios provide more dangerous urban drainage simulation results, and this risk is higher in low return period and future environment.

In this paper, we only use time t as the explanatory variable to analyze the nonstationarity of extreme rainstorm. In fact, there are still many physical factors that can be used as explanatory variables and have achieved good results (Gao et al. 2018; Rashid & Beecham 2019), such as North Pacific Oscillation (NPO), Pacific Decadal Oscillation (PDO), Arctic Oscillation (AO) and El-Niño Southern Oscillation (ENSO). In the future, more physical explanatory variables will be used to make up for the lack of using time as a covariate. In addition, we have established and verified the MIKE flood hydrodynamics model, but the model cannot completely replace the real ponding process. The error between the simulation results and the actual results cannot be avoided and uncertainty analysis is needed. Soil properties and land cover are very important for urban flood, while the underlying surface changes are not considered in this paper (Miller et al. 2014; Hossain Anni et al. 2020). The impact of extreme rainstorm and underlying surface changes on urban drainage systems may be comprehensively considered in the future.

The research results highlight the problems and challenges in the simulation of urban drainage systems under the changing environment: the nonstationarity of extreme rainfall series makes the rainstorm intensity formula based on stationary assumption show obvious deficiencies in the simulation of the urban drainage system, and this defect will be more important in the future. The characteristics of extreme rainstorm change rapidly and the hydrological design standards are insufficient. How to deal with the low design standards caused by the nonstationarity of extreme rainfall series under changing environments is also an important challenge for urban flood control and disaster prevention. The urbanization and sustainable development of China and other countries in the world need to study the current and future rainstorm intensity formula and its impact on urban drainage systems under changing conditions. Through GAMLSS and Mike flood models, the model discusses the influence of nonstationary short-duration rainstorm intensity formula on urban drainage systems, pointing out the risk of continuing to use the original scheme and providing a nonstationary solution considering the changing environment for the local people. These works are important for the coordination of municipal drainage engineering design, alleviating the increasingly serious urban flood disaster, and strengthening the construction of urban waterlogging prevention and drainage system makes sense.

The authors thank the Zhengzhou University logistics office for providing pipe network data and Zhengzhou Meteorological Bureau for rainfall data.

All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Z.W., S.L. and H.W. This revised version of the manuscript was written by S.L. and reviewed by Z.W. and H.W. All authors read and approved the final manuscript.

The study was funded by the Key Project of National Natural Science Foundation of China (No. 51739009); Natural Science Foundation of China (No. 51879242); Science and Technology Innovation Talents Project of Henan Education Department of China (No. 21HASTIT011); Young backbone Teachers Training Fund of Henan Education Department of China (No. 2020GGJS005) and Excellent Youth Fund of Henan Province of China (No. 212300410088).

The authors declare that there are no conflicts of interest.

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