Future Increase in Extreme Precipitation: Historical Data Analysis and Influential Factors

Abstract

With global warming driving an increase in extreme precipitation, the ensuing disasters present an unsustainable scenario for humanity. Consequently, understanding the characteristics of extreme precipitation has become paramount. Analyzing observational data from 1961 to 2020 across 29 meteorological stations in Heilongjiang Province, China, we employed kriging interpolation, the trend-free pre-whitening Mann–Kendall (TFPW–MK) method, and linear trend analysis. These methods allowed us to effectively assess the spatiotemporal features of extreme precipitation. Furthermore, Pearson’s correlation analysis explored the relationship between extreme precipitation indices (EPIs) and geographic factors, while the geodetector quantified the impacts of climate teleconnections. The results revealed the following: (1) There has been a clear trend in increasing extreme precipitation over the last few decades, particularly in the indices of wet day precipitation (PRCPTOT), very wet day precipitation (R95P), and extremely wet day precipitation (R99P), with regional mean trends of 10.4 mm/decade, 5.7 mm/decade, and 3.4 mm/decade, respectively. This spatial trend showed a decrease from south to north. (2) Significant upward trends were observed in both spring and winter for the maximum 1-day precipitation (RX1day) and the maximum 5-day precipitation (RX5day). (3) The latitude and longitude were significantly correlated with the most extreme precipitation indices, while elevation showed a weaker correlation. (4) Extreme precipitation exhibited a nonlinear response to large-scale climate teleconnections, with the combined influence of factors having a greater impact than individual factors. This research provides critical insights into the spatiotemporal dynamics of extreme precipitation, guiding the development of targeted strategies to mitigate risks and enhance resilience. It offers essential support for addressing regional climate challenges and promoting agricultural development in Heilongjiang Province.

1. Introduction

According to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change [1], the frequency of extreme climate events has increased in recent years due to climate change [1]. Extreme precipitation poses significant threats to human society, the ecological environment, and economic development. At present, it is generally believed that precipitation and temperature conform to the C–C (Clausius–Clapeyron) relationship, indicating that a 1 °C rise in surface temperature increases the atmospheric water vapor by approximately 6% to 7% [2]. It is anticipated that the frequency and intensity of extreme precipitation will likely escalate as climate change intensifies in the foreseeable future [3,4,5,6,7]. Therefore, understanding the dynamics and influencing factors of extreme precipitation is crucial for effective disaster management and predicting future extreme climate events.

Recently, the study of extreme climates has gained widespread attention, with researchers actively investigating these phenomena. Increased climate change affects the atmospheric cycle, with large amounts of water evaporating into the air as temperatures rise. Consequently, a warmer atmosphere holds more water vapor, triggering extreme events like floods, droughts, typhoons, and storms, seriously impacting agriculture, social stability, economic development, and residents’ lives [8,9,10]. Extreme precipitation, a prominent type of climate event, is influenced by factors such as sea and land position, atmospheric–oceanic circulation patterns, topography, monsoons, subsurface conditions, and human activities [11,12,13,14,15]. This leads to complex temporal and spatial variability of extreme precipitation. Numerous studies have shown that large-scale climate factors significantly impact extreme precipitation through a “teleconnection” mechanism [15,16,17,18,19]. Zhang [20] suggested that comprehending the variability of regional climate extremes requires circulation patterns to be accounted for. The East Asian Summer Monsoon (EASM) constitutes a vital element within the Asian monsoon complex, exerting profound influences on the meteorological conditions and climatic patterns across China and neighboring countries [21]. The El Nino/Southern Oscillation (ENSO) serves as a pivotal driver of climatic fluctuations intrinsic to the climate system [22]. Moreover, the Arctic polar vortex, East Asian trough, and atmospheric circulation at 850 hPa play significant roles in shaping the occurrence of extreme winter precipitation in Northeast China [23]. We are well aware that the contribution of large-scale climate remote correlation factors to extreme precipitation shows significant differences in different regions. Despite advancements, there remain challenges in fully understanding these complex interactions, particularly at regional scales.

Given its impact on the human living environment, scholars are deeply invested in studying the spatiotemporal distribution and causes of extreme precipitation to better understand these variations and effectively cope with future extreme climates. Global studies, including those in the United States [24], Canada [25], Europe [26], India [27], Japan [19], and the Asia monsoon region [28], have revealed varying trends in extreme precipitation at different regional scales. In China, the overall trend in extreme precipitation is consistent with global changes [29]. At the same time, various regions in China, such as the non-monsoon zone, the Hengduan Mountains, the Northeast, the Southwest, the Tibetan Plateau, the Yangtze River Basin, the Huanghuaihai Basin, the Pearl River, the Weihe River, and coastal areas, showed different patterns of extreme precipitation variability [30,31,32,33,34]. These studies demonstrate that the temporal and spatial variability of extreme precipitation is complicated by China’s vast geography and diverse climates. This underscores the necessity for more detailed, localized studies to capture regional peculiarities. It is more meaningful to conduct regional studies, which can better guide us to rationally respond to future climate change and take countermeasures in advance through region-specific studies. Although scholars have extensively studied extreme precipitation changes in northern China [35,36,37,38,39], the following shortcomings still exist in the Heilongjiang region. Previous studies have focused on specific regions or seasons, such as the correlation between winter extreme precipitation and atmospheric circulation in Heilongjiang [23], and variations in the Songhua River Basin [40], but lack long-term, region-wide analyses. Therefore, a comprehensive analysis of extreme precipitation in Heilongjiang Province is required to reveal the spatiotemporal distribution characteristics of extreme precipitation and causal mechanisms within the context of climate variability. Heilongjiang Province is an important commercial grain base in China, and dry farming is mostly based on rain-fed agriculture, with strong dependence on natural precipitation for crop growth. The inhomogeneity distribution of precipitation makes it have a greater impact on crop yields [41,42]. Understanding precipitation variations and mechanisms helps address water scarcity in rain-fed agriculture, as it is the primary water source. Addressing these gaps will provide valuable insights for sustainable agricultural practices and water resource management.

