Spatiotemporal Variation in the Meteorological Drought Comprehensive Index in the Beijing–Tianjin–Hebei Region during 1961–2023

Abstract

It is crucial to investigate the characteristics of meteorological drought in the Beijing–Tianjin–Hebei (BTH) region to improve the accuracy of agriculture and water resource monitoring and management. In this study, using instrumental observation data from 85 meteorological stations in the BTH region during 1961–2023 derived from the National Meteorological Information Center, we first calculated the meteorological drought comprehensive index (MCI) and analyzed the spatiotemporal characteristics of the MCI. In the BTH region, the MCI intensity from May to June was the most severe in the intraseasonal variation. The trend of the decreasing drought intensity in May–June has occurred for the past 60 years. The southern region in the BTH region was more likely to experience droughts. Next, the spatial patterns of the top two EOF modes of the May–June MCI were depicted. The primary spatial pattern of the BTH, which was characterized by consistent changes in the MCI throughout the entire BTH region, could be represented by the first mode’s R² of 69.01%. Then, we compared the spatial pattern of the MCI intensity under different return periods. Using the May–June MCI of 1961–2023, the drought intensity gradually increased from northwest to southeast for the 10-, 20-, and 50-year return periods. However, the drought intensity decreased and then increased from northwest to southeast based on the 1991–2023 MCI. Notably, a 20-year return period of severe drought affected Beijing and northern Hebei in 2023. Finally, we discussed the linkages of drought in the BTH region and atmospheric circulation/sea surface temperature (SST) anomalies, which were calculated using the Pearson correlation coefficient and wavelet coherence. We suggest that the MCI variations in the BTH region may be related to the SST anomalies of the Indian Ocean in 1961–1990 and the Pacific Ocean in 1991–2023, respectively. The abovementioned studies have enlightened us to focus on predicting the Pacific SST for drought, which will facilitate agricultural production and water resource management in the BTH region.

1. Introduction

The frequency and intensity of extreme droughts are widespread and increasing with rapid global warming, significantly impacting ecosystems and society [1,2,3,4]. Extreme drought, characterized by a period of below-normal precipitation, is economically, ecologically, and socially disruptive, leading to food and water shortages [5,6]. To improve the ability to predict extreme droughts, it is valuable to understand the linkages between droughts and atmospheric circulation/sea surface temperature (SST) anomalies. Many studies have focused on the frequency and strength variations in regional droughts using instrumental meteorological data. Dai (2011) found a general global increase in drought, although with substantial regional variation and individual events dominating trend signatures in some regions [7], such as the extreme droughts that occurred in the Sahel during 1960–1990 [8] and continental extreme droughts in central Australia and the central USA from the 1950s to 1960s [9]. We compared the sc-PDSI calculated by van der Schrier et al. (2013) [10] and Dai (2011) [7], and found that the dominant mode of global drought variability is very similar, with a temporal evolution suggesting a trend towards drying. However, the other analysis for the 1950–2009 period shows an initial increase in drying in the Van der Schrier et al. [10] dataset, followed by a decrease from the mid-1980s onwards, while the Dai data show a continuing increase until 2000. Donat et al. (2013) [11] found that the annual maximum number of consecutive dry days has declined since the 1950s in more regions than it has increased. However, only regions in Russia and the USA indicate significant changes. Most studies on global dryness find decadal drought dominating longer-term trends [11,12,13]. Giorgi et al. (2011) [14] indicate that positive trends (reflecting an increase in the length of drought and/or extreme precipitation events) are most marked in Europe, India, and parts of South America and East Asia. However, trends appear to have decreased (reflecting a decrease in the length of drought and/or extreme precipitation events) in Australia and northern South America. Furthermore, as a result of global warming, intense drought occurrences will become more frequent [15], and the risk of megadroughts will rise as well [7]. Therefore, extreme droughts are hot topics in the climate science community.

