Precipitation Trends Analysis Using Gridded Dynamic Sampling Zones: Case Study Yangtze Delta Megalopolis

Introduction

The impact of climate change and associated anthropogenic activities on the alteration of the urban hydrological regime and meteorological characteristics are an emerging and challenging subject of current urban hydrological research. The World Bank reported that the global air temperature in the middle of this century will rise by 2°C and reach 4°C by the end of the century, compared to the pre-industrial level in the 18th century (Potsdam Institute for Climate Impact Research and Climate Analytics, 2013). This will lead to an increase in the frequency and intensity of extreme rainfall, heat waves, and consequently, of droughts and floods. There is scientific consensus attributing global warming and more frequent climate extremes to human activities (IPCC et al., 2013a; IPCC et al., 2021b; Zhang et al., 2022). In the past 100 years (1908–2007), the climate change pattern observed in China was similar to the global one (Ding et al., 2006). Its average surface temperature increased by 1.1°C, slightly higher than the global average for the same period. The average annual rainfall showed an upward trend in the past 50 years, especially in the middle and lower reaches of the Yangtze River basin and the southeast of China. Both the frequency and intensity of climate extreme events have increased (Zhang et al., 2022; Ding et al., 2006; The State Council Information Office of the People’s Republic of China, 2008; National Development and Reform Commission (NDRC), 2013).

Since the 1980s, China’s economy has grown sharply, and the level of urbanization has risen drastically (Liu et al., 2021). The nominal gross domestic product (GDP) expanded 274 times from 367.87 billion CNY in 1978 to 100.88 trillion CNY in 2020, while the urban population percentage climbed from 17.92 to 63.89%. Although the level of urbanization is a critical indicator of social and economic development, the rapid urbanization in the last two decades has triggered a series of social-ecological-environmental processes, resulting amongst others into the occurrence of urban rain islands, urban heat islands, urban dust domes, and urban canyon effects (Central Committee of the Communist Party of China (CPC), 2014; Ding, 2018; Zhang et al., 2014). Compared with natural vegetation coverage, the urban heat island effect and dust cover effect aggravate variations in the regional spatial distribution of rainfall, resulting in a marked increase in rainfall peaks in built zones and downwind areas (Shepherd, 2005; Goddard Space Flight Center, 2012). The increase of impervious surfaces reduces evapotranspiration and infiltration, shortens the time of concentration, and increases the runoff coefficient (Zhang et al., 2014; Schueler, 1995; Federal Interagency Stream Restoration Working Group (FISRWG), 1989; Fletcher et al., 2013; Zhao et al., 2013; He et al., 2003). The changes of urbanization to the underlying surface have exacerbated the occurrence of extreme rainfalls, floods, and droughts (Research Council, 2008; Yu et al., 2014).

The Yangtze River Delta (YRD) has experienced rapid urbanization in the past decades and has become one of the most developed regions in China, forming the world-renowned Yangtze Delta Megalopolis (YDM) (Fang, 2019; Yu et al., 2021). While rapid development brought prosperity to the YDM, it also caused drawbacks challenging its livability, such as foggy haze (Meng et al., 2015; Liu et al., 2020), heat waves (Huang et al., 2008; Huang and Lu, 2015), and heavy precipitation (Fu et al., 2017; Ding, 2018). In the past decade, the intensity and frequency of extreme rainfall events and the resulting flood disasters exceeded the stationary reference standards for both urban drainage and flood control systems in the region (Wang et al., 2015; Xu et al., 2019; Xu et al., 2020). Extreme rainfall and snow destroy infrastructure services such as housing, transportation, communications, and electricity, threatening city safety. The impact of urbanization on regional and extreme rainfall has received extensive attention. Most of the studies on the impact of urbanization on urban rainfall are based on ground rain gauge stations (Table 1). The statistical rigor and hence the effectiveness of the methods applied are dependent on the availability of reliable and representative urban precipitation data. These are often limited as they require a dense gauge station network to achieve true urban coverage. In addition, these studies rarely consider urbanization as a dynamic phenomenon influencing the environment of the observation sites, such as a change of the underlying surface conditions due to urban densification and expansion. This hinders the exploration of the temporal variation of rainfall characteristics, which makes such an analysis inconsistent in time and space.

