Identiﬁcation and Spatiotemporal Migration Analysis of Groundwater Drought Events in the North China Plain

: Groundwater droughts can explain developments and changes in groundwater from a climatological perspective. The North China Plain (NCP) is a typical underground funnel area. Therefore, groundwater drought studies in the NCP can provide better understanding of the local hydrogeological characteristics from new perspectives. In this paper, the GRACE groundwater drought index (GGDI) was used to evaluate groundwater drought events in the NCP. Additionally, a new method was proposed in this study for investigating groundwater drought events at the spatiotemporal scale. On this basis, the centroid theory was used to construct an appropriate groundwater drought migration model for the NCP. The results showed that (1) the groundwater drought frequency in the NCP was 24.54%. In addition, the most severe groundwater drought events in the study occurred in March 2020. (2) In total, 49 groundwater drought events occurred in the NCP over the 2003–2020 period. The most intense groundwater drought event occurred over the June 2018– December 2020 period (DE.49), covering the entire study area. DE.29 was the second most intense groundwater drought event over the August 2012–September 2013 period (14 months), resulting in a maximum arid area of 75.57% of the entire study area. (3) The migration of the groundwater drought events was in the southwest–northeast and northeast–southwest directions, which was consistent with the terrain inclination, while most of the groundwater drought centroids were concentrated in Area II. The groundwater drought event identiﬁcation method and the groundwater drought migration model were effective and reliable for assessing groundwater drought events in the NCP and provided a better understanding of developments and changes in groundwater droughts, which is of great practical signiﬁcance and theoretical value for the rational development and use of groundwater resources, as well as for guiding industrial and agricultural activities.