This research centers on Heilongjiang Province, an exemplar frigid zone within Northeast China, to analyze the spatiotemporal distribution attributes of extreme precipitation and its drivers. Employing the RClimDex model facilitated the calculation of extreme precipitation indices (EPIs), aiming to furnish a conceptual framework and technical support for the effective utilization of rainwater resources and the sustainable advancement of agricultural practices in the frigid zone of Northeast China. In addition, the outcomes of this study enrich our comprehension of regional mechanisms governing extreme precipitation occurrences and their dynamics.

2. Materials and Methods

2.1. Study Area

Heilongjiang Province lies in the northernmost and highest latitude of China (121°11′ E–135°05′ E, 43°26′ N–53°33′ N), bordered by Russia to the north, Inner Mongolia to the west, and Jilin Province to the south. Its total land area measures approximately 473,000 km², constituting 4.9% of China’s total land area, ranking it as the sixth largest province in China. The topography of Heilongjiang is diverse, with elevations higher in the northwest, north, and southeast, and lower in the Northeast and Southwest, encompassing mountains, plateaus, plains, and bodies of water. The province is traversed by the Heilongjiang, Wusuli, Songhua, and Suifen Rivers, and experiences a temperate continental monsoon climate. Annual precipitation ranges from 500 to 600 mm, with most rainfall occurring in the summer, and higher precipitation levels observed in the eastern regions compared to the western areas. Figure 1 shows the distribution of the study area and the meteorological stations.

2.2. Data Source