TABLE 1. Summary of the characteristics of the case study on precipitation.

It is common practice to examine the temporal and spatial distribution of urban precipitation using precipitation data from visible ground weather stations. Hu (2015) examined the characteristics of rainstorm-induced hazards in Beijing City and the effect of urbanization based on precipitation data derived from 20 precipitation gauge stations. Using precipitation data from 170 precipitation stations of Shanghai City, Chow and Chang (1984) plotted the isohyets of average annual precipitation, average wet season precipitation, and average dry season precipitation to analyze the relationship between precipitation distribution and city development. Significant urban rainfall islands were detected in the central region of Shanghai City based on the hourly precipitation of 11 precipitation stations by Ding (2018), who found that the maximum hourly precipitation occurred in the city center and gradually decreased from the city center to the suburbs. Fang et al. (2012) analyzed the daily precipitation at Baoshan Station in Shanghai and showed that the number of rainy days gradually decreased from 1971 to 2010, whilst the intensity and frequency of short-duration rainfall increased. The domain size, number of precipitation stations, and average density of research concentrating on differences in urban precipitation at different spatial scales are shown in Table 1. By dividing the domain area by the number of precipitation stations, the average density of precipitation stations was computed. The spatial scale of the study areas varied widely, ranging from cities to countries, natural river basins to the worldwide river system. Station network density ranged from 10.91 km² per station to 51,612.9 km² per station, with an average of 7,182.56 km² per station. Some scholars were overly optimistic in their assessment of site precipitation trends as regional trends, while others cautioned that site precipitation characteristics merely represent gauge station precipitation trends (Fang et al., 2012; Şen, 2014). The spatial distribution of precipitation stations employed in most research was excessively sparse when compared to the World Meteorological Organization’s (WMO) suggested minimum densities of precipitation stations (10–20 km² per station in urban areas) (World Meteorological Organization, 2008a; Chacon-Hurtado et al., 2017). When the density of precipitation stations is excessively high, the results of statistical analysis may be equivocal according to the large sample principle of statistical analysis, which states that the larger the sample size, the more they tend to their true value.

In a study spanning several decades, it is not uncommon that initially non-urban areas, or upwind/downwind zones, gradually become urban areas as a result of urban expansion. As a result, the weather stations in these shifting areas have transitioned from rural to urban stations. Many studies have used the existing underlying surface to differentiate between urban and non-urban areas when detecting precipitation trends without taking into account the passive changes in the site background environment generated by urbanization: for example, in Hu (2015), Hand and Shepherd (2009), Chen et al. (2016), and Kishtawal et al. (2010). Some researchers have noted that when cities grow, the classification of certain weather stations changes from non-urban to urban. To assess the impact of small and fragmented urban areas on precipitation along the Dutch West coast, Daniels et al. (2016) used land use data and geostrophic wind direction with a record length of 59 years for every 10 years to dynamically divide precipitation sites into rural and urban categories. By using remotely sensed image analysis, Zhang et al. (2021) discovered that land use and anthropogenic heat flow produced by urbanization had a considerable impact on the consistency of surface weather station data. Yu et al. (2022) identified urban and rural stations by using a method that combined a circular buffer and the percentage of urban area within the buffer to explore the shift pathways of the light and heavy precipitation associated with urbanization in the Yangtze River Delta region. Due to the availability of LULC data, there was no extra station identification for precipitation before 1980 and after 2015. In addition, the categories of stations were considered static in a few analyses, for example, the case analyzing the effect on sample size. Once the category of a rain gauge station being affected by urban expansion alters, it is often accompanied by a modification of the features of the domain that the rain station represents. This may lead to a loss of consistency in the rainfall data observed before and after the change (Dahmen and Hall, 1990; World Meteorological Organization, 2008b). If a trend analysis is performed directly on those precipitation stations while ignoring the change in their categories, the applicability of the precipitation characteristic data will be considerably reduced.

The purpose of this study is to describe and demonstrate a novel method based on a statistical analysis of remotely sensed data sets to overcome the sampling sufficiency challenge in rainfall trend analysis in terms of size and representativeness which acknowledges the dynamic nature of the surface area due to urbanization. The YDM region was selected as a case study as it experienced rapid expansion of its built area and frequent heavy precipitation events in recent decades.