Introduction
Groundwater is an important water source for human beings and has accompanied the development of human civilization. In arid and semi-arid regions, the development of human society depends on depleted groundwater resources [1]. However, the development of human society has disrupted the natural cycle of groundwater, causing the depletion of aquifers around the world at an alarming rate [2][3][4]. The UN World Water Development Report 2022 focused on groundwater issues, with the conference theme of "Groundwater: making the invisible visible" [5], the same theme as that of World Water Day 2022, highlighting groundwater problems. Climate change is a new challenge facing human civilization. The fifth assessment report of the Intergovernmental Panel on Climate Change (IPCC) highlighted a potential increase in the frequency of extreme climate events [6]. In addition, several studies have demonstrated significant increases in the temperatures of arid regions worldwide over the past 100 years, resulting in drought areas almost doubling, which may potentially continue to expand and cause the further intensification of drought degree [7][8][9]. Measures to combat climate warming are required urgently [10,11].
There is still no clear and unified definition of the groundwater drought concept. Indeed, scholars have interpreted groundwater droughts differently according to different assessment methods. The "groundwater drought" term was first proposed by Rutulis [12] to distinguish between droughts that occur in groundwater aquifers and other drought types (e.g., meteorological, hydrological, and agricultural droughts). He defined a groundwater drought as a significant natural drop in groundwater levels, causing substantial aquifer dewatering, thus resulting in severe water supply problems and emphasizing natural conditions, non-human interference, and the transmission of meteorological droughts to groundwater droughts. Groundwater droughts under natural conditions are often cyclical and caused by reduced groundwater recharge and storage over periods of time, while other drought types are due to insufficient precipitation [13]. However, groundwater droughts can also be caused by anthropogenic factors, such as groundwater extraction rates that are greater than recharge and groundwater storage amounts, which thus intensifies groundwater droughts to a certain extent. In either case, groundwater droughts are characterized by low groundwater levels and low groundwater abstraction rates (or even dry wells). Therefore, groundwater droughts may be defined as a lack of groundwater recharge or the lack of groundwater in terms of reserves or groundwater heads over a specific period of time [13]. However, groundwater droughts cannot only be evaluated from the perspectives of natural and human factors. In fact, Marchant and Bloomfield [14] defined a groundwater drought as a decrease in groundwater level below the normal level or a reduction in spring flow. Van Lanen and Peter [13] believed that groundwater droughts occur in three consecutive phases: temporal decreases in groundwater recharge, followed by decreases in groundwater levels and decreases in discharge rates. Indeed, the emergence of groundwater droughts can first affect groundwater ecosystems and then spread to other ecosystems. In addition, a groundwater drought is a type of drought that occurs following a meteorological drought that has caused hydrological and agricultural drought events [6]. Drought propagates into groundwater systems through mechanisms such as pooling in river catchments, hysteresis, and the prolongation of drought signals [15]. In the absence of long-term precipitation deficiency, variations in groundwater drought characteristics could be affected by anthropogenic warming activities [16]. Therefore, the occurrence of groundwater droughts is an indicator of the significant impacts of natural and/or anthropogenic factors on the water cycle (e.g., severe climatic influences and groundwater overuse).
At present, the most commonly used methods for evaluating groundwater droughts are the standardized precipitation index (SPI) [17], the standardized groundwater level index (SGI) [18,19], the standardized water level index (SWI) [20], and the GRACE groundwater drought index (GGDI) [21,22]. In addition, other methods have also been used in some studies to assess groundwater droughts, such as the groundwater drought index (GWI) [23], the groundwater resource index (GRI) [24], and the groundwater recharge drought index (GRDI) [25]. Among them, the GGDI has been extensively used to assess groundwater droughts because it is based on remote sensing data, which are highly available, thereby overcoming the limitations of spatiotemporal groundwater drought monitoring [26,27].
Some previous studies have assessed the temporal and spatial evolution of drought events separately, resulting in inaccurate descriptions of the evolution process of drought events [28]. Therefore, the evolution of drought dynamics at the spatiotemporal scale has attracted considerable attention from some scholars. Zhou et al. [29] proposed a new method for studying meteorological drought migration based on topological spatial relationships to accurately describe the spatiotemporal dynamic behavior of droughts. Guo et al. [30] used an improved three-dimensional clustering algorithm (longitude-latitude-time and space) to identify drought episodes and described the characteristics of drought events according to several indicators, including drought duration, severity, intensity, area, centroid, and trajectory. Herrera-Estrada et al. [31] used Lagrangian methods to monitor worldwide spatiotemporal drought events over the 1979-2009 period, analyzing their characteristics and behaviors. Wen et al. [28] proposed a new 3D drought structure extraction method to assess the long-term spatiotemporal distribution patterns of droughts during their development. Han et al. [32] constructed an agricultural drought migration model using the migration of drought centroids to describe the spatiotemporal evolution of agricultural droughts. Understanding the spatiotemporal variations in groundwater droughts is of great significance for the comprehensive assessment of groundwater droughts. However, other studies on groundwater drought events have been carried out from the perspective of time [12,33,34]. Therefore, it is necessary to assess groundwater drought events at both the temporal and spatial scales.
The North China Plain (NCP) is one of three plains in China. Indeed, it has attracted much attention because of its relevant political and economic characteristics, its importance in food production, and its severe water shortage. The large groundwater funnel in this area makes it an interesting area for hydrogeological studies, thus contributing to sustainable groundwater development [35]. However, in recent years, relevant departments have implemented certain measures in the NCP, such as groundwater overexploitation controls, the South-to-North Water Diversion, and ecological water replenishment, to prevent the continuous decrease in the groundwater level in the NCP, especially in urban areas [36]. However, although these measures have resulted in slight increases in groundwater levels, the magnitude of the groundwater level drop created by long-term groundwater depletion is still huge [37]. In the context of climate change, groundwater drought research in the NCP needs to consider climate change aspects and discuss groundwater changes from a climatological perspective. Although groundwater management policies have been effective in slightly increasing groundwater levels in the NCP, the effects of climate warming on groundwater droughts are still unclear. In addition, the spatiotemporal variations in groundwater drought events in the NCP need to be assessed. Indeed, Wang et al. [34] used GGDI to evaluate groundwater drought events in the NCP over the 2003-2015 period, identifying the time of the most severe groundwater drought events in the region and analyzing the influence of teleconnection factors on groundwater drought events. The results demonstrated the significant impact of ENSO on groundwater drought events in the NCP. Meanwhile, Zhang et al. [38] studied groundwater drought trends in the NCP using the GWI and singular spectrum analysis. However, they assessed the spatial heterogeneity of groundwater drought events in a fixed time without analyzing temporal variations in the groundwater drought events. Therefore, further studies on the spatiotemporal changes in groundwater drought events in the NCP may provide a better understanding of the evolution characteristics of groundwater droughts, which is crucial for ensuring the sustainable use of groundwater resources.
The main purpose of this study was (1) to quantitatively assess groundwater drought events in the NCP using the GGDI, (2) propose a new method for investigating the spatiotemporal patterns of groundwater drought events in the NCP from 2003 to 2020, and (3) construct a specific groundwater drought migration model to analyze the characteristics of groundwater drought migration in the NCP. Therefore, this study interpreted developments and changes in groundwater from the perspective of climatology and examined the spatial and temporal evolution of groundwater drought events from a new perspective, which is of great practical significance and theoretical value for the rational exploitation of groundwater resources and the guidance of industrial and agricultural production.

