Impact of Urbanization on Precipitation in North Haihe Basin, China

: The inﬂuence of urbanization on regional precipitation is one of the most important issues in hydrology. This paper selected the North Haihe Basin as the study area and explored the connection between summer precipitation and urbanization index (URBI) based on satellite precipitation and nighttime lights datasets. A moving spatial percentage anomaly method (MSP) was used to signify the local induced enhancement of precipitation (LIEP). Local indicators of spatial association (LISA) clustering for LIEP and URBI based on Bivariate Moran’s I coupled with digital elevation model (DEM) segmentation was used to separate the impacts caused by urbanization and terrain uplifting, and thus, the regions of interest (ROI) where the urban induced enhancement of precipitation (UIEP) plays a dominant role were located. A geographically weighted regression (GWR) model was used to reveal the spatial variation of the sensitivity of UIEP to URBI within the ROI. Pearson correlation and cross-wavelet analyses (XWT) were adopted to investigate the teleconnection between UIEP and climate anomalies, using Niño3.4 and SOI as indices. The results indicated that: (1) Urbanization e ﬀ ects on precipitation spatial variation in the upstream mountainous area would be hidden away by topographic factors. (2) From the perspective of the basin, the highest urbanization level areas have the Lowest LIEP, while the suburban areas have the highest LIEP, and the rural areas are in the middle. (3) The UIEP and the URBI are generally negatively correlated within ROI. (4) UIEP versus Niño3.4 and SOI both show a signiﬁcant common high-power period at a time scale of 2 years. This work can help comprehensively understand the hydrological response to urbanization.


Introduction
Precipitation is one of the most basic hydro-meteorological elements. Researchers have been aware of the urban influence on precipitation for many years [1][2][3]. However, many different viewpoints on the possible mechanisms related to urbanization-induced precipitation alterations still exist [3].
In the past few decades, human society has experienced a rapid urbanization process, and this trend is predicted to continue [4,5]. Thus, the influence of urbanization on regional precipitation has aroused many people's interests [6][7][8].
At present, studies of urban impacts on precipitation are commonly based on either observational data or numerical models. Observational data mainly involve rain-gauge data [9,10], remote sensing precipitation products [11], and radar quantitative precipitation estimation products [12]. Compared to the urban heat island (UHI) effect, urban effects on precipitation are not always concentrated in the city but change with various geographical conditions and climate backgrounds (a dry or wet climate situation). To evaluate urbanization impacts on zonal precipitation, a statistical analysis based on long-term observations of rainfall data is an acceptable way to distinguish likely influenced areas,

Study Area
This study takes the north part of the Haihe Basin as the study area ( Figure 1). The study area is located between 112°~120° E longitude and 38°~43° N latitude, facing the Bohai sea in the east, and the Inner Mongolian Plateau in the north. To the west is the Taihang Mountain and to the south is the North China plain. The area's climate covers semi-humid and semiarid continental monsoons. Under the influence of the East Asian summer monsoon (EASM) and the terrain conditions, the spatial and temporal distribution of precipitation shows clear zonal, seasonal, and inter-annual variation [27][28][29]. The average precipitation is 550 mm, most of which is concentrated in summer (June, July, and August, JJA). The study area is densely populated with large and medium-sized cities, among which Beijing, Tianjin and their surrounding cities are the core area of the Beijing-Tianjin-Hebei city cluster, one of the three major metropolitan areas of China.