Study Area and Data

The Yangtze River Delta is one of the most developed regions in China. The growth rate of its economic and urban development is ahead of the national average (Yu et al., 2021). The “Development Plan for the Yangtze Delta Megalopolis” was issued by the Chinese central government as one of the national strategies in 2016, in order to strengthen the coordination and cooperation in cross-city/regional development and enhance the international competitiveness of Chinese cities. The Yangtze Delta Megalopolis, located downstream of the Yangtze River Basin, has a land area of approximately 2.12 ×10⁵ km², and accommodates 26 cities, i.e. Shanghai, nine cities in Jiangsu Province (Nanjing, Wuxi, Changzhou, Suzhou, Nantong, Yancheng, Yangzhou, Zhenjiang, and Taizhou), eight cities in Zhejiang Province (Hangzhou, Ningbo, Jiaxing, Huzhou, Shaoxing, Jinhua, Zhoushan, and Taizhou), and eight cities in Anhui Province (Hefei, Wuhu, Maanshan, Tongling, Anqing, Chuzhou, Chizhou and Xuancheng) (Figure 1) (National Development and Reform Commission of the People’s Republic of China (NDRC), 2016). At the end of 2019, its registered population had reached 132.6 million, accounting for 9.47% of China’s population and 1.73% of the global population. Its nominal GDP in 2019 was 19.71 trillion CNY, accounting for 19.98% of China’s GDP and 3.55% of the global GDP, respectively.

FIGURE 1. The study area of the Yangtze Delta Megalopolis, China. Terrain tiles by Stamen Design (Terrain Classic Contributors, 2021) are used as a background.

Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) was jointly released by the University of California, Santa Barbara (UCSB) Climate Hazards Center and the United States Agency for International Development (USAID) in 2015 (Funk et al., 2015). This data set is based on the fusion of ground station and satellite-based remotely sensed data. Ranging from 1981 to near-present, it completely covers the YDM with a resolution of 0.05°×0.05° (approximately 5.6 km at the equator, less than 5 km in the YDM). Bai et al. (2018) assessed the accuracy of CHIRPS compared with rain gauges over mainland China and concluded that CHIRPS had a better performance in rainy areas than in arid or semi-arid areas. Tang et al. (2019) conducted a series of hydrological simulations using the Soil and Water Assessment Tool (SWAT) forced with different precipitation products in the Lancang River Basin. CHIRPS performed well in both the whole basin and sub-basins on the daily and monthly scales. Catholic University of Louvain (UCLouvain) has produced Global Land Cover (ESA-CCI-LC), a global land cover data product based on data obtained from the European Space Agency’s (ESA) Climate Change Initiative (CCI) (Defourny et al., 2017; Copernicus Climate Change Service, 2020; European Space Agency (ESA), 2015). With a horizontal resolution of 300 m and a time resolution of a year, this data collection follows the land cover categorization system (LCCS) defined by the United Nations (UN) Food and Agriculture Organization (FAO) (Gregorio and Jansen, 2000). On the premise of ensuring the availability of data in the YDM region, and after comparing the time series lengths and temporal and spatial resolution of the current available precipitation/land cover data, the CHIRPS data set as the daily precipitation data source was selected for this study. This allowed taking advantage of the high spatial resolution of 5 km as well as the long time coverage from 1981 to the present. ESA-CCI-LC was chosen as the basic land cover data in this study since it has a horizontal resolution of 300 m and a yearly time resolution from 1992 to the present.

Methodology

Generation of Dynamic Regions

The original boundary of the YDM was extended about 10 km outwards. Then, the extended YDM was divided into grids with the same size and location as the land use and precipitation data grids. The spatial coordinate reference system (CRS) used for grid creation was EPSG:32650; the affine transformation was set by defining the limits (west, south, east, north), length, and width, with parameters of 340,000, 3,065,000, 1,060,000, 3,875,000, 5,000, and 5,000, respectively. In the study area YDM, a total of 9,698 effective grids with a horizontal resolution of 5 km × 5 km were created (Figure 2). The gridding method provides the advantages of a large sample size, uniform distribution, and the same cell area, as shown in Figure 2.