Study Area
The NCP (113 • 10 -119 • 25 E; 34 • 52 -40 • 29 N) is one of the three major plains in China ( Figure 1). It is located in the southern part of the Taihang Mountains and the northern part of the Yellow River, covering a total area of about 13.9 × 104 km 2 [39]. The geomorphological features of the NCP can be divided into three main areas, namely, the piedmont alluvial-pluvial inclined plain (Area I), the central and alluvial lacustrine plain (Area II), and the eastern alluvial marine plain (Area III). According to the type of groundwater, the pre-mountain plain in North China can be divided into carbonate rock karst water, bedrock fracture water, and loose rock pore water. The available groundwater is mainly stored in loose pore space, and the regional flow direction is northwest to southeast. Locally, due to the formation of landing funnels by groundwater overexploitation, the groundwater dynamic field has changed and groundwater no longer flows in one direction, forming a seepage field from the edge of the landing funnels to the center of each funnel [40]. Groundwater in the NCP circulates mainly in the quaternary pores of underground aquifer rock series. According to burial conditions, circulation characteristics, and retention time, groundwater in the NCP can be divided into two types, namely, shallow and deep groundwater. The shallow groundwater level varies from 50 to 210 m, while the depth of the deep groundwater level ranges from 100 m in Area I to 600 m in Area III. On the other hand, the aquifer lithology in the study area consisted of gravel, medium-coarse sand, and mediumfine sand in Area I, medium-fine sand, fine sand, and silty sand in Area II, and silty fine sand and silty sand in Area III [41]. The NCP has a semi-arid and semi-humid continental monsoon climate, with an average annual rainfall of about 500-600 mm, occurring mainly between July and September [42]. In addition, the evaporation rates increase with air temperature and decrease with latitude. Indeed, the uneven spatiotemporal distributions of precipitation and evaporation have directly impacted groundwater resources in the NCP [37].
part of the Yellow River, covering a total area of about 13.9 × 104 km 2 [39]. The geomorphological features of the NCP can be divided into three main areas, namely, the piedmont alluvial-pluvial inclined plain (Area I), the central and alluvial lacustrine plain (Area II), and the eastern alluvial marine plain (Area III). According to the type of groundwater, the pre-mountain plain in North China can be divided into carbonate rock karst water, bedrock fracture water, and loose rock pore water. The available groundwater is mainly stored in loose pore space, and the regional flow direction is northwest to southeast. Locally, due to the formation of landing funnels by groundwater overexploitation, the groundwater dynamic field has changed and groundwater no longer flows in one direction, forming a seepage field from the edge of the landing funnels to the center of each funnel [40]. Groundwater in the NCP circulates mainly in the quaternary pores of underground aquifer rock series. According to burial conditions, circulation characteristics, and retention time, groundwater in the NCP can be divided into two types, namely, shallow and deep groundwater. The shallow groundwater level varies from 50 to 210 m, while the depth of the deep groundwater level ranges from 100 m in Area I to 600 m in Area III. On the other hand, the aquifer lithology in the study area consisted of gravel, medium-coarse sand, and medium-fine sand in Area I, medium-fine sand, fine sand, and silty sand in Area II, and silty fine sand and silty sand in Area III [41]. The NCP has a semi-arid and semihumid continental monsoon climate, with an average annual rainfall of about 500-600 mm, occurring mainly between July and September [42]. In addition, the evaporation rates increase with air temperature and decrease with latitude. Indeed, the uneven spatiotemporal distributions of precipitation and evaporation have directly impacted groundwater resources in the NCP [37]. The exploitation of groundwater in the North China Plain began in the mid-20th century, and the addition of human activities has led to significant variability in groundwater systems. Humans have influenced the evolution of groundwater systems through ground- The exploitation of groundwater in the North China Plain began in the mid-20th century, and the addition of human activities has led to significant variability in groundwater systems. Humans have influenced the evolution of groundwater systems through groundwater extraction, river management, sewage discharge and management, and the construction of water conservancy. In recent decades, groundwater extraction in the North China Plain has intensified, and the form of the groundwater cycle has been greatly altered. According to statistics, until the 1960s, the smaller scale of groundwater extraction left the groundwater cycle in its natural state, with precipitation as its main source of recharge and evaporation as its main mode of discharge. After the turn of the century, precipitation remained the main factor in maintaining the dynamic balance of groundwater, but groundwater was subjected to higher intensity extraction, resulting in the thickening of the bale zone in the pre-mountain plain area, an increase in the depth of groundwater burial, and a serious lag in the infiltration of precipitation for recharge. In parts of the eastern plain, although the infiltration of precipitation for recharge has increased, changes in the substrata due to urbanization have hindered the infiltration of precipitation and the main mode of groundwater discharge has changed from evaporation to artificial extraction. Deep groundwater in the North China Plain is buried deeper, circulation is slower, and it is a non-renewable resource. In its natural state, the recharge of the deep groundwater mainly comes from runoff and vertical infiltration from alluvial fans and hidden carbonate karst water in front of the mountains.