The daily meteorological records spanning from 1961 to 2020 were obtained from 37 meteorological stations situated within the research region. To ensure data consistency and integrity, stations with over 30 days of missing records were excluded, while those with fewer than 30 days of missing data were supplemented using data from nearby stations deemed most relevant. Following rigorous quality control measures, a total of 29 stations with comprehensive daily meteorological records were identified. The dataset was acquired from the China Meteorological Administration (http://data.cma.cn/, last accessed on 18 January 2021).

Taking into account the geographical and regional atmospheric circulation background [43,44,45], we opted to explore the relationship between regional EPIs and the corresponding large-scale climate teleconnections by selecting 14 atmospheric teleconnection indices. These indices include the Arctic Oscillation (AO), North Atlantic Oscillation (NAO), Atlantic Multidecadal Oscillation [27], ENSO Modoki Index (EMI), Multivariate ENSO Index (MEI), Pacific/North American Pattern (PNA), Tropic Indian Ocean Dipole (IOD), NINO 3.4 area sea surface temperature anomaly (NINO 3.4), NINO A area sea surface temperature anomaly (NINO A), 30 hPa Zonal Wind Index (30ZWI), 50 hPa Zonal Wind Index (50ZWI), Northern Hemisphere Subtropical High Intensity (NSI), Western Pacific Subtropical High Intensity (WPSH), East Asian Summer Monsoon Indices (EASMI). Data on these were obtained from the National Climate Center of the China Meteorological Administration (http://cmdp.ncc-cma.net/cn/index.htm, accessed on 1 May 2024).

2.3. Study Methods

2.3.1. Definition of Extreme Indices

The Expert Team on Climate Change Detection and Indices (ETCCDI) recommended 27 extreme climate indices, including 16 extreme temperature indices and 11 extreme precipitation indices (http://etccdi.pacificclimate.org, accessed on 1 May 2024). Compared with traditional methods such as anomaly values and precipitation anomaly percentages, the extreme climate indices defined by the ETCCDI are more reflective of the changes in extreme climate events; have weak extremeness, strong significance, low noise; and easily compare data from different regions. They are useful for decision makers to identify different types of extreme climate events. Therefore, this indicator system has been widely used by scholars in various countries to study the characteristics of extreme temperature and precipitation events. Considering the geographical and climatic conditions of the study area, 11 indices were chosen (see

Table 1. Definitions of extreme precipitation indices.

Index	Descriptive Name	Definition	Unit
R10mm	Number of heavy precipitation days	Annual count of days when RR ≥ 10 mm	days
R20mm	Number of very heavy precipitation days	Annual count of days when RR ≥ 20 mm	days
R25mm	Number of extremely heavy precipitation days	Annual count of days when RR ≥ 25 mm	days
R95P	Very wet day precipitation	Annual total precipitation when RR > 95th percentile	mm
R99P	Extremely wet day precipitation	Annual total precipitation when RR > 99th percentile	mm
CDD	Consecutive dry days	Maximum number of consecutive dry days	days
CWD	Consecutive wet days	Maximum number of consecutive wet days	days
PRCPTOT	Wet day precipitation	Annual total PRCP in wet days	mm
RX1day	Maximum 1-day precipitation	Annual maximum 1-day precipitation	mm
RX5day	Maximum 5-day precipitation	Annual maximum consecutive 5-day precipitation	mm
SDII	Simple daily intensity index	Average precipitation on wet days	mm/day

Notes: RR denotes daily precipitation amount. A wet day is defined when RR ≥ 1 mm and a dry day when RR < 1 mm.

) to depict extreme precipitation in Heilongjiang Province. These indices have been widely used in other studies evaluating extreme precipitation changes, ensuring that our analysis aligns with established methodologies and allows for regional comparisons. These specified indices reflect alterations in the occurrence, strength, and length of precipitation events and have gained extensive application in evaluating changes in extreme precipitation [20,46,47].

2.3.2. Line Tendency Test and TFPW–MK Trend Analysis

The trend in the EPIs was calculated using the ordinary least square (OLS) [48,49], in which the OLS uses a linear model to estimate the magnitude of the trend slope.

When analyzing the long time series of meteorological data in hydrometeorology, parametric and nonparametric methods are usually used to test their trend changes. Among them, the Mann–Kendall (MK) test [50,51], widely recommended for time series analysis, is extensively used. However, as research has progressed continuously, the method has also revealed some shortcomings [52]; it is believed that the autocorrelation of the original data could lead to the enhancement or degradation of the detection results to a certain extent, so the PW–MK (pre-whitening–MK) method has been proposed. Refs. [53,54] suggested that the pre-whitening treatment could reduce the ability of MK to detect the significance results; ref. [55] found that by proving for an autocorrelation sequence with trend terms, the pre-whitening treatment removes some of the trend terms, leading to the acceptance of the null hypothesis, and then proposed another method of de-trending pre-whitening, namely, the TFPW–MK (trend-free pre-whitening MK) method. In this study, the TFPW method is used to preprocess the EPIs to remove autocorrelation in the original time series [56], and the trend test is performed on the new series after TFPW processing. The specific procedure of the TFPW method is as follows:

(1)

Calculate the linear trend $\beta$ of each index series $X_t$(t = 1, 2, …, N):

$$\beta =Median\left({\displaystyle \frac{X_j-X_i}{j-i}}\right)\forall i<j$$

(1)

(2)

Form a new series $Y_t$ by removing the trend component from the series:

$$Y_t=X_t-\beta \cdot t$$

(2)

(3)

Calculate the first-order autocorrelation coefficient $r$ of the new series $Y_t$ and perform a significance test on $r$ (the significance test $p$ is 0.1). If the test is successful, directly apply the Mann–Kendall (MK) method to test the original series $X_t$. If not, proceed to step (4) for preprocessing.

(4)

Form a new series $Y_t^{\prime }$ by removing the autocorrelation terms from the series, and then reintroduce the trend component to create a new series $Y_t^{″}$ that is not affected by autocorrelation interference:

$$Y_t^{\prime }=Y_t-r\cdot Y_{t-1}$$

(3)

$$Y_t^{″}=Y_t^{\prime }+\beta \cdot t$$

(4)

(5)

Substitute the new series $Y_t^{″}$ into the MK test, which involves the following calculation steps:

For the time series variables $X=\left\{x_1,x_2,\cdots ,x_n\right\}$, where $n$ is the length of the time series, the test statistic $S$ computed as

$$S={\displaystyle \sum _{i=1}^{n-1}{\displaystyle \sum _{j=i+1}^{n}sign\left(x_j-x_i\right)}}$$

(5)

$$sign\left(x_j-x_i\right)=\left\{\begin{array}{l}1,x_j-x_i>0\hfill \\ 0,x_j-x_i=0\hfill \\ -1,x_j-x_i<0\hfill \end{array}\right.$$

(6)

where $x_i$ and $x_j$ denote the values of year $i$ and year $j$ of the time series, respectively. $sign$ is the sign function, and $S$ the variance is calculated as

$$Var={\displaystyle \frac{n\left(n-1\right)\left(2n+5\right)}{18}}$$

(7)

A standardized test statistic $Z$ can be constructed according to the following equation:

$$Z=\left\{\begin{array}{l}{\displaystyle \frac{S-1}{\sqrt{Var\left(S\right)}}},S>0\hfill \\ 0,S=0\hfill \\ {\displaystyle \frac{S+1}{\sqrt{Var\left(S\right)}}},S<0\hfill \end{array}\right.$$

(8)

The significance of a trend is measured by calculating a statistical parameter. When positive, the trend increases, and when negative, the trend decreases. Statistical significance was tested at the 95% confidence level (p < 0.05) for ≥1.96 and at the 99% confidence level (p < 0.01) for ≥2.58 [57].

2.3.3. Spatial Interpolation

Meteorological data are characterized by discontinuity and inhomogeneity because they are experimentally collected and organized from meteorological stations throughout the country. In scientific research, the method of spatial interpolation is often needed to obtain data from regions other than stations [58,59,60]. This study used kriging interpolation to spatially interpolate the EPIs calculated from 29 meteorological stations in Heilongjiang Province from 1961 to 2020.

2.3.4. Correlation Analysis

Pearson’s correlation analysis is a statistical method utilized to quantify the strength and direction of the linear relationship between two continuous variables [48,57,61]. It provides valuable insights into the relationships between variables [35,62,63,64]. The Pearson correlation coefficient, denoted as $r$, ranges from −1 to 1, where $r=1$ indicates a perfect positive linear relationship, $r=-1$ indicates a perfect negative linear relationship, and $r=0$ indicates no linear relationship. Pearson’s correlation analysis investigated the relationship between EPIs and revealed its relationship with the latitude, longitude, and altitude of the study area.

2.3.5. Influencing Factor Quantification

The geodetector method (GDM) is a comprehensive measurement technique grounded in the principles of analysis of variance (ANOVA) for identifying stratified heterogeneity (SH) and elucidating its driving factors. The fundamental premise of geodetector is that if an independent variable significantly impacts a dependent variable, their spatial patterns should exhibit similarity [65]. The method is described as follows:

(1)

Factor detection: Factor detection mainly examines the strength of the driving force using a statistical measure denoted as $q$. The calculation is as follows:

$$q=1-{\displaystyle \frac{1}{N{\sigma }^2}}{\displaystyle \sum _{h=1}^L{N}_h{\sigma }_h^2=1-\frac{SSW}{SST}}$$

(9)

$$SSW={\displaystyle \sum _{h=1}^L{N}_h{\sigma }_h^2}$$

(10)

$$SST=N{\sigma }^2$$

(11)

where $h$ represents the number of categories or partitions of $X$ or $Y$; $N$ and $N_h$ denote the number of units in the entire area and within layer $h$, respectively. The variances ${\sigma }^2$ and ${\sigma }_h^2$ are for the entire area and layer $h$ regarding the dependent variable $Y$, $SSW$ is the sum of variances within layers, while $SST$ is the total variance of the entire area. The value $q$ indicates the contribution rate of the independent variable $X$ to the dependent variable $Y$, ranging from 0 to 1.

(2)

Interaction detection: Interaction detection can be used to determine the type of relationship between two factors and their impact on EPIs. The specific interactions are detailed in

Table 2. The manner in which two factors interact.

Comparative Result	Interaction Type
$q\left(X_1\cap X_2\right)<Min\left(q\left(X_1\right),q\left(X_2\right)\right)$	Nonlinear weakening
$Min\left(q\left(X_1\right),q\left(X_2\right)\right)<q\left(X_1\cap X_2\right)<Max\left(q\left(X_1\right),q\left(X_2\right)\right)$	Univariate weakening
$q\left(X_1\cap X_2\right)>Max\left(q\left(X_1\right),q\left(X_2\right)\right)$	Two-factor enhancement
$q\left(X_1\cap X_2\right)=q\left(X_1\right)+q\left(X_2\right)$	Mutually independent
$q\left(X_1\cap X_2\right)>q\left(X_1\right)+q\left(X_2\right)$	Nonlinear enhancement

.

3. Results and Analysis

3.1. Spatiotemporal Variation in EPIs

The results of the linear trend and statistical analyses of the selected 11 EPIs are shown in

Table 3. Trend analysis of EPIs in Heilongjiang Province from 1961 to 2020.

Indices	Units	Range of Regional Trends (Mean) (Decades−1)	Decreasing Trend (SS)	Increasing Trend (SS)	No Trend
R10	days	0.01~0.88 (0.3)	0 (0)	29 (1)	0
R20	days	−0.14~0.57 (0.2)	3 (2)	26 (1)	0
R25	days	−0.13~0.57 (0.1)	4 (4)	25 (2)	0
R95P	mm	−2.62~17.14 (5.7)	2 (2)	27 (0)	0
R99P	mm	−5.15~11.14 (3.4)	4 (4)	25 (0)	0
CDD	days	−6.99~1.29 (0.2)	25 (25)	4 (0)	0
CWD	days	−0.31~0.30 (0)	10 (8)	15 (8)	4
PRCPTOT	mm	1.10~32.36 (10.4)	0 (0)	29 (0)	0
RX1day	mm	−1.24~6.48 (0.5)	5 (5)	24 (1)	0
RX5day	mm	−1.31~6.80 (0.8)	6 (6)	23 (1)	0
SDII	mm/day	−0.10~0.34 (0.1)	4 (3)	25 (0)	0

Note: The numbers in columns 4, 5, and 6 indicate the number of stations; SS denotes the number of stations with trend in extreme precipitation indices passing the 0.05 confidence level.

and Figure 2. The spatial distribution pattern of the decadal rate of change in the 11 indices in Heilongjiang Province from 1961 to 2020 is shown in Figure 3. The spatial variations in the indices are also analyzed, as shown in Figure 4.

Except for the CWD, which did not change significantly, all 10 EPIs showed an increasing trend, with R10mm, R20mm, and R25mm increasing at a rate of 0.3 days decade⁻¹, 0.2 days decade⁻¹, and 0.1 days decade⁻¹, respectively (Figure 2a–c). Linear trends for regional means ranged from 0.01 to 0.88 days decade⁻¹, −0.14 to 0.57 days decade⁻¹, and −0.13 to 0.57 days decade⁻¹, respectively (Table 3). The 5-year smoothing average showed a small upward trend in the three indices after the 1980s (Figure 2a–c). One of the spatial trends in R10mm showed an increasing trend at 29 stations, but only one showed a significant increasing trend (p < 0.05) (Table 3), and none of the stations showed a decreasing trend at the 95% confidence band. These stations with increasing trends in R10mm were mainly located in the central part of the study area, and stations with higher rates of increasing R10mm were found in the southeastern corner (Figure 3a). The spatial pattern of R10mm indicated that the eastern region had longer days of heavy rainfall than the western region.

The spatial trend of R20mm showed that 26 stations showed an increasing trend and only 1 station showed a significant increasing trend (p < 0.05), while 3 stations showed a decreasing trend and 2 stations showed a significant decreasing trend (Table 3). Stations showing an increasing trend were predominantly located in the eastern region, and stations showing a decreasing trend were predominantly located in the central region (Figure 3b). In addition, the results of the R25mm linear regression and MK significance test showed that 4 stations showed a significant decreasing trend (p < 0.05) and 2 stations showed a significant increasing trend out of 25 stations with an increasing trend (Table 3). Both stations with significant increasing and decreasing trends were distributed in the central location (Figure 3c).

From (Figure 2d,e), it can be observed that the regional average trends in the R95P and R99P from 1961 to 2020 were 5.7 mm decade⁻¹ and 3.4 mm decade⁻¹, respectively, which showed a large increasing trend, and the range of the linear trend in R95P was −2.62 to 17.14 mm decade⁻¹ and for R99P was −5.15 to 11.14 mm decade⁻¹, respectively. The 5-year smoothing average showed a sharp upward trend after 2010 (Figure 2d,e). However, the results of the linear trend and MK significance test showed that R95P did not show a significant trend at any of the 27 increasing trend stations, 4 stations showed a decreasing trend, and all of them showed a significant trend (p < 0.05). There were two significant decreasing trend stations for R99P, and none of the 27 increasing trend stations showed a significant increasing trend (Table 3). For the spatial distribution of the R95P index (Figure 3d), most of the stations showed a greater increasing trend in the eastern and western regions of Heilongjiang Province. In addition, two stations with significant trends were in the northwest and southwest corners. The spatial distribution of the R99P index was similar to that of the R95P, and it was also found that stations with larger upward trends were mainly concentrated in the western region, whereas those with larger downward trends were located in the central part of the area, and three stations with significant decreases were located in the southern corner (Figure 3e).

CDD exhibited a slight upward trajectory, with a regional average increase of 0.2 days decade⁻¹ and a linear trend of −6.99–1.29 days decade⁻¹ (Table 3). Following the 1990s, there was a marginal increase in CDD as depicted by the 5-year smoothed average (Figure 2f). Notably, the spatial distribution trend was more evident, as only 4 stations displayed an ascending trend, while none showed a statistically significant positive trend. Conversely, a decreasing trend was observed in 25 stations, all of which demonstrated statistical significance (p < 0.05) (Table 3). These declining trends in CDD encompassed nearly the entire study area, with stations exhibiting larger ranges primarily situated in the northwestern and southern regions (Figure 3f). The spatial configuration of CDD indicated that the eastern portion of the study area experiences lengthier periods of consecutive dry days compared to the western.

Compared with the changes in the CDD, the long-term regional average changes in the CWD index were inconspicuous, with the linear trend extending from −0.31 to 0.30 days decade⁻¹ (Table 3). As shown in Figure 2g, the CWD showed a rapid and sustained downward trend from the 1960s to the 1980s, followed by an upward trend. Ten stations showed a negative trend whereas fifteen stations showed a positive trend (Table 3), with eight stations showing a significant downward trend and eight stations showing a significant upward trend (p < 0.05) within the 95% confidence level. At the same time, four stations showed no trend. The statistics showed only stations where the CWD index showed no trend, while the rest of the indices showed increasing or decreasing trends (Table 3). The spatial distribution pattern map of the CWD changes (Figure 3g) revealed that the stations exhibiting an upward trend were predominantly situated in the western region, while those displaying a downward trend were mainly found in the eastern and northern areas. This indicated that the western portion experiences longer periods of consecutive wet days compared to the eastern part.

The findings indicated that the linear trend in the PRCPTOT varied between 1.10 and 32.36 mm decade⁻¹ over the past 60 years, with a local mean of 10.4 mm decade⁻¹ (Table 3 and Figure 2h). With reference to the 5-year smoothing average (Figure 2h), a minor decline in the PRCPTOT was observed from the early 1960s to the 1980s, which then showed a rapid increasing trend with large variations, followed by a slow decreasing trend. The 29 stations exhibited an upward trend, yet according to the MK test, none of these stations displayed a statistically significant increasing trend (p < 0.05). As shown in Figure 3h, notable increases were observed in both the western and eastern regions, whereas the central region depicted a comparatively minor uptrend.

The linear trends for RX1day and RX5day were in the range −1.24–6.48 mm decade⁻¹ and −1.31–6.80 mm decade⁻¹, respectively. Despite the relatively large variations in the different regions, the regional averages of RX1day and RX5day showed an increasing trend of 0.5 mm decade⁻¹ and 0.8 mm decade⁻¹ over the period 1961–2020 (Table 3). In terms of the regional average annual trends, the two indices showed similar trends (Figure 2i,j). The 5-year smoothing average line showed a slight long-term downward change from the 1960s to the 1980s, followed by a sharp increase after 2010 (Figure 2i,j). Table 3 additionally displays the statistically significant upward (positive), downward (negative) trends and without discernible trends. Out of all of the RX1day stations, five exhibited a significant negative trend (p < 0.05), whereas another twenty-four showed positive trends and only one of all the stations was significant. Similarly to RX1day, 6 of the stations studied for RX5day showed a negative trend, and all of the stations were significant (p < 0.05), whereas another 23 showed a positive trend, with only 1 station being significant (Table 3). The spatial distribution of the linear trends for RX1day and RX5day are shown in Figure 3i,j. In terms of the spatial distribution, RX1day and RX5day also had similar spatial distribution characteristics, and in general, the rate of change in RX1day and RX5day was larger in the western regions, and most of the stations were significantly downward. Meanwhile, the stations with larger ranges of increasing trends in RX1day and RX5day were also distributed in the western region.

Similarly to most other indices, the extended-term variation in SDII during the preceding 60 years indicated a rise (Figure 2k), spanning regionally from −0.10 to 0.34 mm day⁻¹ decade⁻¹ (averaging 0.1 mm day⁻¹ decade⁻¹). Concerning the 5-year smoothing average trend, SDII fluctuated similarly to the RX1day and RX5day regional trends, declining continuously from the 1960s until the 1980s, with little change in the middle of 1980–2010, followed by a sharp upward trend. Concerning trends at individual stations (Table 3), 25 stations exhibited an upward trajectory, yet none reached the significance threshold at the 0.05 level. Conversely, four stations displayed a downward trend, with statistical significance noted at three stations (p < 0.05). Notably, the stations experiencing substantial upward trends were concentrated in the western and eastern regions, while those with minor negative trends predominated in the central region (Figure 3k).

In addition, the results of the kriging interpolation analysis based on the calculated EPIs indicated that the multi-year mean extreme precipitation in the study area was spatially heterogeneous (Figure 4a–k). Higher values of R10, R20, and R25 were observed in the central part of Heilongjiang Province, roughly 15–18 days, 5–7 days, and 3–5 days, respectively, whereas lower values of about 13 days, 4 days, and 2 days, respectively, were observed in the Daxing’anling region, where the main extreme precipitation occurred in the Songhua River basin (the specific location can be referred to in [40]). The two indices, R95P and R99P, had similar distributional characteristics, with the highest intense precipitation mainly in the central and eastern regions, amounting to 130–150 mm and 40–50 mm, respectively, and 110 mm and 35 mm in the Daxing’anling region, where the intense precipitation was lower than that in the other regions of the study area.

The interpolated results of CDD and CWD showed that it is driest in parts of the Songnen Plain and wettest in the Sanjiang Plain in the central and eastern parts of the study area. Then, for the total precipitation PRCPTOT results, it was still the central and eastern regions that showed the highest precipitation, between about 480 and 590 mm, especially in Harbin, where the multi-year average precipitation total is more than 550 mm, and in the Daxing’anling region and part of the Songnen Plain, the total precipitation is the smallest, at only 410–450 mm, and this part of the region is semiarid. The RX1day and RX5day values are low in the Daxing’anling region, being 45 mm and 75 mm, respectively, and the maximum daily precipitation and 5-day maximum precipitation in the Songhua River basin are 61.28 mm and 99.78 mm, respectively. SDII was approximately 8–9.5 mm/day in the Songnen Plain region with an average daily precipitation intensity, and was 6.45–7.06 mm/day in the Daxing’anling region.

3.2. Seasonal Variation in EPIs

Due to the RClimdex providing solely monthly data for RX1day and RX5day, this study exclusively examines the seasonal fluctuations in these two indices. Figure 5 illustrates the seasonal variation outcomes for RX1day and RX5day. RX1day illustrated a clear trend of growth in all the seasons, with the largest trend being 2.4 mm decade⁻¹ in summer and 2.0 mm decade⁻¹ in autumn. The values for spring and winter were 1.3 mm decade⁻¹ and 0.8 mm decade⁻¹, respectively (Figure 4). For RX5day, there was also an upward trend in all seasons, with a 3.8 mm decade⁻¹ upward trend in summer, 2.1 mm decade⁻¹ in spring, and 1.6 mm decade⁻¹ in autumn, and a minimum upward trend of 1.4 mm decade⁻¹ in winter. The trend analysis of the two indices revealed significant changes at the 0.05 significance level in spring and the 0.01 significance level in winter, indicating that the extreme precipitation indices changed significantly on a seasonal basis. The summer extreme precipitation in Heilongjiang Province has shown a pronounced upward trend in the RX1day and RX5day indices, while the other seasons display only modest increases. This uneven distribution of seasonal precipitation suggests an increase in summer rainfall and a decrease in rainfall during other seasons. Consequently, this pattern may elevate the risk of summer floods and droughts in other seasons [66]. Heilongjiang Province resides in the chilly domain of Northeast China, frequently experiencing snowstorms. Therefore, research on seasonal changes in EPIs can provide some basis for local agricultural development.

3.3. Relationship Between Extreme Precipitation Indices and Total Precipitation

The precipitation intensity, frequency, and duration are important factors affecting the total precipitation [67]. We analyzed the main EPIs that influenced the increasing precipitation trend in Heilongjiang Province from 1961 to 2020.

The correlation between the EPIs and total precipitation (PRCPTOT) is depicted in Figure 6. Apart from CDD, the remaining indices exhibited a positive correlation with the total precipitation (p < 0.01). Notably, the correlation coefficients between R10mm, R20mm, R25mm, and PRCPTOT were 0.97, 0.97, and 0.93, respectively, surpassing those of the other indices. These findings highlight that the augmentation of R10mm, R20mm, and R25mm is the primary reason for the escalation in the regional total precipitation in Heilongjiang Province from 1961 to 2020. Other indices are also related to some extent, but have no correlation with the change in CDD.

3.4. Relationship Between EPIs and the Regional Geographic Factors

Pearson’s correlation analysis was employed to further investigate the influence of regional geographic elements (longitude, latitude, and elevation) on the changes in EPIs in Heilongjiang Province (

Table 4. Correlation between extreme precipitation indices and longitude, latitude, and altitude.

Indices	Longitude	Latitude	Altitude
R10mm	0.58 b	−0.37 a	0.10
R20mm	0.29	−0.34	−0.02
R25mm	0.22	−0.35	−0.06
R95P	0.58 b	−0.45 a	0.08
R99P	0.57 b	−0.49 b	0.16
CDD	−0.70 b	0.05	−0.16
CWD	0.38 a	−0.13	0.41 a
PRCPTOT	0.63 b	−0.39 a	0.10
RX1day	−0.01	−0.58 b	−0.09
RX5day	−0.08	−0.39 a	0.01
SDII	−0.37 a	−0.32	−0.21

Note: ^a exhibits significance at 0.05, while ^b demonstrates significance at 0.01.

).

The correlation analysis revealed a significant association between the EPIs and latitude as well as longitude. While R20mm, R25mm, RX1day, and RX5day did not exhibit a significant correlation with the longitude, other variables demonstrated a significant correlation. Specifically, CWD displayed a significant positive correlation (p < 0.05), whereas SDII showed a significant negative correlation (p < 0.05). Moreover, R10mm, R95P, R95P, and PRCPTOT exhibited a significant positive correlation with the longitude (p < 0.01), whereas CDD showed a significant negative correlation (p < 0.01). These findings suggested that the longitude is an important factor influencing extreme precipitation in Heilongjiang, China.

Generally, higher latitudes correspond to lower solar angles and reduced solar energy reaching the Earth’s surface. From 1961 to 2020, there was a positive correlation between CDD and latitude, while the other exhibited negative correlations. Moreover, with the exception of R20mm, R25mm, CWD, and SDII, most of the EPIs showed significant correlations, with R99P and RX1day demonstrating significant negative correlations (p < 0.01). Figure 3 illustrated that stations exhibiting substantial variation trends are predominantly situated in the south, while those with minimal gradient trends are clustered in the north. These findings indicated a close association between the latitude and extreme precipitation, with the trend in extreme precipitation gradually decreasing from south to north. Notably, Heilongjiang Province is situated in the region with the highest latitude in China. The solar radiation energy is low because of the higher latitude, which determines the heat belt of climate and the high and low distribution of temperature, and has a greater impact on precipitation.

As shown in Table 4, the majority of the indices displayed weak correlations with the altitude, with the exception of CWD. CWD exhibited a significant positive correlation, with a correlation coefficient of 0.028 (p < 0.05). R10mm, R95P, R99P, PRCPTOT, and RX5day were positively correlated with the altitude, whereas R20mm, R25mm, CDD, RX1day, and SDII were negatively correlated. The results showed that extreme precipitation mainly occurs in low-altitude areas.

3.5. Influence of Large-Scale Climate Factors on Extreme Precipitation

3.5.1. Contribution of Single Large-Scale Climate Factor

The study calculated q-values to delineate the extent of each factor’s contribution to the EPIs, with q acting as detectors in GDM. The magnitude of the contribution of the large-scale climate factors to the EPIs is shown in Figure 7. While variations exist in the impact of large-scale climatic teleconnections factors on the EPIs, notable contributors to the EPIs include AMO, MEI, IOD, Nino 3.4, NSI, WPSH, and EASMI. Among these, AMO and NSI exerted the most significant influence, with an average contribution rate of 15%, while MEI, WPSH, EASMI, IOD, and Nino 3.4 averaged 14%, 14%, 12%, 11%, and 8%, respectively. Notably, for the duration indices, AMO remained a primary factor, exerting the largest effects on R20mm (18%) and R25mm (18%). MEI emerged as the crucial climate factor influencing R10mm (21%). CDD and CWD were predominantly influenced by Nino A (21%) and EASMI (16%). Regarding the frequency indices (R95p, R99p, and PRCPTOT), AMO (17%), MEI (21%), and NSI (21%) have the largest shares. MEI (19%), IOD (19%), and AMO (17%) exhibited the greatest influence on the intensity indices (RX1day, RX5day, and SDII).

3.5.2. Influence of the Interaction of Large-Scale Climate Factors

Extreme precipitation is actually influenced by a combination of multiple factors rather than a single factor alone. To evaluate the influence of interactions between atmospheric circulation factors on the EPIs in this region, the “Interaction detection” function of GDM was used to gauge the extent of the impact from two-factor interactions. It was found that both two-factor enhancement (q(X1 ∩ X2) > Max(q(X1),q(X2))) and nonlinear enhancement (q(X1 ∩ X2) > q(X1) + q(X2)) were observed.

The nonlinear enhancement predominated in the interactions between factors. However, according to the analysis of interaction detection, it was found that the NSI ∩ WPSH interactions exhibited two-factor enhancement for R10, R20, R25, CDD, R95p, R99p, PRCPTOT, RX1day, RX5day, and SDII, while the Nino 3.4 ∩ EASMI interactions displayed bivariate enhancement for R10, R20, R25, R95p, R99p, PRCPTOT, RX5day, and SDII. The same interaction produces different contributions to different EPIs. Figure 8 illustrates that NSI ∩ MEI had explanatory powers of 59%, 62%, 62%, and 63% for R20, R25, R95p, and PRCPTOT, respectively. However, the explanatory capacity of NSI ∩ MEI regarding CWD stood at 38%. Even though NSI had the most substantial influence on the EPIs, its mean contribution rate amounted to merely 15%. This underscored the significance of the combined effect of multiple factors in influencing the EPIs compared to any single factor alone.

It was observed that Nino A emerged as the predominant factor in the interaction, contrasting with the previous conclusion that AMO and NSI were the primary factors in the univariate interaction. The majority of factors interacting with Nino A exhibited higher contribution for changes in the EPIs.

4. Discussion

4.1. Evident Spatiotemporal Variability in Extreme Precipitation

The temporal and spatial distribution patterns of extreme precipitation vary regionally, and numerous factors contribute to the occurrence of such events, including geographic factors, atmospheric circulation, surface cover, climate change, and human activities [43,68,69,70]. Researchers have selected different regions and climate models to explore the spatiotemporal distribution of extreme precipitation and quantitatively analyzed various influencing factors to reveal their intrinsic connection with extreme precipitation. Our investigation was conducted using a lengthy temporal dataset spanning six decades, which indicated a rising tendency in the EPIs within Heilongjiang Province. The exception was the CDD index, where the interannual trend remained inconspicuous; however, all the other indices exhibited an upward trajectory. Among them, the trends in RX1day and RX5day were significant on the seasonal scale. Zhang [23] conducted an analysis on the alterations in wintertime extreme precipitation within Heilongjiang Province. The findings revealed that the average annual precipitation was recorded at 13.7 mm during the period 1961–2018, with a corresponding increase rate of 1.8 mm/10a. Previous research has indicated a rise in the intensity of RX5day in the majority of regions within the Songhua River basin, with approximately 76.9% and 61.5% of observation stations showing noteworthy (p < 0.05) increments in April and October, respectively [40]. However, within the scope of China, the precipitation distribution and its heterogeneity exhibit notable variations, while extreme precipitation manifests distinct spatial distribution patterns across different regions. An examination of the temporal and spatial variations in extreme climate events in northeastern China revealed that stations exhibiting CDD decline rates exceeding 2 day/decade are exclusively situated in Heilongjiang Province. Furthermore, significant variations were observed in PRCPTOT, with over 70% of stations demonstrating an increasing trend, paralleled by similar trends in R95P, which dominates precipitation in these regions. Notably, the upward trend in R95P is primarily concentrated in the southwest of Liaoning and eastern Heilongjiang [71], consistent with our findings of a declining trend in 25 stations, all of which displayed significant decreases (p < 0.05). Additionally, our analysis highlighted similar distribution characteristics between the indices R95P and R99P, with the most intense precipitation occurring predominantly in the central and eastern regions. The assessment of the spatiotemporal distribution of extreme climate hazards, notably droughts and floods, revealed a pronounced upward trend from 2006 to 2019. The spatial analysis of precipitation in the Northeast indicated a decreasing trend from south to north, with relatively minimal variation from east to west [38].

Extreme precipitation’s regional traits can be gauged by its intensity, frequency, and duration. In the context of global warming, there has been a general tendency for both the frequency and intensity of extreme precipitation to rise, yet significant regional variations exist [72]. Research findings based on observational data and CMIP5 model simulations indicated a sharp increasing trend in both total precipitation and extreme precipitation in arid regions worldwide, as well as extreme precipitation in humid regions, over the past 60 years [3]. For every 9 °C increase in surface warming in China, precipitation and extreme precipitation increase by 7.22% and 6.1%, respectively [67], a result larger than the global average, indicating that China’s regional rainfall is particularly responsive to climate change, underscoring the necessity of investigating how different regions respond to climate warming. Global warming is potentially leading to a decrease in snowfall in mountainous regions, which could shift to extreme precipitation [7]. The increase in snowfall and other precipitation types suggests that while global warming may reduce snowfall in some regions, it can also intensify precipitation in others, leading to more frequent and intense extreme precipitation events. Heilongjiang Province, located in the cold Northeast, is showing an increase in snowfall, which also contains other forms of water, such as snowfall and hail, and thus an increase in the precipitation intensity. Amid the ongoing large-scale warming trend, long-term precipitation patterns exhibit a significant rise in Heilongjiang Province, aligning with findings from prior studies in different areas [34,38,73,74,75]. Within the arid Northwest China region, heightened precipitation intensity and frequency contribute to a marked escalation in extreme precipitation occurrences [43,76]. The findings from the study indicated a rising trend in the total precipitation within Heilongjiang Province. The correlation analysis revealed that R10mm, R20mm, and R25mm exhibited the most significant correlations with PRCPTOT (p < 0.01) (Figure 6). The heightened total precipitation primarily resulted from increased daily precipitation. Over the course of continuous climate change, the northeastern region of China has experienced gradual warming and increased humidity, leading to a steady upward trend in precipitation.

4.2. Potential Drivers of Extreme Precipitation Variability

The spatial variability of extreme precipitation also exhibits variance owing to the impact of East Asian monsoons. The effect of topography on extreme precipitation in Northeast China has become evident in recent years [39]. National-scale investigations in China regarding the relationship between extreme precipitation and terrain reveal a consistent observation: altitude, latitude, and longitude demonstrate significant correlations with extreme precipitation [33,74,77,78]. The EPIs showed a significant correlation with both latitude and longitude, indicating a clear influence of these geographical factors on extreme precipitation events. This phenomenon could be attributed to the warming climate, which accentuates the contrast in thermal characteristics between land and sea, thereby amplifying extreme precipitation [1]. However, the correlation between the EPIs and altitude lacks significance; it is noteworthy that extreme precipitation predominantly occurs at lower altitudes. Therefore, we believe that extreme precipitation in low-lying regions exhibits greater sensitivity to climate change. The results of the correlation analysis were not significant because the altitude range of the selected stations (66.4–567.8 m) spanned a small area. In contrast, in the non-monsoon region of China, researchers analyzed the significant positive correlation between stations distributed above 3500 m in altitude and most indices in the Tibetan Plateau and surrounding areas from 1961 to 2017 [43]. In the Pearl River Basin, except for CDD, all other indices were negatively correlated with the altitude, and extreme precipitation events mainly occurred at low altitudes [79]. The distinct cyclical patterns of individual atmospheric circulations can potentially influence the characteristics of regional precipitation variations [80]. We came to a conclusion consistent with those of the majority of studies: the combined influence of two factors contributes more significantly to EPIs than individual factors [20,81,82]. Furthermore, a closer examination of the mechanisms driving the impact of atmospheric circulation on extreme precipitation reveals the significant influence of various phenomena, including AMO, MEI, IOD, Nino 3.4, NAO, WPSH, and EASMI. Research has suggested that these interactions can lead to significant changes in weather patterns, influencing both the intensity and frequency of extreme events. For instance, the NAO is known to affect storm tracks and the precipitation distribution in the Northern Hemisphere [83]. Similarly, the ENSO phenomena, captured by MEI and Nino 3.4, are crucial in modulating precipitation patterns globally [84]. The interplay and modulation of regional atmospheric dynamics by these large-scale factors dictate the variability in extreme precipitation. However, the complex interplay between large-scale atmospheric circulation patterns and regional climate dynamics serves as a fundamental driver of extreme precipitation variability. Understanding these complex mechanisms is essential for improving our capacity to forecast and alleviate the effects of extreme precipitation in a changing climate landscape.

4.3. Exploring Uncertainty and Future Prospects

In water resource management and climate change research, meteorological station-based datasets are prioritized [85]. However, to comprehensively analyze extreme precipitation in specific regions, it is essential to employ multiple datasets or models, particularly those that offer long-term, high-precision, and high-resolution reanalysis grid data. This facilitates a deeper understanding of precipitation patterns from a global perspective [86]. Future studies should consider integrating comparative results from different data sources to more accurately highlight the characteristics of extreme precipitation changes in Heilongjiang Province. The driving force behind extreme precipitation change is remarkably intricate due to a multitude of factors influencing regional precipitation variability. These factors include fluctuations in the climate system, atmospheric circulation, regional environmental characteristics, human activities, and more. In this study, we concentrate exclusively on the impacts of geographic factors and interconnected large-scale climate teleconnections. Therefore, future investigations may entail a comprehensive analysis of physical mechanisms. In summary, with the intensification of climate warming, researchers have thoroughly examined the mechanisms behind extreme precipitation from diverse perspectives. These findings can provide valuable insights for both human adaptation and future predictions of extreme precipitation. The study’s focus on Heilongjiang Province, as an important commercial grain base in China, offers guidance for local agricultural production and life. In addition, as a typical representative of the cold region in Northeast China, the research methods and results provide a foundation for understanding extreme precipitation in cold regions, contributing to the broader academic discourse on climate change impacts.

5. Conclusions

This study aimed to explore the spatiotemporal patterns of extreme precipitation, using 11 selected indices to potentially uncover mechanisms underlying the impact of extreme precipitation. The following main conclusions were drawn from this work:

(1)

Most of the extreme precipitation indices showed increasing trends and these increasing trends account for about 80% of the stations. However, only CWD was significantly increased in about half of the stations (p < 0.05). Furthermore, RX1day and RX5day exhibited obvious increasing trends across all four seasons, with the most pronounced trends observed during spring and winter.

(2)

The increasing trends in extreme precipitation are mainly distributed in the south, whereas the downward trends are distributed in the northern part of the study area; extreme precipitation is mainly distributed in the central and eastern regions, with a higher concentration of precipitation in the Songhua River basin.

(3)

The increase in the total precipitation in Heilongjiang Province from 1961 to 2020 was primarily driven by R10mm, R20mm, and R25mm, whose correlation coefficients reached 0.97, 0.97, and 0.93, respectively. In contrast, CDD showed no correlation with the total precipitation.

(4)

The altitude exhibited a weak correlation with the EPIs, while the latitude and longitude demonstrated significant correlations. Meanwhile, an analysis identified that the AMO, MEI, IOD, Nino 3.4, NSI, WPSH, and EASMI were the most important factors affecting extreme precipitation. The importance of considering the combined impact of multiple factors on the EPIs, rather than focusing solely on any single factor, was highlighted.