FIGURE 2. Fishnets in the YDM.

To compare the impact of urbanization on regional precipitation, the ESA-CCI-LC data set was binarized/converted into land use maps with two classes of urban/built-up and non-urban/built-up based on the class value of the land cover being equal to 190 or not (190 for urban areas). All grids of the urban/built-up class were designated as urban region, the others were deemed non-urban region in the YDM region (Supplementary Figure S1). To ensure that the land use maps matched the precipitation data grid, they were resampled at a horizontal resolution of 5 km × 5 km. Supplementary Figure S1 depicts the 2020 binarized land use map.

The annual land use maps for urban and non-urban regions from 1992 to 2020 were generalized by the binarization approach to characterize the evolution dynamics of urbanization. After binarization, the grid number of urban and non-urban regions was counted year by year using the land use map. Figure 3 presents the trend of the grid number over time. Supplementary Figure S2 depicts the temporal and spatial distributions of urban construction in the YDM from 1995 to 2015.

FIGURE 3. Temporal variations of the sampling grid number for urban and non-urban regions in the YDM. The sum of the grid numbers of urban region and non-urban region always equals to 9,698.

Statistics of Precipitation

Using the yearly MDP as an example, the following statistical approach was developed for the above-mentioned year-by-year change of the gridded precipitation. Assuming that the dynamic area has n grids in the $x^{th}$ year, the following matrix can be used to express all of the precipitation sample values for that year:

where $p_{i,j}$ is the precipitation of the $j^{th}$ grid on the $i^{th}$ day of the $x^{th}$ year, in millimeters; m is the total number of days in a year, which is 365 for an average year and 366 for a leap year.

The first statistic was taken on the above matrix $P_x$ by extracting the maximum value of each column. These maxima represent the annual MDP of each grid in the $x^{th}$ year, which can be expressed as the vector:

where $p_{max,j}$ is the annual MDP of the $j^{th}$ grid the $x^{th}$ year, in millimeters. The extraction was conducted using the Google Earth Engine platform, the outputs of which were exported locally for subsequent use (Gorelick et al., 2017).

Then, in $P_{max,cell}^x$, the 0 values were removed as they represent a cell without precipitation, and the second statistic can be performed on the remaining. A variety of characteristic values can be chosen, such as the maximum value, the average value, the median, and so on, to meet specific demands. In this study the maximum, 75th percentile, 50th percentile (median), 25th percentile, and minimum as characteristic values, have been used which are recorded as $P_{max,max}^x$, $P_{max,75\text{\%}}^x$, $P_{max,50\text{\%}}^x$, $P_{max,25\text{\%}}^x$, and $P_{max,min}^x$, respectively.

Assuming $t$ – the total number of years, the time series of the characteristic values of the MDP can be constructed as follows:

Taking out the $max$:

The time series constructed above has been used to conduct a precipitation trend analysis by adopting a probabilistic analysis strategy with percentile statistics (Stensrud and Yussouf, 2007; Du and Chen, 2010; IPCC, 2021b). Meanwhile, a decadal average was performed to smooth out the fluctuations of the yearly statistics and highlight their trends (Akrami et al., 2014). To investigate the linear correlation for the time series between urban and non-urban regions, the Pearson correlation coefficient (PCC) and corresponding significance p-value were calculated by using the linregress function provided by the SciPy library (Taylor, 1990; Virtanen et al., 2020). The overall agreement was quantified using modified Kling-Gupta efficiency (KGE) (Gupta et al., 2009; Kling et al., 2012). The Mann-Kendall test supported by pyMannKendall was implemented to determine the trends of statistics (Yue and Wang, 2004; Hussain and Mahmud, 2019).

The statistical method comparable to annual MAP was made for finding total annual precipitation. Here, in Eq. 2 $max$ was substituted to calculate the sum precipitation on a grid,

where $p_{sum,j}$ is the total annual precipitation of the $j^{th}$ grid in the $x^{th}$ year, in millimeters.

Results and Discussion

Annual Maximum Daily Precipitation