Data Acuisition
In this paper, The Tropical Rainfall Measuring Mission (TRMM) multi-Satellite precipitation analysis products 3B43 V7 (Huffman et al., 2010) [30] from 1998 to 2018 were used. The resolution is 0.25° × 0.25° (https://pmm.nasa.gov) [31]. It has been validated that the TRMM 3B43 data are of high adaptability in the Haihe basin (Wei et al., 2017) [23]. The Digital Elevation Model (DEM) dataset was provided by the Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences (http://www.gscloud.cn) [32], with a 90 m × 90 m resolution. The Defense Meteorological Satellite Program/Operational Linescan System (DMSP/OLS, https://ngdc.noaa.gov/eog/dmsp/) [33] nighttime lights time series V4 acquired from NOAA's National Geophysical Data Center was applied to represent the urbanization level in the study area, which has a resolution of 1 km × 1 km [34,35]. The DMSP/OLS nighttime lights data have been widely adopted to quantize the urbanization level in studies related to urban issues [35][36][37][38]. These data can detect low-intensity lights in urban areas and distinguish them from dark rural backgrounds, thus, they are suitable to be used as an urbanization index (URBI). Considering that urban sprawl is relatively insignificant against the entire basin's background, and since the study period is not too long, the DMSP/OLS nighttime lights data of 2010 were used over the entire study period for simplicity. In this study, all the above data are resampled into 0.25° × 0.25°. The NCEP monthly mean wind reanalysis data provided by the NOAA Earth System Research Laboratory,

Data Acuisition
In this paper, The Tropical Rainfall Measuring Mission (TRMM) multi-Satellite precipitation analysis products 3B43 V7 [30] from 1998 to 2018 were used. The resolution is 0.25 • × 0.25 • [31]. It has been validated that the TRMM 3B43 data are of high adaptability in the Haihe basin [23]. The Digital Elevation Model (DEM) dataset was provided by the Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences [32], with a 90 m × 90 m resolution. The Defense Meteorological Satellite Program/Operational Linescan System [33] nighttime lights time series V4 acquired from NOAA's National Geophysical Data Center was applied to represent the urbanization level in the study area, which has a resolution of 1 km × 1 km [34,35]. The DMSP/OLS nighttime lights data have been widely adopted to quantize the urbanization level in studies related to urban issues [35][36][37][38]. These data can detect low-intensity lights in urban areas and distinguish them from dark rural backgrounds, thus, they are suitable to be used as an urbanization index (URBI). Considering that urban sprawl is relatively insignificant against the entire basin's background, and since the study period is not too long, the DMSP/OLS nighttime lights data of 2010 were used over the entire study period for simplicity. In this study, all the above data are resampled into 0.25 • × 0.25 • . The NCEP monthly mean wind reanalysis data provided by the NOAA Earth System Research Laboratory, Physical Sciences Division (NOAA/OAR/ESRL PSD), Boulder, Colorado, USA, from their Web site Atmosphere 2020, 11, 16 4 of 16 at https://www.esrl.noaa.gov/ [39], were used to obtain the summer mean wind direction, with a resolution of 2.5 • × 2.5 • .

Empirical Orthogonal Function (EOF) Analysis
The empirical orthogonal function (EOF) and the principal components analysis (PCA) were sometimes used to describe the same method [40], but EOF can emphasize both temporal coefficient and spatial pattern [41]. The EOF was applied to realize the decomposition of the summer precipitation, and both spatial patterns and temporal coefficients were obtained. The spatial patterns represent the main spatial feature of summer precipitation, and the temporal coefficients are the principal components.
The EOF is based on the estimation of the eigenvalues and eigenvectors of an anomaly covariance matrix of the precipitation field. The derived eigenvalues provide a measure of the percent variance explained by each spatial pattern. North's test [42] is selected to determine if a particular eigenvalue (pattern) is distinct from its nearest neighbor.
The purpose of the EOF analysis is to compare the obtained main spatial pattern of summer precipitation with urbanization spatial distribution delineated by nighttime lights data. The result can offer a preliminary assessment of urbanization influence on the spatial rainfall patterns at the basin scale. The EOF analysis was also used in the decomposition of the urban induced enhancement of the precipitation index (UIEP), and temporal coefficients of the first EOF pattern were used as a substitution for the UIEP field in the correlation analysis with climate anomalies.

Local Induced Enhancement of Precipitation
This article uses a moving spatial percentage anomaly (MSP) to signify the local high/low points; this is a filtering process [43,44]. Assuming that the climate change beyond the moving-window scale is caused by large-scale natural climate forces, then the uneven spatial distribution after filtering can be attributed to the effects of local climate forces (such as urbanization, land use change, abrupt terrain change, and anthropogenic aerosol emissions). The value of MSP was taken as the LIEP index. The steps for calculating LIEP using MSP are as follows. First, a rectangular moving window with certain size (larger than 0.25 • × 0.25 • ) centered on each cell of the TRMM 3B43 data was introduced. Second, for every year, the mean summer rainfall P mean within the whole window and the mean rainfall P c on the central cell were calculated. Then, the percentage difference between the central cell value P c and the P mean can be easily obtained as the LIEP: If the LIEP value is positive, then the local elements enhance the precipitation, where the corresponding region tends to have a "wet island effect". On the contrary, if the LIEP value is negative, then the local elements inhibited the precipitation, where the corresponding region tends to show a "dry island effect". If the LIEP value is zero, the local elements' effects are not obvious.
In the method mentioned above, the moving window's size would definitely impact the results of the rainfall's spatial heterogeneity. In order to reasonably determine the scope of the moving window, a spatial auto-correlation analysis of the precipitation in different years was carried out. The spatial autocorrelation of the precipitation would remain positive and strong within a certain distance, which suggested that the precipitation within a limited range tended to be affected by the same local climate force. The spacing within which the autocorrelation of precipitation remains positive was thus determined to be in the range of the moving window.  [45,46] for the LIEP versus Urbanization index was adopted to extract the region where the urban effect may play a dominant role in influencing precipitation.
The bivariate Moran's I is extended from the traditional method of Moran's I and can identify the clustering of extreme values between two variants. The bivariate Moran's I is defined as in Equation (2): where x i k indicates the value of feature k on space unit i; x j l indicates the value of feature l on space unit j; x k and x l are the mean value of feature k and l, respectively; w ij is the spatial weight between these two features; and S 2 is the covariance of features k and l.
The bivariate Local Moran's I has been used to produce local indicators of spatial association (LISA) measure. The calculation of bivariate Local Moran's I is shown in Equation (3): where S 2 k and S 2 l are the observed variance of feature k and l, respectively. It is worth noting that the bivariate Moran's I does not suggest an in-place inherent correlation between x i k and x i l , but instead between x i k and Σ j w ij x j l . Through the LISA cluster map, prepared with Geoda [47,48] (http://geodacenter.github.io), locations with significant local Moran's statistics can be obtained. We are interested in the cluster types with a High LIEP value versus a high urban level (H-H) and a Low LIEP value versus a high urban level (L-H). Coupled with the DEM data, areas with large topographic uplift need to be removed to exclude the interference of topography alteration elements. Then, the region of interest (ROI) can be obtained, where the LIEPs are considered to be dominantly caused by the urbanization factor (land use change, anthropogenic aerosol emissions, etc.) and are thus expressed as the Urban induced enhancement of precipitation (UIEP).

Geographically Weighted Regression, GWR
The geographic weighted regression model (GWR) is an extension of ordinary least squares regression (OLS) and adds a level of modeling sophistication by allowing the relation between the independent and dependent variables to vary by locality. GWR constructs a separate OLS equation for every location in the dataset, which incorporates the dependent and explanatory variables of the locations falling within the bandwidth of each target location [49][50][51]. The bandwidth is the key controlling parameter in the weighting scheme, which can be optimized by the Akaike Information Criterion (AIC) [52]. This paper uses the GWR model to interpret the spatial variation of the sensitivity of UIEP to URBI.

Pearson Correlation and Cross-Wavelet Analysis
The Pearson correlation is used to investigate a simple linear relation. In order to further reveal the specific details of the non-linear relation from both the time domain and frequency domain, a cross-wavelet analysis (XWT) was adopted. XWT is an effective tool for exploring the correlations between two time series, as it combines the characteristics of wavelet transform and cross spectrum analysis. If the wavelet transforms of the x t and y t series are W x n (s, t) and W y n (s, t), respectively, the cross-wavelet spectrum can be expressed as Equation (4) [53]: where s and t, respectively, denote the scale expansion factor and the time shift factor, and the * refers to complex conjugation. The modulus of the cross-wavelet transform W xy n (s, t) indicates the cross-amplitude, where a larger value indicates a higher degree of correlation.
The phase difference between these two series for a given delay time and a specific scale can be defined as Equation (5): where S is the smooth operator and R and I represent the real and imaginary parts of the cross-wavelet spectrum.

The Summer Rainfall Spatial Distribution Characteristics in North Haihe Basin
For the purpose of understanding the spatial and temporal distribution characteristics of the summer rainfall, an EOF analysis was performed on the 3B43 precipitation data. The first three spatial patterns and their corresponding temporal coefficients are shown in Figure 2. Their cumulative variance contribution rate is 0.7 and passes North's test [42], which indicates that these eigenvalues are significantly separated. The phase difference between these two series for a given delay time and a specific scale can be defined as Equation (5): where S is the smooth operator and R and I represent the real and imaginary parts of the cross-wavelet spectrum.

The Summer Rainfall Spatial Distribution Characteristics in North Haihe Basin
For the purpose of understanding the spatial and temporal distribution characteristics of the summer rainfall, an EOF analysis was performed on the 3B43 precipitation data. The first three spatial patterns and their corresponding temporal coefficients are shown in Figure 2. Their cumulative variance contribution rate is 0.7 and passes North's test (North et al., 1982) [42], which indicates that these eigenvalues are significantly separated. As shown in Figure 2, the spatial pattern of the first leading spatial pattern (EOF1) mainly reflects the gradual precipitation amount changes from the coastal to inland and also from plain to mountainous areas. These gradual changes eventually lead to the distinct difference in precipitation As shown in Figure 2, the spatial pattern of the first leading spatial pattern (EOF1) mainly reflects the gradual precipitation amount changes from the coastal to inland and also from plain to mountainous areas. These gradual changes eventually lead to the distinct difference in precipitation distribution between the mountainous upstream area and the flat middle and downstream area. The correlation coefficient between the temporal coefficient of the first pattern (PC1) and the average summer precipitation amount of the basin reaches −0.57 (at a significance level of 0.01). This may be related to the control of large-scale weather systems in the summer. Since the PC1 is entirely negative, the negative loading region in EOF1 from middle and downstream areas to Bohai gulf contributes to large portions of the basin's precipitation.
The second EOF pattern (EOF2) shows an obvious south-to-north reverse change. Combined with the temporal coefficient (PC2), this pattern reveals the fact that there are more chances for the south region to have greater precipitation. Both the distributions of EOF1 and EOF2 are related to the large-scale climate background and geographical environment.
However, in the third EOF pattern (EOF3), heavier negative loading is located at the north of Tianjin and the northeast of Beijing. This location is coincident with the windward slope of the mountains in the north Haihe basin, and this area is also the downwind area of Beijing-Tianjin (shown in Figure 3; the downwind area is located in the northeast of these two cities). This might be the consequence of the bifactor functioning of topographic forces and urban effects. Within the plain area at the southeast of the study area, the loading in both the Beijing-Tianjin urban and peri-urban areas shows a slight spatial difference with the exterior zone. Compared to prior patterns, the EOF3 reflects the spatial difference of precipitation induced by local climate forces. Taken together, the contribution rate of urban effects alone to the spatial differentiation of rainfall is relatively small at the basin scale, and the effects are always mixed up with local terrain influence. Thus, this non-urban factors' influence should be removed based on regional spatial and temporal variation characteristics when assessing urban impacts on precipitation.
Atmosphere 2019, 10, x FOR PEER REVIEW 7 of 16 distribution between the mountainous upstream area and the flat middle and downstream area. The correlation coefficient between the temporal coefficient of the first pattern (PC1) and the average summer precipitation amount of the basin reaches −0.57 (at a significance level of 0.01). This may be related to the control of large-scale weather systems in the summer. Since the PC1 is entirely negative, the negative loading region in EOF1 from middle and downstream areas to Bohai gulf contributes to large portions of the basin's precipitation. The second EOF pattern (EOF2) shows an obvious south-to-north reverse change. Combined with the temporal coefficient (PC2), this pattern reveals the fact that there are more chances for the south region to have greater precipitation. Both the distributions of EOF1 and EOF2 are related to the large-scale climate background and geographical environment.
However, in the third EOF pattern (EOF3), heavier negative loading is located at the north of Tianjin and the northeast of Beijing. This location is coincident with the windward slope of the mountains in the north Haihe basin, and this area is also the downwind area of Beijing-Tianjin (shown in Figure 3; the downwind area is located in the northeast of these two cities). This might be the consequence of the bifactor functioning of topographic forces and urban effects. Within the plain area at the southeast of the study area, the loading in both the Beijing-Tianjin urban and peri-urban areas shows a slight spatial difference with the exterior zone. Compared to prior patterns, the EOF3 reflects the spatial difference of precipitation induced by local climate forces. Taken together, the contribution rate of urban effects alone to the spatial differentiation of rainfall is relatively small at the basin scale, and the effects are always mixed up with local terrain influence. Thus, this non-urban factors' influence should be removed based on regional spatial and temporal variation characteristics when assessing urban impacts on precipitation.

Local Induced Enhancement of Precipitation and Urban Distribution
To reasonably determine the size of the moving window in MSP, the spatial autocorrelation of rainfall in different level years were estimated. As described in the previous section, the precipitation in the study area showed an obvious trend (from east to west and from south to north). The summer precipitation anomaly percentage was calculated in every cell through the TRMM 3b43 monthly data. In this way, the trends can be eliminated. As shown in Figure 4, the precipitation anomaly percentage shows a positive spatial correlation until the distance is larger than 1.5 to 2.5 degrees. This suggests that when the distance between two rainfall cells is shorter, the yearly precipitation change tends to be more clearly synchronous. Thus, rainfall within 1.5 to 2.5 degrees could be considered impacted by the same local climate forcing factors. This paper set the size of the moving window in MSP as 1.75 degrees.

Local Induced Enhancement of Precipitation and Urban Distribution
To reasonably determine the size of the moving window in MSP, the spatial autocorrelation of rainfall in different level years were estimated. As described in the previous section, the precipitation in the study area showed an obvious trend (from east to west and from south to north). In this way, the trends can be eliminated. As shown in Figure 4, the precipitation anomaly percentage shows a positive spatial correlation until the distance is larger than 1.5 to 2.5 degrees. This suggests that when the distance between two rainfall cells is shorter, the yearly precipitation change tends to  The MSP was performed, and the LIEP distribution is shown in Figure 5a. As shown in Figure  5a, for LIEP, the most obvious high value area is located downwind of Beijing-Tianjin-Tangshan. This location happens to be on the windward slope of the northern mountains, which overlaps considerably with the high loading area of the EOF2 mode of summer rainfall (Figure 2c). Meanwhile, Zhang et al. (2015) [20] validated that the change of land use in Beijing-Tianjin-Tangshan practically led to an increase in rainfall volume and frequency in this area based on simulation by WRF coupled with the UCM model. Thus, the LIEP of this area should be influenced by both the terrain uplift and urbanization effects. According to the basin topography variation and urban distribution patterns, the areal LIEP distribution above the 500 m contour shows that the boundary of high and low LIEP is mostly concentrated in the zigzag part of the topographic contour but without an obvious relation to urban distribution. That means that precipitation here is mostly influenced by complex local terrain deformation. Even if urban impacts exist, they would be covered up by the counterparts precipitated by the terrain factor. In the plain area, however, where the elevation is below the 50 m contour, the terrain is flat. It can be seen that several important urban centers (such as Beijing, Tianjin, Tangshan, and Baoding city) and their surrounding areas show obviously low LIEP concentrations, which can be considered as a response to urbanization. For the purpose of comparison with LIEP, the areal urbanization level distribution is provided in Figure 5b, where L1 to L5 represent the urbanization level from low to high. In this article, the The MSP was performed, and the LIEP distribution is shown in Figure 5a. As shown in Figure 5a, for LIEP, the most obvious high value area is located downwind of Beijing-Tianjin-Tangshan. This location happens to be on the windward slope of the northern mountains, which overlaps considerably with the high loading area of the EOF2 mode of summer rainfall (Figure 2c). Meanwhile, Ref. [20] validated that the change of land use in Beijing-Tianjin-Tangshan practically led to an increase in rainfall volume and frequency in this area based on simulation by WRF coupled with the UCM model. Thus, the LIEP of this area should be influenced by both the terrain uplift and urbanization effects. According to the basin topography variation and urban distribution patterns, the areal LIEP distribution above the 500 m contour shows that the boundary of high and low LIEP is mostly concentrated in the zigzag part of the topographic contour but without an obvious relation to urban distribution. That means that precipitation here is mostly influenced by complex local terrain deformation. Even if urban impacts exist, they would be covered up by the counterparts precipitated by the terrain factor. In the plain area, however, where the elevation is below the 50 m contour, the terrain is flat. It can be seen that several important urban centers (such as Beijing, Tianjin, Tangshan, and Baoding city) and their surrounding areas show obviously low LIEP concentrations, which can be considered as a response to urbanization.  The MSP was performed, and the LIEP distribution is shown in Figure 5a. As shown in Figure  5a, for LIEP, the most obvious high value area is located downwind of Beijing-Tianjin-Tangshan. This location happens to be on the windward slope of the northern mountains, which overlaps considerably with the high loading area of the EOF2 mode of summer rainfall (Figure 2c). Meanwhile, Zhang et al. (2015) [20] validated that the change of land use in Beijing-Tianjin-Tangshan practically led to an increase in rainfall volume and frequency in this area based on simulation by WRF coupled with the UCM model. Thus, the LIEP of this area should be influenced by both the terrain uplift and urbanization effects. According to the basin topography variation and urban distribution patterns, the areal LIEP distribution above the 500 m contour shows that the boundary of high and low LIEP is mostly concentrated in the zigzag part of the topographic contour but without an obvious relation to urban distribution. That means that precipitation here is mostly influenced by complex local terrain deformation. Even if urban impacts exist, they would be covered up by the counterparts precipitated by the terrain factor. In the plain area, however, where the elevation is below the 50 m contour, the terrain is flat. It can be seen that several important urban centers (such as Beijing, Tianjin, Tangshan, and Baoding city) and their surrounding areas show obviously low LIEP concentrations, which can be considered as a response to urbanization. For the purpose of comparison with LIEP, the areal urbanization level distribution is provided in Figure 5b, where L1 to L5 represent the urbanization level from low to high. In this article, the  For the purpose of comparison with LIEP, the areal urbanization level distribution is provided in Figure 5b, where L1 to L5 represent the urbanization level from low to high. In this article, the urbanization level was derived from DMSP/OLS nightlight data, which were classified into five levels using the Jenks Natural Breaks method. Larger cities (areas with L4 or above and an extent larger than 2 cells) are mostly concentrated in the region below the 500m contour (the thick purple contour). Therefore, the spatial variation of precipitation caused by urbanization at the basin scale would mainly occur in the middle and downstream areas.
Statistical analysis was carried out on the LIEP in different urbanization level areas, and the boxplot result is shown in Figure 6. These results suggest that in the area where the urbanization level is high (L4 and L5), the summer LIEPs are relatively lower; even the upper quartile is less than zero. This means that the multi-year mean rainfall in a metropolis tends to be less than that in a peri-urban area. In the modest urbanization level (L3) areas located in the urban and rural junction zones, the LIEPs seem to be a bit larger, with a median greater than zero. However, the LIEPs in these areas show a larger deviation. In the rural area (L2, L1), the LIEP distribution is more even compared to the above areas.
Atmosphere 2019, 10, x FOR PEER REVIEW 9 of 16 urbanization level was derived from DMSP/OLS nightlight data, which were classified into five levels using the Jenks Natural Breaks method. Larger cities (areas with L4 or above and an extent larger than 2 cells) are mostly concentrated in the region below the 500m contour (the thick purple contour). Therefore, the spatial variation of precipitation caused by urbanization at the basin scale would mainly occur in the middle and downstream areas.
Statistical analysis was carried out on the LIEP in different urbanization level areas, and the boxplot result is shown in Figure 6. These results suggest that in the area where the urbanization level is high (L4 and L5), the summer LIEPs are relatively lower; even the upper quartile is less than zero. This means that the multi-year mean rainfall in a metropolis tends to be less than that in a peri-urban area. In the modest urbanization level (L3) areas located in the urban and rural junction zones, the LIEPs seem to be a bit larger, with a median greater than zero. However, the LIEPs in these areas show a larger deviation. In the rural area (L2, L1), the LIEP distribution is more even compared to the above areas.

Identifying the Precipitation Response to Urbanization
By using bivariate Moran's I, LIEP was taken as the first variable, and URBI-represented by nighttime light data-was taken as the second variable for the spatial correlation analysis. The corresponding bivariate Moran's I was estimated to be −0.11 (p < 0.001). This indicates that LIEP has a negative correlation with the surrounding area's URBI. The LISA cluster map for LIEP and URBI is shown in Figure 7a.

Identifying the Precipitation Response to Urbanization
By using bivariate Moran's I, LIEP was taken as the first variable, and URBI-represented by nighttime light data-was taken as the second variable for the spatial correlation analysis. The corresponding bivariate Moran's I was estimated to be −0.11 (p < 0.001). This indicates that LIEP has a negative correlation with the surrounding area's URBI. The LISA cluster map for LIEP and URBI is shown in Figure 7a.
The upstream region above the 500 m contour mainly features mountains and hills, with a low urban density and low urbanization level. Significant cluster types (significance level < 0.05, shown in Figure 7b) are defined as the high LIEP versus low URBI (H-L) and low LIEP versus low URBI (L-L), which suggests that the variance in LIEP is mainly influenced by the terrain factor, and the urban effects are obfuscated.
The middle and downstream region below the 500 m contour (the thick one) includes high urban density and a high urbanization level. The topography is a flat plain, thus, the terrain deformation effects can be excluded in the analysis of precipitation response. Significant cluster types (significance level < 0.05) consist of high LIEP versus High URBI (H-H) and Low LIEP versus High URBI (L-H). The inner parts of the city agglomeration show an L-H cluster, which suggests that highly urbanized areas induced a decrease of the summer precipitation compared to the surrounding areas. This result is consistent with the conclusions of [20] and [54]. The edge of the urban agglomeration shows significant H-H clustering, indicating that the suburban area can induce a certain increase in precipitation, which further confirms the results in Figure 6. The upstream region above the 500 m contour mainly features mountains and hills, with a low urban density and low urbanization level. Significant cluster types (significance level < 0.05, shown in Figure 7b) are defined as the high LIEP versus low URBI (H-L) and low LIEP versus low URBI (L-L), which suggests that the variance in LIEP is mainly influenced by the terrain factor, and the urban effects are obfuscated.
The middle and downstream region below the 500 m contour (the thick one) includes high urban density and a high urbanization level. The topography is a flat plain, thus, the terrain deformation effects can be excluded in the analysis of precipitation response. Significant cluster types (significance level < 0.05) consist of high LIEP versus High URBI (H-H) and Low LIEP versus High URBI (L-H). The inner parts of the city agglomeration show an L-H cluster, which suggests that highly urbanized areas induced a decrease of the summer precipitation compared to the surrounding areas. This result is consistent with the conclusions of Zhang et al. (2015) [20] and Zhao et al. (2018) [54]. The edge of the urban agglomeration shows significant H-H clustering, indicating that the suburban area can induce a certain increase in precipitation, which further confirms the results in Figure 6.
From the above, the LIEP of the region below the 500 m contour clumped in the L-H and H-H clusters embodies the influence caused by urban effects, with the terrain effects removed. Thus, the LIEP in these parts can be restated as the UIEP. In general, a bivariate Moran's I for LIEP versus URBI coupled with DEM segmentation could effectively extract the urban agglomeration's area of influence on precipitation form the north Haihe basin. Thus, these areas were taken as the ROI for the next analysis on the correlation between UIEP and potential influencing factors.

The Spatial Correlation between UIEP and URBI
It can be seen that the inner parts of ROI show a "dry island effect", while the edge parts show a "wet island effect" (Figures 6 and 7). Since the precipitation response to the urbanization level in different parts of the urban agglomeration might vary, the spatial non-stationarity of their relation needs to be investigated. For this purpose, a GWR model was established to analyze the spatial correlation between the multi-year averages of UIEP and URBI. The R 2 and Akaike information criteria (AIC) are, respectively, 0.82 and −283.71. The estimated regression coefficients are shown as Figure 8. In the summer, "dry island effect" dominated inner parts, the estimated coefficients of URBI are negative but infinitesimal. This means that higher URBI may lead to a more remarkable "dry" condition but with a slight sensitivity. In the "wet island effects" dominated edge parts, where spots are marked with a red border, the estimated coefficients of URBI are still negative but larger. From the above, the LIEP of the region below the 500 m contour clumped in the L-H and H-H clusters embodies the influence caused by urban effects, with the terrain effects removed. Thus, the LIEP in these parts can be restated as the UIEP. In general, a bivariate Moran's I for LIEP versus URBI coupled with DEM segmentation could effectively extract the urban agglomeration's area of influence on precipitation form the north Haihe basin. Thus, these areas were taken as the ROI for the next analysis on the correlation between UIEP and potential influencing factors.

The Spatial Correlation between UIEP and URBI
It can be seen that the inner parts of ROI show a "dry island effect", while the edge parts show a "wet island effect" (Figures 6 and 7). Since the precipitation response to the urbanization level in different parts of the urban agglomeration might vary, the spatial non-stationarity of their relation needs to be investigated. For this purpose, a GWR model was established to analyze the spatial correlation between the multi-year averages of UIEP and URBI. The R 2 and Akaike information criteria (AIC) are, respectively, 0.82 and −283.71. The estimated regression coefficients are shown as Figure 8. In the summer, "dry island effect" dominated inner parts, the estimated coefficients of URBI are negative but infinitesimal. This means that higher URBI may lead to a more remarkable "dry" condition but with a slight sensitivity. In the "wet island effects" dominated edge parts, where spots are marked with a red border, the estimated coefficients of URBI are still negative but larger. This indicates that the "wet island effect" is inhibited by the urbanization to some extent and is more sensitive to the URBI than the inner areas.
Atmosphere 2019, 10, x FOR PEER REVIEW 11 of 16 This indicates that the "wet island effect" is inhibited by the urbanization to some extent and is more sensitive to the URBI than the inner areas. Squares outlined in red represent the "wet island effects" dominated area.

Lagged Correlation of UIEP and Climate Anomalies
As mentioned above, URBI has a significant impact on the spatial distribution of the multi-year average UIEP, but the varying strength of UIEP in each year is considered to be influenced by the climatic background (Zhang et al., 2018) [25]. The ROI area is located in the east Asian monsoon

Lagged Correlation of UIEP and Climate Anomalies
As mentioned above, URBI has a significant impact on the spatial distribution of the multi-year average UIEP, but the varying strength of UIEP in each year is considered to be influenced by the climatic background [25]. The ROI area is located in the east Asian monsoon (EASM) area, and the summer precipitation is controlled by the EASM [55]. The UIEP is thus expected to be influenced by the EASM. Since the ENSO was considered to have a great impact on EASM [56], the potential association between UIEP and ENSO is worthy of studying. We calculated the 1-12-month lag correlation fields of UIEP versus Niño3.4 and SOI from 1998 to 2018. Given the space limitations, only 0-5-month lag correlation fields are shown in Figures 9 and 10.   Figure 7. The distribution of the estimated regression coefficient of urbanization index (UBRI). Squares outlined in red represent the "wet island effects" dominated area.

Lagged Correlation of UIEP and Climate Anomalies
As mentioned above, URBI has a significant impact on the spatial distribution of the multi-year average UIEP, but the varying strength of UIEP in each year is considered to be influenced by the climatic background   [25]. The ROI area is located in the east Asian monsoon (EASM) area, and the summer precipitation is controlled by the EASM (Huang et al., 2005) [55]. The UIEP is thus expected to be influenced by the EASM. Since the ENSO was considered to have a great impact on EASM (Hao and Ding 2012) [56], the potential association between UIEP and ENSO is worthy of studying. We calculated the 1-12-month lag correlation fields of UIEP versus Niño3.4 and SOI from 1998 to 2018. Given the space limitations, only 0-5-month lag correlation fields are shown in Figures 9 and 10.  From Figures 9 and 10, it can be seen that the regional UIEP shows spatial differentiation for a correlation with climate anomalies during 1998-2018. In Figure 9, the sites where the correlation between UIEP and Niño3.4 is strong appear frequently around the urban core area (marked as a square with a dot in the center), while in the urban core, this association is relatively weak. Among them, the junction zone of Beijing and Tianjin and northwest of Beijing show the most intensive correlation. As mentioned above, urbanization generally induced a "dry island" and a "wet island", From Figures 9 and 10, it can be seen that the regional UIEP shows spatial differentiation for a correlation with climate anomalies during 1998-2018. In Figure 9, the sites where the correlation between UIEP and Niño3.4 is strong appear frequently around the urban core area (marked as a square with a dot in the center), while in the urban core, this association is relatively weak. Among them, the junction zone of Beijing and Tianjin and northwest of Beijing show the most intensive correlation. As mentioned above, urbanization generally induced a "dry island" and a "wet island", respectively, in the inner area and the edge part of the urban agglomeration from a multi-year average precipitation perspective. However, the pattern in which the UIEP is impacted by climate anomalies tends to be different. The highest urbanized core area, which is the low value center of UIEP, tends to be weakly related to Niño3.4, while the peri-urban-core area may have a stronger relation with Niño3.4, but the direction of this association is somewhat dependent on the position of the site. Overall, the performances of Niño3.4 and SOI have opposite effects on the UIEP strength of the coming summer. Jiang et al. (2019) reported the similar opposite effects of Niño3.4 and SOI on extreme precipitation indexes (EPIs) [57].
To judge the overall correlation between regional UIEP and climate factors, the spatial patterns and temporal coefficients of UIEP were isolated using the EOF method, and the temporal coefficient corresponding to the first mode (PC1) was taken to perform a correlation analysis with UIEP. For the 1-12-month lagged UIEP versus Niño3.4 and SOI, the analysis results show that the correlation between PC1 and Niño3.4 is largest when the lagged time of UIEP is 5 months, and the correlation between PC1 and SOI is largets when the lagged time of UIEP is 1 month.

Cross-Wavelet Analysis
The cross wavelet transform method (XWT) is further adopted to discuss the nonlinear correlation between the UIEP and climate anomalies, and the periodic characteristics are revealed from the perspective of multiple time scales. According to the EOF analysis results, the cross-wavelet spectrum of the first principal component (PC1) of UIEP against Niño3.4 and SOI is shown in Figure 11.  Figure 8. the cross-wavelet spectrum of TS1 of the 5-month and 1-month lagged UIEP against Niño3.4 (a) and SOI (b). The thin black line represents the cone of influence. The thick enclosed areas denote the 5% significance level using red noise model based on Monte Carlo simulation method. Arrows signify the phase difference of oscillation. The arrow → means the in-phase relationship between the UIEP and climate anomalies, and ← means the anti-phase. The arrow ↓ means climate anomalies leading UIEP by 90°, and vice versa.
The significant sections shown in Figure 11 indicate a large covariance between both Niño3.4 and SOI against UIEP at a time scale of 2 years. The 5-month lagged UIEP versus Niño3.4 presents a significantly high common power period from 2007 to 2011, and Niño3.4 leads UIEP by about 135°. The 1-month lagged UIEP versus SOI shows significantly high common power periods from 2006 to 2012, while they are in the anti-phase.

Limitation and Uncertainties
Many studies have shown that precipitation response to urbanization has clear regional characteristics. At present, the statistical analysis of urban precipitation research based on measured Figure 11. The cross-wavelet spectrum of TS1 of the 5-month and 1-month lagged UIEP against Niño3.4 (a) and SOI (b). The thin black line represents the cone of influence. The thick enclosed areas denote the 5% significance level using red noise model based on Monte Carlo simulation method. Arrows signify the phase difference of oscillation. The arrow → means the in-phase relationship between the UIEP and climate anomalies, and ← means the anti-phase. The arrow ↓ means climate anomalies leading UIEP by 90 • , and vice versa.
The significant sections shown in Figure 11 indicate a large covariance between both Niño3.4 and SOI against UIEP at a time scale of 2 years. The 5-month lagged UIEP versus Niño3.4 presents a significantly high common power period from 2007 to 2011, and Niño3.4 leads UIEP by about 135 • .