Gravity Recovery and Climate Experiment (GRACE)
According to their different inversion principles, the GRACE gravity satellite inversion methods can be divided into two categories, namely, the spherical harmonic coefficient and Mascon methods. Terrestrial water storage change data, provided by the latest RL06Mascon data product released by the University of Texas Space Research Center (http://www2.csr. utexas.edu/grace/RL06_mascons.html, accessed 28 October 2021 ), were used in this study to invert terrestrial water storage anomalies. The data covered the April 2002-September 2021 period, with a monthly scale and a spatial resolution of 0.25 • × 0.25 • . Some missing data were handled using an efficient and widely used linear interpolation method [29,43]. The data from January 2003 to December 2020 were selected as the basic data after inverting the groundwater storage change data to calculate the GWI.

Global Land Data Assimilation System (GLDAS)
The global land data assimilation system (GLDAS), jointly developed by the National Aeronautics and Space Administration and the National Oceanic and Atmospheric Administration, can use data assimilation technology to merge satellite data and ground observation data and generate surface state quantities and flux [44]. In total, four data products were launched based on the four land surface models, namely, VIC, CLM, CLSM, and Noah. Indeed, the Noah10_M 2.1 (https://disc.sci.gsfc.nasa.gov/data-holdings, accessed on 19 November 2021.) data product provides consistent data with the same spatial resolution as the GRACE RL06Mascon. Therefore, the Noah10_M 2.1 data product, covering the April 2003-September 2021 period, was used to invert the GRACE groundwater storage change data.

GGDI
The GRACE satellite data were first used to invert the terrestrial water storage anomalies (TWSAs) and then the GLDAS hydrological model was used to invert the surface water storage anomalies (SWSAs), soil moisture storage anomalies (SMSAs), snow water equivalent anomalies (SWEAs), and canopy water storage anomalies (CWSAs). Afterward, the GGDI was calculated according to Thomas et al. [21], as follows: where Ci is the mean GWSA of the ith month (i = 1, 2, 3, ..., 12) and n is the number of the ith month in the data column.
To remove any seasonal effects on the final results, Ci was deducted from each month's GWSA using the following formula: where GSD is the net groundwater storage deviation.
In addition, to reflect drought conditions, the normalized net groundwater storage deviation (GGDI) was calculated using the following formula: where x GSD is the mean value of GSD and S GSD is the standard deviation of GSD.

Groundwater Drought Event Identification
In previous studies, most groundwater drought events have been determined from a temporal perspective. This study proposed a new method for identifying groundwater drought events from spatial and temporal perspectives, using the following steps: (1) Drought grid identification The study area was divided into 221 grids, with grid cells of 0.25 × 0.25. When the center of gravity of the grid was within the study area, the grid was considered to belong to the study area ( Figure 2a). When the GGDI of the grid was less than −0.5, the grid was considered to belong to the drought area [29].
Atmosphere 2023, 14, x FOR PEER REVIEW 7 of 19 Figure 2. The groundwater drought event identification process (a) shows the drought grid identification, (b) shows the drought event spatial continuity identification, (c,d) shows the drought event temporal continuity identification, and (e-h) shows the drought event spatial-temporal continuity identification.)

Construction of the Drought Migration Model (DMM)
A centroid refers to a hypothetical point where mass is considered to be concentrated within a material system. Centroids were considered in this study to analyze the spatiotemporal migration characteristics of groundwater drought events. The specific steps were as follows: (1) Determination of groundwater drought events Previous studies have shown that only drought events that last for longer than 3 months can reflect the characteristics of drought migration [32]. Therefore, only drought events with a duration greater than or equal to 3 months were considered for the construc- (2) Identification of the spatial continuity of drought events Overlapping drought grid points were observed in certain months, which showed the same drought event (a, c, and d in Figure 2b). Therefore, a drought event was considered absent when the drought area was smaller than the minimum drought area (A 0 ) (b and e in Figure 2b) [45]. Indeed, A 0 was equal to 1.6% of the total surface of the study area [46], representing four grids in this study.
(3) Identification of the temporal continuity of drought events After identifying the arid areas, the temporal continuity of the drought events was determined. Indeed, a temporal drought event was considered continuous, i.e., belonging to the same drought event, when the overlapping area between the arid areas of two adjacent months was greater than or equal to A 0 (Figure 2c,d).
(4) Identification of the spatiotemporal continuity of drought events The drought events m (DE-m) and n (DE-n) in month t were merged into a single drought event in month t + 1. The merged event was called DE-n since the DE-m grid in the merged drought event was smaller than that of the DE-n grid. The original drought events m and n were two sub-events of the merged drought event n, denoted by (n-1) and (n-2), respectively. (Figure 2e,f).
As shown in Figure 2g, DE-(n-1) and DE-(n-2) belonged to DE-n, even though DE-n was divided into several spatially disconnected drought events in the next month ( Figure 2h) and the overlapping DE-(n-1), DE-(n-2), and DE-n grids were not less than A 0 .

Construction of the Drought Migration Model (DMM)
A centroid refers to a hypothetical point where mass is considered to be concentrated within a material system. Centroids were considered in this study to analyze the spatiotemporal migration characteristics of groundwater drought events. The specific steps were as follows: (1) Determination of groundwater drought events Previous studies have shown that only drought events that last for longer than 3 months can reflect the characteristics of drought migration [32]. Therefore, only drought events with a duration greater than or equal to 3 months were considered for the construction of the DMM.
(2) Determination of the locations of groundwater drought centroids The coordinates of the groundwater drought centroids were determined using the following equations [32]: where X and Y are the longitude and latitude of the groundwater drought centroid, respectively, i represents the ith groundwater drought grid, n denotes the total number of grids contained in the groundwater drought area, X i and Y i are the longitude and latitude of the ith groundwater drought grid, respectively, and P i denotes the GGDI of the ith groundwater drought grid.
(3) Connecting the groundwater drought centroids The trajectory, direction, speed, and other characteristics of groundwater drought migration were determined by connecting the groundwater drought centroids in chronological order.
(4) Calculation of the groundwater drought intensities Centroids can not only represent the spatiotemporal characteristics of points but they can also be used as indicators to measure monthly drought intensity. According to previous studies, the centroid of the GGDI (S am ), representing a drought intensity of DE. (S a ) in month m, can be calculated using the following equations [32]: where k denotes the drought grid number of DE.a in month m and S 0 denotes the groundwater drought threshold (−0.5), and where n denotes the total number of groundwater drought months in DE.a.

GRACE Data Validation
The reliability and authenticity of the GLDAS-based SWSA, SMSA, SWEA, and CWSA data have been confirmed in previous studies; thus, they could be combined with the GRACE data to invert the groundwater storage data [47][48][49]. The GLDAS data were used in this study to validate the GRACE data. The GRACE-and GLDAS-based storage changes showed significant positive correlation coefficients (p < 0.01) in each geomorphic area of the NCP (p < 0.01). The obtained Pearson correlation coefficients ranged from 0.322 to 0.552 (Table 1). The GLDAS-and GRACE-based water storage in the NCP ranged from 13.06 to 6.81 mm and 13.76 to 41.18 mm, respectively, indicating relatively large values. These results could be explained by the fact that the GRACE-based land water storage included both surface water and groundwater, while that of the GLDAS excluded groundwater storage. The GRACE uncertainties in the NCP and its sub-areas were calculated according to the method proposed by Landerer and Swenson [50] ( Table 2). The results showed decreasing trends in the GRACE-based water storage in the NCP and its sub-areas (Table 3). According to the specific situation of the NCP, these decreasing trends could have been due to the ecological restoration projects implemented by the government increasing the consumption of soil water in the study area, as well as groundwater overexploitation. The temporal trends in the GRACE-and GLDAS-water storage in sub-areas of the NCP are shown in Figure 3, indicating consistent changes in water storage from 2003 to 2013. However, gradual increases in the differences between the GRACE-and GLDAS-water storage trends were observed after 2013. Indeed, Area II exhibited the largest mean difference (−19.36), followed by Area I and Area III, with mean differences of −14.80 and −9.30, respectively. These findings suggested increased rates of decline in groundwater storage in the NCP after 2013. According to the obtained results, Area II exhibited the highest decline rate, followed by Area I and Area III. This finding was consistent with those reported by Zhang et al. [38]. In summary, the obtained GRACE data for the NCP were reliable and could be used to effectively investigate groundwater drought events in this study area.  trends were observed after 2013. Indeed, Area II exhibited the largest mean difference (−19.36), followed by Area I and Area III, with mean differences of −14.80 and −9.30, respectively. These findings suggested increased rates of decline in groundwater storage in the NCP after 2013. According to the obtained results, Area II exhibited the highest decline rate, followed by Area I and Area III. This finding was consistent with those reported by Zhang et al. [38]. In summary, the obtained GRACE data for the NCP were reliable and could be used to effectively investigate groundwater drought events in this study area.

Quantitative Assessment of Groundwater Drought Events
The obtained quantitative evaluation results of groundwater drought events in the NCP and its geomorphic sub-areas using the GGDI are shown in Figure 4. The groundwater drought trends in the geomorphic sub-areas were consistent with those observed across the entire NCP, showing gradually decreasing trends. In addition, the most severe groundwater drought event in all areas was observed in March 2020. Indeed, the GGDI values for the NCP and its geomorphic sub-areas in March 2020 were −3.02, −2.76, −3.07,

Quantitative Assessment of Groundwater Drought Events
The obtained quantitative evaluation results of groundwater drought events in the NCP and its geomorphic sub-areas using the GGDI are shown in Figure 4. The groundwater drought trends in the geomorphic sub-areas were consistent with those observed across the entire NCP, showing gradually decreasing trends. In addition, the most severe groundwater drought event in all areas was observed in March 2020. Indeed, the GGDI values for the NCP and its geomorphic sub-areas in March 2020 were −3.02, −2.76, −3.07, and −3.10, respectively. The groundwater drought frequency in the entire NCP was 24.54%, with mild, moderate, severe, and extreme groundwater drought frequencies of 8.33, 8.80, 3.70, and 3.70%, respectively. The groundwater drought frequency in Area I was 26.39%, with mild, moderate, severe, and extreme groundwater drought frequencies of 10.65, 9.72, 3.24, and 2.78%, respectively. The groundwater drought frequency in Area II was 24.54%, with mild, moderate, severe, and extreme groundwater drought frequencies of 7.87, 7.41, 5.09, and 4.17%, respectively. The groundwater drought frequency in Area III was 23.61%, with mild, moderate, severe, and extreme groundwater drought frequencies of 9.72, 6.02, 4.63, and 3.24%, respectively. The highest frequency of groundwater drought events was observed in Area I, followed by Area II and Area III. However, unlike Area I and Area III, which exhibited mainly mild and moderate drought events, Area II mainly experienced severe and extreme groundwater drought events. was 24.54%, with mild, moderate, severe, and extreme groundwater drought frequencies of 7.87, 7.41, 5.09, and 4.17%, respectively. The groundwater drought frequency in Area III was 23.61%, with mild, moderate, severe, and extreme groundwater drought frequencies of 9.72, 6.02, 4.63, and 3.24%, respectively. The highest frequency of groundwater drought events was observed in Area I, followed by Area II and Area III. However, unlike Area I and Area III, which exhibited mainly mild and moderate drought events, Area II mainly experienced severe and extreme groundwater drought events.

Groundwater Drought Event Identification
According to the identification principles reported in Section 2.3, groundwater drought events in the NCP from 2003 to 2021 were identified using GGDI. The results showed that a total of 49 drought events occurred in this region. The duration of each drought event is shown in Figure 5.

Groundwater Drought Event Identification
According to the identification principles reported in Section 2.3, groundwater drought events in the NCP from 2003 to 2021 were identified using GGDI. The results showed that a total of 49 drought events occurred in this region. The duration of each drought event is shown in Figure 5. of 7.87, 7.41, 5.09, and 4.17%, respectively. The groundwater drought frequency in Area III was 23.61%, with mild, moderate, severe, and extreme groundwater drought frequencies of 9.72, 6.02, 4.63, and 3.24%, respectively. The highest frequency of groundwater drought events was observed in Area I, followed by Area II and Area III. However, unlike Area I and Area III, which exhibited mainly mild and moderate drought events, Area II mainly experienced severe and extreme groundwater drought events.

Groundwater Drought Event Identification
According to the identification principles reported in Section 2.3, groundwater drought events in the NCP from 2003 to 2021 were identified using GGDI. The results showed that a total of 49 drought events occurred in this region. The duration of each drought event is shown in Figure 5.  In total, 49 groundwater drought events occurred in the NCP, with the shortest and longest duration of 1 month and 31 months, respectively. In addition, there were 23, 13, and 2 drought events with a duration of 1, 2, and 3 months, accounting for 46.94, 26.53, and 4.08% of the total groundwater drought events, respectively. Whereas 3, 2, 2, and 1 groundwater drought events lasted 4, 5, 7, and 8 months, accounting for 6.12, 4.08, 4.08, and 2.04% of the total groundwater drought events, respectively. In addition, 2 and 1 drought events lasted 14 and 31 months, accounting for 4.08 and 2.04% of the total groundwater drought events, respectively.

Characteristics of the Groundwater Drought Centroid Migration
Based on the DMM construction principles, 11 groundwater drought events from 2003 to 2020 were selected in this study to investigate the groundwater drought centroid migration characteristics in the NCP. The characteristics of these drought events are reported in Table 4. According to Table 4, DE.49 exhibited the highest intensity, longest duration, lowest monthly average minimum GGDI, and highest drought grid number, with 38.17, 31 months, −3.12, and 221, respectively. Therefore, this event was the most serious groundwater drought event over the considered period, followed by DE.29, with drought duration, monthly minimum GGDI, and drought intensity of 14 months, −247, and 11.53, respectively. DE.46 showed the third highest groundwater drought intensity, with the same groundwater drought duration as DE.29 (14 months) and a maximum number of groundwater drought grids of 216, which was second only to DE.49. In addition, it can be seen from Table 4  In this study, the drought centroid of each drought event was obtained based on the calculation method of the above-mentioned drought centroid. In addition, the migration trajectory of each groundwater drought event was determined by connecting them in chronological order ( Figure 6). According to the obtained results, the migration direction of the 11 drought events was mainly southwest -northeast. In contrast, only DE.41 exhibited a southeast-northwest direction. Because the drought migration in the DMM was on a monthly scale, the longer the distance between the drought centroids of two adjacent months, the higher the drought migration rate. Figure 6a shows that the maximum drought migration rate (1.62×10 8 km/month) of all drought events occurred in the fourth to fifth months of DE.49, resulting in the largest overall migration distance. The mid-month drought centroid was more concentrated. In addition, DE.1 exhibited the smallest drought migration distance, followed by DE.48. It can also be seen from Figure 6 that the centroid migration of the 11 drought events was mostly concentrated in Area II. to fifth months of DE.49, resulting in the largest overall migration distance. The midmonth drought centroid was more concentrated. In addition, DE.1 exhibited the smallest drought migration distance, followed by DE.48. It can also be seen from Figure 6 that the centroid migration of the 11 drought events was mostly concentrated in Area II. The drought centroids of DE.7, DE.9, DE.29, and DE.49 were in Area II and Area III, while that of DE.48 was concentrated in Area I. In order to reveal the developments and changes in the drought events more intuitively, the three drought events with the highest drought intensities among the 11 drought events, namely DE.49, DE.29, and DE.46, were visualized. The drought centroids of these three drought events were first connected in a monthly series, then their drought intensities were represented by spheres. The larger the sphere, the more serious the groundwater drought event.
(1) DE.49 DE.49 occurred from June 2018 to December 2020, lasting 31 months. It was the longest drought event with the highest monthly drought intensity. According to Figure 6 and Figure 7, the highest drought intensities of DE.49 were observed mainly in the central part of the NCP (Area II) over the last 10 months (from March 2020 to December 2020) of the entire drought event. The intensity of the entire drought event showed four obvious peaks in the 5th, 15th, 22nd, and 27th months. The highest drought intensity was observed in In order to reveal the developments and changes in the drought events more intuitively, the three drought events with the highest drought intensities among the 11 drought events, namely DE.49, DE.29, and DE.46, were visualized. The drought centroids of these three drought events were first connected in a monthly series, then their drought intensities were represented by spheres. The larger the sphere, the more serious the groundwater drought event.
(1) DE.49 DE.49 occurred from June 2018 to December 2020, lasting 31 months. It was the longest drought event with the highest monthly drought intensity. According to Figures 6 and 7, the highest drought intensities of DE.49 were observed mainly in the central part of the NCP (Area II) over the last 10 months (from March 2020 to December 2020) of the entire drought event. The intensity of the entire drought event showed four obvious peaks in the 5th, 15th, 22nd, and 27th months. The highest drought intensity was observed in the 22nd month (2.62), followed by those in the 27th (2.31), 15th (1.91), and 5th (1.81) months. In addition, the lowest DE.49 drought intensity (0.31) was observed in the 7th month, followed by that in the 4th month (0.51). The entire DE.49 period showed slightly different drought grid numbers and drought intensities, without exhibiting any apparent correlations between them. The number of drought grids fluctuated greatly in the first 9 months, then remained basically at a high value before substantially decreasing in the 21st month. The results showed three minimum drought grid numbers over the DE.49 drought event, which were 71, 159, and 196 in the 4th, 7th, and 21st months, respectively.
months. In addition, the lowest DE.49 drought intensity (0.31) was observed in the 7th month, followed by that in the 4th month (0.51). The entire DE.49 period showed slightly different drought grid numbers and drought intensities, without exhibiting any apparent correlations between them. The number of drought grids fluctuated greatly in the first 9 months, then remained basically at a high value before substantially decreasing in the 21st month. The results showed three minimum drought grid numbers over the DE.49 drought event, which were 71, 159, and 196 in the 4th, 7th, and 21st months, respectively. (2) DE.29 DE.29 occurred from August 2012 to September 2013, lasting 14 months. According to Figures 6 and 8, the highest monthly groundwater drought intensities were observed mainly in the Bohai Bay area of the NCP in the 4-9th months. The drought intensity of the entire drought event experienced two obvious peaks of 1.97 and 0.75 in the 4th and 12th months, respectively. However, the lowest groundwater drought intensity (0.33) was observed in the 14th month, followed by that in the 10th month (0.37). In addition, the results showed that DE.29 had slightly different monthly groundwater drought grid numbers and intensities. The highest drought grid numbers were observed in the 3rd and 13th months, without presenting high drought intensities. The groundwater drought grid number in the month with the highest groundwater drought intensity was 10, while that in the month with the second highest groundwater drought intensity was 18. It can be seen that small and large groundwater drought areas exhibited high and low drought intensities, respectively.  Figures 6 and 8, the highest monthly groundwater drought intensities were observed mainly in the Bohai Bay area of the NCP in the 4-9th months. The drought intensity of the entire drought event experienced two obvious peaks of 1.97 and 0.75 in the 4th and 12th months, respectively. However, the lowest groundwater drought intensity (0.33) was observed in the 14th month, followed by that in the 10th month (0.37). In addition, the results showed that DE.29 had slightly different monthly groundwater drought grid numbers and intensities. The highest drought grid numbers were observed in the 3rd and 13th months, without presenting high drought intensities. The groundwater drought grid number in the month with the highest groundwater drought intensity was 10, while that in the month with the second highest groundwater drought intensity was 18. It can be seen that small and large groundwater drought areas exhibited high and low drought intensities, respectively.

Discussion
According to the results of our drought event assessment of the NCP over the 2003-2020 period, the most intense groundwater drought event occurred in the June 2018-December 2020 period (DE.  [27], which could have been due to two main reasons. The first reason was related to the difference in the scales of the study areas. Indeed, Wang et al. [27] considered the entire Henan Province and Shandong Province in their groundwater drought assessment of the NCP, while only some parts of Henan Province and Shandong Province were considered in the current study. According to the results of Wang et al. [27], the groundwater drought events in the Henan and Shandong regions were more severe than those in other regions during the August 2013-September 2014 period. However, these regions were not considered in the present study. In fact, Wang et al. [27] assessed the severity of GGDI-based groundwater drought events using the average intensity observed across the entire study area, which also explained the different results obtained. The second reason was related to the identification scale of the groundwater drought events. Indeed, the groundwater drought events in the NCP were assessed from a spatiotemporal perspective in this study, while Wang et al. [27] only assessed temporal variations in the groundwater drought events, which explained the different conclusions obtained.
The NCP is one of the largest groundwater funnels in the world. The Chinese government is committed to improving the local hydrogeological conditions by implementing several appropriate policies. The South-to-North Water Diversion Project is a large-scale water conservancy project established by the Chinese government to alleviate pressure on water resources in North China. Indeed, the Central South-North Water Diversion Project opened to the North China Plain in December 2014 [51]. In addition, relevant departments have implemented water-saving pressure mining policies and carried out ecological water replenishment, which has played an important role in raising the groundwater level in the NCP [38]. However, the results of this study revealed that 7 of the 49 groundwater drought events in the NCP occurred after the implementation of the South-to-North Water Diversion Project and 5 of the 11 droughts that lasted over 3 months occurred after the implementation of the Central South-North Water Diversion Project. Although these measures have raised the groundwater level in the NCP [34,38], the groundwater drought events caused by long-term groundwater overexploitation cannot be mitigated in the short term. Indeed, continuous improvements to water replenishment and water-saving pressure mining measures can enhance groundwater resilience to climate warming. A study by Yang Huifeng et al. showed that groundwater levels in the urban areas of the North China Plain rose significantly after the South-North Water Diversion was opened but water levels in agricultural areas still showed a continuous decline [37] and those agricultural areas in the North China Plain accounted for a large proportion of the GGDI and dominated the regional average. Therefore, although a small recovery in groundwater levels has occurred in some areas of the North China Plain, according to reported studies, no immediate effects on groundwater droughts caused by long-term groundwater overexploitation have been observed, and only by continuously improving water replenishment measures and continuing to carry out water conservation and suppression work will the groundwater systems be able to withstand the test of climate warming.
In this study, a new method for identifying groundwater drought events was proposed from a spatiotemporal perspective to improve our understanding of groundwater droughts. Indeed, the constructed drought migration model may have a beneficial effect on our understanding of the developments and changes in groundwater droughts. The results of our drought migration assessment showed that the centroids of the drought events were mostly concentrated in the deep funnel area [32], indicating that the center of gravity of a groundwater drought is in the deep funnel area. The groundwater drought migration was in the southwest-northeast and northeast-southwest directions, which was the same direction as the slope of the terrain in the study area. There is still an academic gap in our knowledge of the mechanisms of drought migration characteristics. Therefore, we could not explain this phenomenon. Indeed, migration mechanisms may play a crucial role in the prediction and early warning of groundwater drought events. Therefore, future studies should assess the migration mechanisms of groundwater drought events.

Conclusions
This paper used GGDI to quantitatively assess groundwater drought events in the North China Plain, proposed a spatiotemporal groundwater drought event identification method, introduced the center of mass theory, constructed a groundwater drought center of gravity migration model, analyzed the migration characteristics of groundwater drought events in the North China Plain, and conducted a visual mapping study of the three groundwater drought events with the greatest drought intensities. The following conclusions were drawn from this study.
(1) The validation results of the GRACE data showed significant Pearson correlation coefficients (p < 0.01) between the changes in the GRACE-and GLDAS-based water storage in the NCP and its geomorphological areas, ranging from 0.322 to 0.552. Therefore, the GRACE-based results were reliable and could be used to effectively investigate groundwater drought events in the NCP. (2) The groundwater drought frequencies in the NCP, Area I, Area II, and Area III over the 2003-2020 period were 24.54, 26.39, 24.54, and 23.61%, respectively. Although Area I showed a higher groundwater drought frequency than the other sub-areas, mild and moderate groundwater drought events were the most prevalent. In addition, Area II showed a high frequency of moderate groundwater drought events, but its severe and extreme groundwater drought frequencies were relatively higher than those in the other sub-areas. (3) According to the new identification principle for groundwater drought events, 49 groundwater drought events were identified in the NCP over the 2003-2020 period. The maximum duration of drought was 31 months and the minimum was 1 month. Drought events with a drought duration of 1 month were the most frequent, accounting for 46.94% of the total number of drought events, followed by those with a drought duration of 2 months, accounting for 26.53% of the total. Drought events with a drought duration of 8 or 31 months were the least frequent, both accounting for 2.04% of the total. The obtained results indicated that DE.49 was the most severe groundwater drought event, with a drought intensity, duration, and grid number of 38.17, 31 months, and 221, respectively. Meanwhile, DE.29 was the second most intense groundwater drought event, with a drought intensity, duration, and grid number of 11.53, 14 months, and 167, respectively. (4) A total of 11 groundwater drought events were selected from the 49 drought events to construct a drought migration model. The migration direction of 10 of the groundwater drought events was southwest-northeast, which was in line with the slope of the NCP. However, only DE.41 exhibited a southeast-northwest migration direction. The centroids of the groundwater drought events were mostly concentrated in Area II. In summary, the principles for identifying groundwater drought events proposed in this study could reveal regional groundwater drought changes from both temporal and spatial perspectives. In addition, the groundwater drought migration model proposed based on these principles could use the center of gravity and migration direction of groundwater drought events to capture their migration characteristics. The principles for defining groundwater drought events, as well as the groundwater drought migration model, were not only applicable to the NCP but could also be used in other regions, thus providing an in-depth understanding of groundwater development laws and promoting the sustainable use of groundwater resources. However, due to the constraints of research time and data, this paper did not conduct a more in-depth study of groundwater drought migration mechanisms or the causes of such migration characteristics. However, the study of groundwater drought migration mechanisms is important to fully grasp groundwater drought development patterns. Therefore, this would be a good direction for future research.