As shown in Figure 3, the number of sampling grids exhibit a gradual increase in the urban region over time with a decrease in the non-urban region. It is to be expected that an expansion of the sampling regions will result in higher maxima and lower minima of daily precipitation than the unchanged sampling cells. To investigate how a change in sample size for the urban region affects the results, we compared the precipitation statistics for the expanding and unchanged urban regions using grids from 1992 to 2020. It follows from the results that these effects are too small (fall within the noise range) to be of particular concern for MDP in the urban regions before 2010 (Supplementary Figure S3). The interquartile ranges show a minor shift downwards from 2010 to 2020. These shifts are also quantitatively reflected in the MK values (Table 2). The effects of the reduced sampling cells are negligible in the non-urban regions (Supplementary Figure S4). Kling-Gupta efficiency coefficient (KGE) values reveal a similar pattern of all five quantiles of the annual MDP between the changing urban and non-urban regions and their unchanging parts, with values greater than 0.902 (Supplementary Table S1). In addition, the difference in sample size between urban and non-urban regions may also affect the results. We resampled an equal number of sampling cells from the non-urban region, equivalent to the number of sampling cells in the urban region, and extracted its five-number summary. After repeating 1,000 times, we found that the difference between the urban region and resampled non-urban regions was positively related to the difference between the urban and non-urban regions (Figures 4A–E). The extreme values of MDP amplified/decreased by the large sampling cells in the non-urban region. This amplification/decrease complicates an analysis of the relationship of precipitation extremes between the urban and non-urban regions, especially minima. As for the other three quartiles, they show great stability, irrespective of the sample size. Although the sample size matters, the difference in sample size disturbed the extrema of MDP slightly as well as strengthened the interquartile law under the context of a large sample strategy adopted in this study. Previous studies were hindered due to the limited number of ground stations and their representativeness (Sang et al., 2013; Han et al., 2015a; Han et al., 2015b; Wang et al., 2016; Cao et al., 2018; Pei et al., 2018; Wang et al., 2020; Yu et al., 2022). This limitation exacerbated after applying complex category identification of these stations (Diem and Mote, 2005; Yu et al., 2022). High-resolution reanalysis data such as CHIRPS provides a promising solution.

TABLE 2. The MK trend test for annual MDP in different percentiles of the YDM from 1992 to 2020.

FIGURE 4. The annual MDP in the YDM from 1992 to 2020. (A) maximum, (B) 75th percentile, (C) median, (D) 25th percentile, and (E) minimum represent the yearly changes. The boxplot of decadal averages is shown in (F).

The distribution of MDP has recognizable spatial characteristics within each year interval, but there is no unified feature between the years (Figure 5). However, Figures 6A,B depict that most of the maxima of MDP are located south of 31°N and most of the minima are located north of 31°N for both urban and non-urban regions, which may be related to topography and landform (Supplementary Figure S7). This is in line with observations of more heavy precipitation in the southern part of YRD reported in previous work (Sang et al., 2013; Han et al., 2015a; Wang et al., 2016).

FIGURE 5. The distribution of annual MDP in the YDM at quinquennial intervals from 1995 to 2020.

FIGURE 6. The distribution of yearly precipitation extrema in the YDM from 1995 to 2020 for (A) maxima of MDP, (B) minima of MDP, (C) maxima of TP, and (D) minima of TP.