Generally, drought is classified into four categories, including meteorological drought, hydrological drought, agricultural drought, and socioeconomic drought, by the American Meteorological Society [16]. Meteorological drought is a shortage of water due to an imbalance between precipitation and evapotranspiration payments and receipts during a certain period of time [17]. The widely used meteorological drought indices include the relative humidity index, standardized precipitation index (SPI) [18,19], standardized precipitation evapotranspiration index (SPEI) [20,21], Palmer drought severity index (PDSI) [22], and meteorological drought comprehensive index (MCI) [23,24]. In detail, the SPI, which is frequently used in drought monitoring, is estimated from the probability of precipitation occurring in a specific time period, but it ignores the effects of temperature and evapotranspiration [25]. The SPEI is constructed based on the SPI by introducing potential evapotranspiration. The monthly SPEI is sensitive to short-term precipitation and temperature changes, while the 6-month-scale SPEI better reflects subsoil moisture and stream runoff [26,27]. The MCI also considers the potential evapotranspiration and water shortages. These different drought indices have been widely used in drought research.

The Beijing–Tianjin–Hebei region (hereafter, the BTH region) is China’s “Capital Economic Circle” with a high population density and economic output, and it includes Beijing City, Tianjin City, and Hebei Province [28,29]. On the one hand, it is bordered by the Yanshan Mountains in the north and the Taihang Mountains in the west, with high vegetation cover. Precipitation in North China is characterized by significant seasonal and interdecadal variations [30,31]. From the 1950s to the 1970s, there was a dipole pattern of “floods in the north and droughts in the south”, with high precipitation in northern and northeastern China and low precipitation in the middle and lower reaches of the Yangtze River. Then, the rain belt moved southwards to southern China, and the precipitation in the eastern part of China often showed the triple-pole pattern of “+–+”. In the 1990s, the rain belt shifted northwards to the Yangtze River basin and changed to the dipole pattern of “floods in the south and droughts in the north” [32,33]. Summer precipitation in China underwent a significant climatic leap in approximately 1976, with North China and the Yellow River Basin experiencing decreased precipitation from 1977 [34,35,36]. Zhang et al. (2004) [37] analyzed the changes in precipitation in North China since 1880 and found that the precipitation in North China was low in the periods 1899–1920 and 1965–1997, in which the precipitation in 1980–1993 continued to be low and the drought phenomenon grew serious. Aerosol pollution caused by winter/spring forest wildfires, induced by extreme drought due to the East Asia winter monsoon system [33,38], may flow to urban agglomerations downstream and affect human health [39,40,41]. On the other hand, the central and southern parts of the BTH region belong to the North China Plain, which is largely made up of cropland and is one of China’s main grain-producing regions. Some studies have found that drought is the most important agrometeorological disaster in the North China Plain, with the largest annual disaster loss rate [42]. Between 1951 and 2010, there was a considerable increase in the intensity of drought disasters, and among all types of disasters, the amount of food lost due to drought and the risk of yield decline were the highest [43]. It is noted that North China experienced persistent hot waves and extreme drought in May–June 2023. Such extreme meteorological conditions have already caused accelerated moisture loss in cropland, hindering the growth of seedlings of some spring-sown crops, and affecting the progress of summer maize emergence in North China [44]. Therefore, we should intensify our research on the characteristics and causes of droughts in the BTH region.

This study aimed to explore the spatiotemporal variation in the MCI in the BTH region and drought intensity at different return periods. This study addressed the following questions: (1) What have the spatiotemporal characteristics of the MCI in the BTH region been over the last 60 years? (2) What were the drought intensities during the 10-, 20-, and 50-year return periods? The remainder of this paper is organized as follows: Section 2 introduces the materials and methods, Section 3 presents the results, Section 4 discusses the potential causes of drought during different time periods, and Section 5 presents the conclusions.

2. Materials and Methods

2.1. Study Area

The BTH region, with an area of 2.18 × 10⁷ hm², is located in the northern region of North China, including Beijing, Tianjin, and Hebei (Figure 1). In general, the BTH region experiences a warm, temperate, continental monsoon climate with cold, dry winters and hot, wet summers. Its average temperature ranges from less than 0 °C in January to approximately 18–27 °C in July (Figure 1c). The average annual precipitation is 484.5 mm, which decreases from southeast to northwest. This region has an average annual precipitation of 484.5 mm, which decreases from the southeast to the northwest (Figure 1d). The seasonal distribution is irregular, with precipitation in July and August prevailing since it is influenced by the East Asia summer monsoon [21,26,45]. The BTH region is distinguished by a high altitude in the northwest and a low altitude in the southeast. Forest predominates in the high-altitude mountainous regions of the north and west. Meanwhile, grassland is located in the western transitional area between forest and cropland. Urban land is also predominantly concentrated in and around large cities.