Figures 4A–E show that the five quantiles of the yearly MDP increased to varying degrees in both the urban and non-urban regions. The percentiles of 25, 50, and 75% exhibit stable characteristics, ignoring the complexities of extrema that are influenced by the differences in sampling cells between urban and non-urban regions. All six time series in the MK trend test revealed a significant increasing trend, with five out of six (25% quantile and median of the annual MDP in the urban region and the 25%, median, and 75% in the non-urban region) meeting the significance level of 0.001 (Table 2). The distribution of difference for 10-year period averages of annual MDP in the YDM between the 1990s and 2010s is illustrated in Figure 7. In the past three decades, most cells have experienced an increase in annual MTP of decadal average. This increasing trend is uniform and significant for both urban and non-urban regions. If the non-urban region that is relatively less affected by human activities is regarded as a baseline, this explicitly indicates that the YDM climate change may be playing a dominant role in the trend variations of annual extreme precipitation. However, the yearly MDP in the YDM’s urban region is rising which may lead to a growing risk of urban flooding (Yu et al., 2022). This is in line with the observed increase in the frequency of flooding of many urban areas in recent years, as well as the higher economic losses caused by floods in the YDM. Figure 8 depicts a significant rise in economic losses due to urban flooding in the YDM region from 2006 to 2020 (Ministry of Water Resources of the People’s Republic of China (MWR), 2021). In addition, these results differ from previous studies on the impact of urbanization on hourly precipitation, which have shown a correlation between rapid urbanization and more extreme precipitation events (sub-daily) in urban areas (Liang and Ding, 2017; Ding, 2018; Jiang et al., 2020; Wang et al., 2021). In the comparison of these studies, this may be affected by data sources (station-/ non-station- based and the number of stations), statistic methods, and sizes and locations of case domains. And/or, it may be recommended that there are critical temporal resolutions greater than 3 h and less than 24 h in YDM, which separate these rules.

FIGURE 7. The distribution of precipitation difference for 10-year period averages in the YDM between the 1990s and 2010s for (A) MDP and (B) TP.

FIGURE 8. Historical flood loss in Jiang-Zhe-Hu-Wan region. Jiang, Zhe, Hu, and Wan are the abbreviations of Jiangsu Province, Zhejiang Province, Shanghai City, and Anhui province, respectively.

The Pearson correlation coefficient (PCC) and its p-value, as well as the Kling-Gupta efficiency coefficient (KGE) were used to assess the correlation between annual MDP of various quantiles across the urban and non-urban regions (Table 3). The line plots fluctuate similarly in the five quantiles of annual MDP in the urban and non-urban regions (Figures 4A–E), with values of PCC and KGE lying between 0.588–0.978 and 0.563–0.964, respectively (Table 3). The non-urban region has a significantly lower annual MDP than the urban region for the lowest, whereas the other quantiles in the non-urban region are generally higher than the urban region, with the maximum being the most significant. In both urban and non-urban regions, the 75th percentile, median, and 25th percentile of annual MDP are very close, and their PCC and KGE values are both above 0.9. According to the findings, there is a significant positive correlation between annual MDP in urban and non-urban areas, with the correlation in the middle three quantiles being significantly stronger than the correlation in the extreme quantiles. The yearly MDP of the urban region matches well with the non-urban region, with the middle three quantiles showing a better fit than the extreme quantiles. The yearly MDP fluctuation features in urban and non-urban locations are similar, with no significant differences, according to the strong correlation coefficients and satisfactory fitting methods.

TABLE 3. Correlation and fitting of the annual MDP in different percentiles between urban and non-urban regions of the YDM during 1992–2020.

The box plot’s interdecadal variations reveal a gradual increase in the 10-year period averages of annual MDP in both urban and non-urban areas from the 1990s to the 2010s (Figure 4F). Except for the lowest percentile, non-urban regions’ yearly MDP percentiles are greater than those in urban areas. Using the median as an example, in the 1990s, the decadal averages in urban and non-urban regions are nearly identical (76.99 and 76.61 mm), but in the 2010s, the averages grow to 91.53 and 95.53 mm, respectively, with percentage increases of 18.89 and 24.70%. Comparing the MK values of the 25th percentile, median, and 75th percentile in Table 2, the increase of the values in the non-urban region are greater than those in the urban region. Regardless of the impact of the sample size, this may be due to land cover and atmospheric humidity changes caused by urban expansion and human activities. Lu et al. (2019) found that the increase in the recurrence levels of annual MDP of non-urban stations is higher than that of urban stations considering the urbanized impacts. The non-urban areas with woodland, rainfed croplands, and grassland as the main underlying surface, which are the dominant land covers of precipitationsheds, may provide more favorable conditions for the generation of stronger heavy precipitation or heavier rainstorm centers at a daily scale (Keys et al., 2012). However, evapotranspiration is weakened in urban areas due to the loss of green vegetation. In comparison with the natural ground cover, the evapotranspiration decreases from 40 to 30% of the precipitation as the impervious surface area increases in urban areas (Federal Interagency Stream Restoration Working Group (FISRWG), 1989; Ellison et al., 2017). This process will result in a drop of atmospheric water vapor content, so the precipitation possibly decreases in the precipitationshed’s sink zone (Keys et al., 2012; Luo and Lau, 2019).

Annual Total Precipitation

An estimate of the impact of sample size changes on TP was conducted. Supplementary Figure S6 shows that the quantiles of TP for the decreasing non-urban regions are almost the same as those extracted from the unchanged region. The good fit reflected by these KGEs close to 1 are presented in Supplementary Table S1. Supplementary Figure S5 depicts that in the growing urban regions with more sampling cells the ranges of TP have expanded in certain years resulting from an increase of their extremes compared to the unchanged urban region. The interquartile ranges are a little bit narrowed due to the downward shift in the upper quartiles in the 2010s. Compared to the unchanged urban region, the MK values for percentiles of 25, 50, and 75% in the changing urban regions decrease as well (Table 4). Considering the large difference in sample size between non-urban and urban regions, Figures 9A–E reflect that this difference has little effect on the 25th percentile, median, and 75th percentile. There is a significant influence on most extrema, but the relative magnitude of TP extrema in urban and non-urban areas is rarely affected. Similar to annual MDP, impacts of sample size changes on annual TP are very limited, in contrast, the characteristics in the interquartile ranges are pretty stable under the large sample strategy.

TABLE 4. The MK trend test for annual TP in different percentiles of the YDM during 1992–2020.

FIGURE 9. The annual TP from 1992 to 2020 in the YDM. Yearly fluctuations of (A) maximum, (B) 75th percentile, (C) median, (D) 25th percentile, and (E) minimum. (F) boxplot of decadal averages.

Figure 10 depicts the spatial characteristics of the coordinated increase in annual TP from northern to southern of the YDM, which is valid for all years. They coincide with the contour map of the average annual precipitation released by the Ministry of Water Resources (Bureau of Hydrology, Ministry of Water Resources and Nanjing Hydraulic Research Institute, 2006). A similar spatial distribution of average annual precipitation was obtained using TRMM 3B43 (Cao et al., 2018). Compared with annual MDP, the extrema of annual TP are more obviously concentrated in a certain small area. The maxima gather in the southern part and most minima are crowded into the northern part. The annual summary consolidates the spatial features and eliminates the randomness of daily precipitation.

FIGURE 10. The distribution of annual TP in the YDM at quinquennial intervals from 1995 to 2020.

The five quantiles of the annual TP increase to varying degrees in both the urban and non-urban regions, as shown in Figures 9A–E. A significant rising trend is observed of the MK values of the lower, median, and upper quartiles of annual TP in both urban and non-urban regions. Except for the upper quartile in rural areas, the other five time series have a confidence level of 99.9% (Table 4). Of the 10 years averages, all sampling grids in YDM have witnessed the growth of TP from the 1990s to the 2010s (Figure 10). The observed rising trend in the annual TP in the YDM region reveals that regional water resources are more abundant in the 2010s than in the 1990s, as stated in China’s “National Assessment Report on Climate Change” (Lin et al., 2006). The significant increase in the amount of water resources in the YDM region will aid in resolving the contradictions between rapid economic development and rising water demands from industrial production, domestic water and ecological water requirements in urban areas, and rainfed agriculture water needs in non-urban areas.

The line charts fluctuate similarly in the five quantiles of annual TP in the urban and non-urban regions (Figures 9A–E). They are highly correlated and have a good fit with PCC and KGE of the middle three quartiles that are greater than or equal to 0.960 and 0.809, respectively (Table 5). Combined with the increasing trend, high correlation, and good agreement between urban and non-urban regions, climate change has impacted these features in the YDM.

TABLE 5. Correlation and fitting of the annual TP in different percentiles between urban and non-urban regions of the YDM during 1992–2020.