2.2. Data

The datasets used in this study mainly consist of the meteorological instrumental observation dataset from atmospheric circulation indices and SST indices.

The daily meteorological instrumental observation dataset from 1961 to 2023, including surface air temperature, precipitation, wind at 2 m, surface net radiation, etc., is provided by the National Meteorological Information Center (http://data.cma.cn/, accessed on 10 July 2023). There are 173 national meteorological stations in the BTH region. The dataset underwent rigorous quality control to ensure reliable data quality, such as eliminating outliers and data homogenization. These errors generally result from the relocation of meteorological station locations or replacement of meteorological instruments. On the one hand, the relocation of meteorological station locations generally results in an abrupt climate change in the time series of meteorological data. We first divided the time series into two segments using the abrupt point as a criterion. Then, we selected the longer series of the two series as the base period. Finally, we calibrated the mean and variance of the short series to a standard consistent with the long series using the idea of data homogenization. On the other hand, the sensitivity error of meteorological instruments generally results in outliers. We removed outliers when we found significantly large or small data by examining the time series. Meanwhile, some meteorological stations began recording data after 1961, so the period of meteorological observation data collection is insufficient in length. Therefore, we finally selected data from 85 national meteorological stations from 1961 to 2023 with no missing data and reliable quality in this study.

The land use dataset in 2015 was extracted from the National Land Use Datasets of China involved in the Resource and Environment Science and Data Center, Chinese Academy of Sciences (RESDC-CAS, https://www.resdc.cn/, accessed on 10 July 2023). These datasets are based on Landsat ETM/TM/OLI, ZY–3, HJ–1A, and other satellite images and rely on a unified standard artificial digital interpretation that comprises six classes (i.e., cropland, forest land, grassland, water, urban and rural construction land, and unused land) as well as 25 sub-classes. Importantly, the average classification accuracy of these datasets is over 90% [46]. Here, we primarily use six classes to discriminate different land types.

Atmospheric circulation and SST indices include 88 atmospheric circulations and 26 SST indices, which are derived from the National Climate Center, China Meteorological Administration [47]. The indices are calculated strictly in accordance with internationally recognized algorithms. Additionally, the data are updated monthly.

2.3. Methods

2.3.1. Meteorological Drought Comprehensive Index

The MCI is the result of the cumulative effect of the combined long- and short-term precipitation deficits [24]. The MCI considered the combined effects of effective precipitation over a 60-day period, evapotranspiration over a 30-day period, precipitation on a seasonal scale (90 days), and precipitation on a near-half-year scale (150 days). The index applies to the monitoring and assessment of daily meteorological drought. The smaller the MCI values, the greater the drought intensity, and vice versa. The calculation processes of the MCI are as follows:

$$\mathrm{M}\mathrm{C}\mathrm{I}=K_a\times \left(a\times SPI{W}_{60}+b\times {MI}_{30}+c\times {SPI}_{90}+d\times {SPI}_{150}\right)$$

(1)

where K_a is the seasonal adjustment coefficient. Generally, K_a is based on how sensitive crops are to soil moisture at growth and development phases through many field experiments. K_a varies with different months and regions. As summarized by Zhang (2017), the K_a values of the agroclimatic zone (including wheat area and corn area) in Beijing, Tianjin, and Hebei are consistent in each month, as shown in

Table 1. The values of the seasonal adjustment coefficient (K_a) for different months in the BTH region.

Month	1	2	3	4	5	6	7	8	9	10	11	12
Ka	0.4	0.8	1.0	1.2	1.2	1.2	1.2	1.0	1.0	0.8	0.6	0.4

. Furthermore, according to Zhang et al. (2017) [24], the impacts of the SPI and MI at different scales on the MCI varied with climate zone. Through many experimental and observational data, the weightings of SPI and MI at different scales are parameterized, thus obtaining a comprehensive index: the MCI. Through a large number of field experiments conducted by previous scholars, Zhang et al. (2017) [24] summarized that the weighting coefficients of a, b, c, and d of the weighting coefficients of SPIW₆₀, MI₃₀, SPI₉₀, SPI₁₅₀ in North China are constant, which are defined as 0.3, 0.5, 0.3, and 0.2, respectively.

SPIW₆₀ is the standardized weighted precipitation index for the last 60 days, which is defined by Equation (2).

$${SPIW}_{60}=SPI\left({\sum }_{n=0}^N{m}^n{P}_n\right)$$

(2)

where N is 1, 2, …, 60, and m is the contribution parameter in North China, defined as 0.85. P_n is the precipitation on the nth day.

MI₃₀ is the relative humidity index for the last 30 days, which is defined by Equation (3).

$${MI}_{30}=\frac{P-PET}{PET}$$

(3)

where P is the precipitation, and PET is the potential evapotranspiration, which is calculated by the functions of FAO Penman–Monteith, defined by Equation (4).

$$PET=\frac{0.408∆\left(R_n-G\right)+\gamma \frac{900}{T_{mean}+273}{u}_2(e_s-e_a)}{∆+\gamma (1+0.34u_2)}$$

(4)

where R_n is the surface net radiation, G is the soil heat flux, T_mean is the daily mean temperature, u₂ is the wind at 2 m, e_s is the saturated water vapor pressure, e_a is the actual water vapor pressure, Δ is the slope of saturated water vapor pressure curve, and γ is the dry and wet table constants.

SPI₉₀ and SPI₁₅₀ are the standardized precipitation indices for the last 90 days and 150 days, respectively, found in McKee et al. (1993) [48]. a, b, c, and d are the weight coefficients of SPIW₆₀, MI₃₀, SPI_90, and SPI₁₅₀, defined as 0.3, 0.5, 0.3, and 0.2 in the BTH region, respectively. According to the definition of Zhang et al. (2017) [24], the daily MCI below −0.5 and above −1 is considered as light drought, below −1 and above −1.5 as moderate drought, below −1.5 and above −2 as severe drought, and below −2 as extreme drought.

Based on the abovementioned calculation processes, we found that the MCI has been used as an indicator for drought monitoring in the National Climate Center, China Meteorological Administration, since it may capture monthly and seasonal precipitation anomalies and monthly scale water shortages [24]. The indicator is updated daily. Therefore, considering the applicability of indicators, the MCI was used in this study to estimate the drought characteristics in the BTH region.

2.3.2. Empirical Orthogonal Function on the MCI

The empirical orthogonal function (EOF) is usually used for single-variable time-varying characteristic decomposition [49]. It can decompose climate information into spatial and temporal dimensions (function). The spatial function characterizes the spatial distribution of the variable field, which does not change with time. The time function consists of a linear combination of spatial points. In this study, the MCI matrix is X_m×n, where m is the number of stations and n is the length of time. The covariance matrix C_m×n is defined as Equation (5).

$$C_{m\times n}=\frac{1}{n}(X_{m\times n}\times X_{m\times n}^T)$$

(5)

The eigenvalues (λ₁, λ₂, …, λ_m) and eigenvectors (V_m×m) of matrix C_m×n satisfy Equations (6) and (7).

$$C_{m\times n}\times V_{m\times n}=V_{m\times n}{\times E}_{m\times n}$$

(6)

$$E_{m\times n}=\left[\begin{array}{ccc}{\lambda }_1& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & {\lambda }_m\end{array}\right]$$

(7)

The time coefficient matrix T_m×n is defined by Equation (8).

$$T_{m\times n}=V_{m\times n}^T{\times X}_{m\times n}$$

(8)

To circumvent the possibility of random numbers in the actual data or any anomalies or false data mixed into the EOF analysis process, it is necessary to determine the eigenvalue of errors. The North test is used in this study [50].

$$∆\lambda =\lambda \sqrt{2/n}$$

(9)

where Δλ is the error range of the eigenvalue λ. The North significance test is passed when the error ranges of two neighboring eigen roots do not overlap via the sequential test and the error range labelling of λ, which means that the results of the empirical orthogonal function decomposition corresponding to the eigenvalues are physically meaningful.

2.3.3. Wavelet Coherence between the MCI and Atmospheric Circulation/SST Indices

Wavelet coherence is generally employed to examine the coherence of two time series in the time–frequency domain [51]. It is defined as follows:

$$R^2\left(s,t\right)=\frac{{\left|S\left[s^{-1}{W}_{xy}\right(s,t\left)\right]\right|}^2}{S[s^{-1}{\left|W_x(s,t)\right|}^2·s^{-1}{\left|W_y(s,t)\right|}^2]}$$

(10)

where s is the scale parameter, t is the time parameter, and S is the smoothing operator. W_x(s,t) and W_y(s,t) are the wavelet transforms of the time series x and y, respectively. W_xy(s,t) is the cross-wavelet transform. The wavelet coherence value R² is between 0 and 1, which is similar to the Pearson correlation coefficient. The values could be regarded as local correlation coefficients in time frequency space, with larger values indicating stronger correlations. When R² is close to 1, it can be assumed that the MCI and sea surface temperature index series are strongly correlated. When R² is close to 0, it can be assumed that the two series are basically irrelevant. The significance test algorithm is the Monte Carlo method [52].

To research the potential impact of atmospheric circulation/SST on the MCI in the BTH region, we calculated the wavelet coherence between the MCI and atmospheric circulation/SST indices. Wavelet coherence values passing the 95% significance test indicate that the MCI and atmospheric circulation/SST indices have the same period. The period variations in the MCI and atmospheric circulation/SST indices are synchronized and the two series are in phase when the arrows are right. The period variations in the MCI and atmospheric circulation/SST indices differ by half a cycle and the two series are in anti-phase when the arrows are left.

3. Results

3.1. Temporal Variations in MCI in the Beijing–Tianjin–Hebei Region

A boxplot of the monthly mean MCI in the BTH region during 1961–2023 is shown in Figure 2. The intraseasonal variation in the MCI followed a normal distribution with a mean value of −0.35. The median MCI rapidly decreased in spring at a rate of 0.10/month. The MCI for April to June were much lower than those for the other months throughout the year, ranging from −0.30 to −0.38. Afterwards, the MCI increased quickly, peaking at 0.01 in August. Notably, April to June (spring to early summer) was the driest time of the year, which was mainly due to the increase in temperature and evaporation during this period, and the drought reached its peak before the rainy season began. However, the monthly MCI for April had significantly more outliers which were approximately close to the median MCI values in other months (e.g., July, September, October). No significant outliers of the MCI were evident in May and June. Meanwhile, there was a significant (p < 0.1) difference between the medians of the MCI between April and May–June using the chi-square test. Whether there were outliers or medians, these suggested more consistent characteristics of the changes in the MCI for May and June, clearly different to the April MCI. Therefore, May–June was used as the target analysis period for a comprehensive meteorological drought in the Beijing–Tianjin–Hebei region.

3.2. Spatial Patterns of the MCI in the Beijing–Tianjin–Hebei Region

Figure 3 illustrates the spatial characteristics of the annual and May–June MCI during 1961–2022. Notably, only the differences between the annual and May–June MCI variations in 1961–2022 are compared because 2023 has not yet passed. The regional mean annual MCI for the BTH region was −0.11, ranging from −0.18 to −0.06. In the BTH region, the annual MCI pattern was described as “flood in the northeast, drought in the southwest”, where the drought intensity in the northeast was minimum and in the southwest was maximum with an MCI of approximately −0.2. The regional mean May–June MCI for the BTH region was −0.36, ranging from −0.43 to −0.22. The spatial distribution of the May–June MCI was drought throughout the whole BTH region. These results demonstrated that the trend of drought intensity increased from the northeast to the southwest in both the annual MCI and May–June MCI. In the BTH region, the standard deviation of the annual mean MCI and May–June MCI displayed a spatial pattern of “small in the north, large in the south”, which indicated that drought variability was low/high in the northern/southern portion. However, the drought variability of the May–June MCI was significantly higher than that of the annual MCI, suggesting that the BTH region experienced extreme droughts more frequently in May–June. Although the spatial pattern of “drought in the southeast, flood in the northwest” was the dominant feature, the annual MCI trend in the BTH region was insignificant. The drought intensity trend of the annual MCI was insignificant. The trend of the decreasing drought intensity of the May–June MCI reached ~1/10a. Meanwhile, the most significant changes in the May–June MCI were in the northern mountainous areas and southern cultivated areas. This indicated that the wetter trend of the BTH region was more significant in May–June than the annual trend.

The abovementioned findings suggested that the MCI changes in May–June were more obvious than the annual changes. The mean MCI was mainly dominated by drought, with a trend of decreasing drought intensity. Furthermore, there was a higher probability of extreme droughts in the southern region, where there was greater drought variability.

We used the data in the period of 1961–2023 to analyze the characteristics of the May–June MCI in the following study. Figure 4 depicts the top two dominant spatial patterns of the EOF mode for the May–June MCI in the BTH. We found that the R² of EOF1 was 69.01%, and its spatial pattern was characterized by largely consistent MCI changes in the BTH region. The R² of EOF2 was 10.23%, and its spatial pattern was a “northeast–southwest” dipole. The cumulative variance explained by the first two modes reached 79.24%, which indicated that the first two modes could essentially reflect the main spatial pattern of the May–June MCI in the BTH region. The R² of the first mode was close to 70%, indicating that it could represent the predominant spatial pattern of the BTH, i.e., the consistent change in the MCI throughout the entire region. Hence, using the results of EOF1, we analyzed the differences in the probability density function pattern of the MCI series between different time periods (i.e., 1961–1990 and 1991–2020) and the wavelet coherence between the MCI series of the whole Beijing–Tianjin–Hebei region and atmospheric circulation/SST indices series to research the potential impact of atmospheric circulation/SST on the MCI in the BTH region.

3.3. Comparison of the MCI Characteristics between 1961–1990 and 1991–2023

Figure 5 illustrates the probability density function pattern of the May–June MCI over different time periods. The May–June MCI was found to have a normal distribution in the time periods of 1961–2023, 1961–1990, and 1991–2023. The mean value of the MCI was approximately −0.35, which indicated that the BTH region had primarily experienced drought over the past 60 years. The probability of extreme events was highest during 1961–1990, followed by 1961 to 2023, and then 1991 to 2023, while the probability of average climatology was reversed. These results showed that the meteorological drought conditions in the BTH region over the past 30 years had gradually changed from extreme to common droughts, and that the extremeness had been reduced.

3.4. Different Return Periods of the MCI in the Beijing–Tianjin–Hebei Region

It is reported by the Ministry of Emergency Management that North China experienced several persistent heat waves and extreme drought events from May to June in 2023, which have resulted in decreasing crop yields and increasing forest fire danger [20]. To analyze the historical status of 2023 persistent extreme drought events in the BTH region over the past several decades, we calculated the BTH return periods based on different time periods and compared them with the May–June MCI intensity in 2023. Figure 6 shows the spatial distribution of the May–June MCI intensity during the 10-, 20-, and 50-year return periods and the May–June MCI intensity in 2023. We found that the drought intensity of the MCI gradually increased from northwest to southeast based on the calculated 10-year return period from 1961 to 2023. The drought intensity of the MCI in 2023 in Beijing and northern Hebei exceeded the drought intensity of the 10-year return period. The drought intensity of MCI first decreased and then increased from northwest to southeast based on the 10-year return period from 1991 to 2023.The drought intensity of the MCI in 2023 in Beijing and northern Hebei was also higher than that of the 10-year return period. The drought intensity of the MCI gradually increased from northwest to southeast based on the 20-year return period calculated from 1961 to 2023. The drought intensity of the MCI in 2023 in southern Beijing, northwestern Beijing, and northern Hebei exceeded the drought intensity of the 20-year return period. The drought intensity of the MCI first decreased and then increased from northwest to southeast based on the 20-year return period from 1991 to 2023. The drought intensity of the MCI in 2023 in Beijing and northern Hebei was also higher than that of the 20-year return period. The extent of 2023 exceeding the 20-year return period calculated from 1991–2023 was greater when compared to 1961–2023. The drought intensity of the MCI gradually increased from northwest to southeast based on the 50-year return period calculated from 1961 to 2023. The drought intensity of the MCI in 2023 only in the northern portions of Hebei exceeded the drought intensity of the 50-year return period. The drought intensity of the MCI first decreased and then increased from northwest to southeast based on the 50-year return period from 1991 to 2023. The drought intensity of the MCI in 2023 only in the northern portions of Hebei was higher than that of the 50-year return period.

The abovementioned findings suggested that the spatial pattern of the MCI drought intensity varied with different time periods, although the same return period was calculated. The spatial pattern of the MCI intensity in the BTH region increased gradually from northwest to southeast for the 10-, 20-, and 50-year return periods calculated based on 1961–2023. However, the spatial pattern of the MCI intensity first decreased and then increased from northwest to southeast for the 10-year, 20-year, and 50-year return periods calculated based on 1991–2023. Whether based on 1961–2023 or 1991–2023, the drought intensity of the MCI in 2023 in Beijing and northern Hebei reached the 20-year return period. However, drought extents greater than the 20-year return period were more extensive based on the period of 1991–2023. Meanwhile, the drought intensity in 2023 reached a 50-year return period only in some parts of northern Hebei. These results suggested that the 2023 drought in northern China (NC) was likely to spread from northern NC to southern NC, which may be due to the influence of a high-pressure system in Northeast China.

4. Discussion

To discuss the linkage of the MCI in the BTH region and atmospheric circulation/SST anomalies, we calculated the correlation coefficients between the regional mean May–June MCI series of the BTH region and a total of 114 indices, including 88 atmospheric circulation indices series and 26 SST indices series.

Table 2. Top five correlation coefficients of the May–June MCI and atmospheric circulation, SST indices.

Atmospheric Circulation/SST Indices	Correlation Coefficient
Pacific Polar Vortex Area Index	−0.392 **
NINO B SSTA Index	0.368 *
Indian Ocean Warm Pool Area Index	0.389 **
Indian Ocean Basin-Wide Index	0.390 **
Indian Ocean Warm Pool Strength Index	0.398 **

Note: ** and * passed significance at a confidence level of 0.01 and 0.05, respectively.

shows the top five correlation coefficients of the May–June MCI series and atmospheric circulation indices/SST indices series. We found that the top five indices with the largest correlation coefficients with the MCI variations were, respectively, the Pacific Polar Vortex Area Index (PPVAI), Indian Ocean Warm Pool Area Index (IOWPAI), Indian Ocean Basin-Wide Index (IOBWI), Indian Ocean Warm Pool Strength Index (IOWPSI), and NINO B SSTA Index. In particular, the PPVAI is defined as the area of the sector enclosed by the polar vortex north of the characteristic contour of the southern boundary of the vortex in the Northern Hemisphere in the 500 hPa height field, 150° E–120° W region. The IOWPAI is defined as a spherical area of areas where the sea surface temperature exceeds 28 °C in the regions 7° S–30° N, 41° E–98° E and 30° S–8° N, 41° E–120° E. The IOBWI is defined as the regional average of the sea surface temperature distance levels in the region 20° S–20° N, 40° E–110° E. The IOWPSI is defined as a spherical area of areas where the sea surface temperature exceeds 28 °C in the regions 7° S–30° N, 41° E–98° E and 30° S–8° N, 41° E–120° E. The NINO B SSTA Index is defined as the regional average sea surface temperature anomaly in the region 0°–10° N, 50° E–90° E. The MCI had a negative correlation with the PPVAI, which had a correlation coefficient of −0.392 (p < 0.01). However, the MCI had positive correlations with the other four indices, which had correlation values of 0.368 (p < 0.05), 0.389 (p < 0.01), 0.390 (p < 0.01), and 0.398 (p < 0.01). Since the correlation coefficients between the MCI and the four indices (i.e., PPVAI, IOWPAI, IOBWWI, and IOWPSI) were the largest, passing the 99% significance test, and significantly higher than that with the NINO B SSTA Index, the wavelet coherences between the May–June MCI series and SST indices, except the NINO B SSTA index, were examined (Figure 7). We discovered that the MCI and PPVAI had significant quasi-4–6-year oscillations since the 1990s–2010s. Additionally, the anti-phase relationship between the MCI and PPVAI was confirmed by the negative correlation between the two indices. During the 1980s, there were significant quasi-decadal oscillations between the MCI and PPVAI, and the two indices exhibited an anti-phase variation relationship. However, the MCI and IOWPAI, IOBWI, and IOWPSI all had quasi-2–4-year oscillations between 1960 and 1990, and there were phase relationships between the MCI and the three SST indices. Moreover, there was quasi-decadal oscillation between the MCI and the IOWPAI and IOWPSI. Numerous studies have examined how ocean systems affect drought/flood variations in North China. The Pacific SST in the mid-high latitude, especially the Pacific Decadal Oscillation (PDO) [53] and the Indian SST [54], has a good correspondence with precipitation in North China. For instance, when the PDO is in the warm phase, summer precipitation in North China is anomalously low, and when the PDO is in the cold phase, summer precipitation in the region is anomalously high, and the PDO mainly influences summer precipitation in North China through the East Asian summer monsoon [55]. Hao et al. (2012) [54] showed that the consistently high sea surface temperature in the tropical Indian Ocean leads to more precipitation and weaker winter winds in East Asia, which further causes weaker thermal contrasts between the land and sea in the spring and summer seasons and weak southerly winds in the eastern part of China and, correspondingly, less precipitation during the flood season in the northern part of the country.

These findings suggested that there was a quasi-2–4-year oscillation between the MCI, IOWPAI, IOBWI, and IOWPSI from 1961 to 1990. The three SST indices over the Indian Ocean may be related to the higher probability of extreme droughts in the BTH region. There was a quasi-4–6-year oscillation between the MCI and IOWPAI over the Pacific from 1991 to 2023, which may have declined the probability of extreme droughts in the BTH region. However, there are some weaknesses or limitations of our study. On the one hand, we calculated the wavelet coherence between the MCI and atmospheric circulation/SST indices to discuss the linkage of drought in the BTH region and atmospheric circulation/SST anomalies. However, this is just a statistical relationship. Studying how SST anomalies affect atmospheric circulation and alter local climate change, thereby triggering local droughts in North China, is equally important. In future research, we may consider analyzing the relationship and potential mechanisms between the three oceans and drought/flood in North China, which could help us to better predict the drought/flood in North China. On the other hand, the main focus has been on the climate states of monthly MCI in this study. In the future, more attention must be paid to the characteristics of continuous extreme drought, such as the number of consecutive MCI days.

5. Conclusions

The abovementioned results demonstrated that in the BTH region, the meteorological drought intensity from May to June was the most severe in the intraseasonal variation. The mean MCI in the BTH region was mainly dominated by drought. The drought intensity trend of the annual MCI was insignificant. The drought intensity trend of the MCI in May–June generally decreased throughout the whole BTH region. Meanwhile, the most significant changes of the May–June MCI were in the northern mountainous areas and southern Plain. In addition, there was a higher probability of extreme droughts in the southern area of the BTH region, where there was greater drought variability. The MCI changes in May–June were more obvious than those in the year. The cumulative variance explained by the first two modes reached 79.24%. The primary spatial pattern of MCI was characterized by consistent MCI changes throughout the whole BTH region, with the R2 reaching 69.01%. Compared with 1961–1990, meteorological drought conditions in the BTH region over the past 30 years gradually changed from extreme to common droughts. The spatial pattern of the MCI intensity in the BTH region increased gradually from northwest to southeast for the 10-, 20-, and 50-year return periods calculated based on 1961–2023. However, the spatial pattern of the MCI intensity was more pronounced in the northwest and softer in the southeast based on the 1991–2023 MCI. A 20-year return period of severe drought affected Beijing and northern Hebei in 2023. The drought intensity of the MCI in 2023 in northern Hebei alone reached the 50-year return period. These results suggested that the 2023 drought in NC was likely to spread from northern to southern NC. Moreover, the droughts in the BTH region may be related to IOWPAI, IOBWI, and IOWPSI, with 2–4 years of oscillation over the Indian Ocean, and IOWPAI, with 4–6 years of oscillation over the Pacific during 1961–1990 and 1991–2023, respectively. Moreover, our study mainly relied on an observational meteorology dataset. In the future, we need to research future projections so that government officials can precisely identify extreme drought events and put effective preventive measures in place to effectively improve local agricultural production and water resource management.