The box plot’s interdecadal changes reveal that yearly TP in urban and non-urban areas grows from the 1990s to the 2010s, after remaining relatively consistent during the 1990s and 2000s (Figure 9F). Except for the lowest percentile, the non-urban percentiles of the annual TP are greater than the urban percentiles. Using the median as an example, from 1990 to 2010, the decadal averages of annual TP in urban and non-urban regions grows from 1,214.39 mm and 1,251.09 mm to 1,379.18 mm and 1,456.61 mm, with increases of 13.57 and 16.43%, respectively. MK values of the 25th percentile, median, and 75th percentile of annual TP in the non-urban region are greater than those in the urban region. This is similar to the annual MDP with a slightly smaller difference between urban and non-urban regions. The non-urban region, with woods, rainfed croplands, and grassland as the primary land uses of precipitationsheds, is more likely to provide more favorable formation conditions for large-scale precipitation than the urban region.

Conclusion

This study provides a statistical method using grid-based precipitation characteristics to demonstrate a novel approach of dynamic sampling area to effectively address the limited representative space of ground precipitation sites in dynamic urban areas. The temporal and spatial inconsistency of static sampling, which is caused by urban growth, has effectively been addressed by this novel approach resulting in figures which are based on hundreds of thousands of sampling grids. In non-urban areas, the density of their fundamental data has surpassed the WMO suggested density, while in urban areas, it is very close to the required value (World Meteorological Organization, 2008a; Chacon-Hurtado et al., 2017). The statistics gathered by the newly created approach in this study are more in line with the requirements of large samples for statistical analysis as compared to traditional site-based research. What’s more, this method is flexible in the sense that the resolution of the reanalyzed product can be adjusted according to requirements.

Applying the foregoing methodologies for the analysis of precipitation features and trends in the YDM region revealed that:

(i) In the context of the large-sample strategy, the spatiotemporal changes of the sampling cells interfere with the limited extrema and have little effect on the laws within the interquartile ranges in the YDM.

(ii) The annual TP and MDP in the YDM region’s urban and non-urban areas both show a considerable upward trend, mainly due to climate change. They demonstrate that the region’s water supplies are expanding and that the risk of urban waterlogging is increasing.

(iii) Most maxima of annual MDP and TP are located in the south of YDM, while the minima are located in the north. Geographical differences affect the spatial distribution characteristics of precipitation in the YDM.

(iv) The increase in annual MDP and TP in the non-urban regions is greater than in the urban regions. It indicates that the land cover categories woods, grassland, and rainfed agricultural land may play a positive role in the formation of precipitationsheds’ sink zones.

Despite the rapid improvement of precipitation reanalysis products, their spatial and temporal resolution remains a limitation. Both horizontal and temporal resolution of existing products are insufficient for a more precise and accurate analysis of precipitation characteristics in urban areas.

To improve the current study, future work should concentrate on the following aspects:

(i) The impact of urbanization on its downwind areas. The contribution of the influence of urbanization at a scale of a megalopolis on precipitation on its downwind areas to the precipitation of non-urban areas needs further investigation.

(ii) Due to a lack of attribution analysis of precipitation change, it is difficult to differentiate between the increase of precipitation in the YDM region which is affected by anthropogenic activities and by global climate change. Therefore, an attribution analysis based on a statistical or deterministic study still needs to be considered in future work.

Data Availability Statement

The binarized land use images with a resolution of 5 km, as well as the extracted time series of percentile statistics in the YDM from 1992 to 2020, are available on Zenodo with a doi of https://doi.org/10.5281/zenodo.6640436 and are free to use and distribute under the CC BY 4.0 license.

Author Contributions

Conceptualization, XL; methodology, XL; software, XL; validation, XL; formal analysis, XL; investigation, XL; resources, XL; data curation, XL; writing—original draft preparation, XL and MY; writing—review and editing, CZ, AK, and MY; visualization, XL; supervision, DF and CZ; project administration, DF; funding acquisition, DF and MY. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Key R&D Program of China (2018YFC0809904) and National Natural Science Foundation of China (51909058).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Many thanks to Python community, Google Earth Engine, and QGIS. Land cover images from 1992 to 2020 in the YDM courtesy of the European Space Agency. The original boundary of the YDM courtesy of OpenStreetMap. Terrain tiles used as background in Figure 1 courtesy of Stamen Design. CHIRPS courtesy of UC Santa Barbara and USAID.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.917069/full#supplementary-material

References