1.Introduction

As a significant supply source of freshwater resources, groundwater plays a crucial role in social production and human life ^[1][2]. Globally, it provides drinking water for approximately two billion people ^[3] and irrigation for roughly 40% of areas equipped for irrigation ^[4]. Due to extreme climate episodes and anthropogenic actives (e.g., drought and overuse of irrigation water), groundwater resources are seriously over-exploited in some typical regions ^[5][6], leading to a series of environmental issues, such as land subsidence and seawater intrusion ^[7][8]. Therefore, understanding the dynamics changes in groundwater is necessary for the effective utilization and sustainable management of water resources ^[2][9].

As a significant supply source of freshwater resources, groundwater plays a crucial role in social production and human life[1,2].Globally, it provides drinking water for approximately two billion people [3] and irrigation for roughly 40% of areas equipped for irrigation[4].Due to extreme climate episodes and anthropogenic actives (e.g., drought and overuse of irrigation water), groundwater resources are seriously over-exploited in some typical regions [5,6]，leading to a series of environmental issues, such as land subsidence and seawater intrusion [7,8]. Therefore, understanding the dynamics changes in groundwater is necessary for the effective utilization and sustainable management of water resources [2,9].

The traditional method to monitor groundwater levels is based on ground-based measurements ^[10][11]. However, it is not applicable to large-scale and remote regions restricted by national policies, limited stations, and instrument accuracy. The Gravity Recovery and Climate Experiment (GRACE) satellites, successfully launched in March 2002 ^[12], provides a kind of new method for monitoring the global time-variable gravity field with unprecedented accuracy ^[13]. Additionally, it can provide continuous terrestrial water storage anomalies (TWSA) and cover most parts of the world, which is especially beneficial for areas lacking ground-based measurements. By integrating auxiliary information from hydrological models, groundwater storage anomalies (GWSA) can be further isolated from GRACE observations. Previous studies have demonstrated that GRACE missions show great potential in various fields, e.g., detection of terrestrial water storage ^[14][15][16] and groundwater storage ^[17][18]. Swenson et al. ^[14] derived region-scale water storage by applying averaging kernels to a realistic synthetic GRACE gravity signal within North American river basins. Results indicated that the accuracy might be better than 1 cm for regions with 400,000 km

The traditional method to monitor groundwater levels is based on ground-based measurements [10,11]. However, it is not applicable to large-scale and remote regions restricted by national policies, limited stations, and instrument accuracy. The Gravity Recovery and Climate Experiment (GRACE) satellites, successfully launched in March 2002 [12]，provides a kind of new method for monitoring the global time-variable gravity field with unprecedented accuracy [13]. Additionally, it can provide continuous terrestrial water storage anomalies (TWSA) and cover most parts of the world, which is especially beneficial for areas lacking ground-based measurements. By integrating auxiliary information from hydrological models, groundwater storage anomalies (GWSA) can be further isolated from GRACE observations. Previous studies have demonstrated that GRACE missions show great potential in various fields, e.g., detection of terrestrial water storage[14,15,16] and groundwater storage [17,18]。Swenson et al. [14]derived region-scale water storage by applying averaging kernels to a realistic synthetic GRACE gravity signal within North American river basins. Results indicated that the accuracy might be better than 1 cm for regions with 400,000 km

²

or larger areas. Rodell et al. [17] simulated groundwater storage based on GRACE and hydrological modeling system, and the results showed that it was being depleted at a mean rate of 4.0 ± 1.0 cm/yr equivalent height of water over the Indian states of Rajasthan, Punjab and Haryana (including Delhi).

Although remarkable achievements have been made in large-scale areas, the application of GRACE observations in local areas is very limited due to the coarse spatial resolution (~200,000 km

²) ^[19]. Consequently, some downscaling methods have been proposed for providing high-resolution GRACE products, which are mainly divided into two categories: dynamic downscaling and statistical downscaling ^[20][21], respectively. Normally, dynamic downscaling tends to achieve regional downscaling by using the initial boundary conditions of global climate models ^[22][23] directly. For example, Eicker et al. ^[24] assimilated GRACE-derived TWSA into the WaterGAP Global Hydrology Model by introducing a new Kalman filter method, which can provide reasonable results in the Mississippi river basin. Although data assimilation methods remain consistent in the physical process, some shortcomings still require to be considered ^[25]. The implementation of data assimilation is relatively complicated ^[26], and its accuracy is subject to the full error covariance matrix of GRACE observations and hydrological models ^[27][28].

) [19]，Consequently, some downscaling methods have been proposed for providing high-resolution GRACE products, which are mainly divided into two categories: dynamic downscaling and statistical downscaling [20,21], respectively. Normally, dynamic downscaling tends to achieve regional downscaling by using the initial boundary conditions of global climate models [22,23], directly. For example, Eicker et al. [24] assimilated GRACE-derived TWSA into the WaterGAP Global Hydrology Model by introducing a new Kalman filter method, which can provide reasonable results in the Mississippi river basin. Although data assimilation methods remain consistent in the physical process, some shortcomings still require to be considered [25]. The implementation of data assimilation is relatively complicated [26], and its accuracy is subject to the full error covariance matrix of GRACE observations and hydrological models [27,28].

Compared to dynamic downscaling, statistical downscaling usually establishes the linear or non-linear relationships between input and target variables, aiming to produce local-scale information ^[29][30][29,30]. Initially, linear regression models are employed to downscale GRACE products ^[31][32][33][31,32,33]. For example, Ning et al. [31] achieved the downscaling of GRACE data in parts of Yunnan by constructing an empirical regression model based on the water balance equation, and the results indicated the feasibility for downscaling GRACE data. Yin et al. [33] proposed a new statistical downscaling algorithm by building the relationship between multi-source evapotranspiration (ET) products and GWSA in the North China Plain, which obtained desirable downscaled results but limited by the strong correlation between TWS and ET. Practically, the relationships between predictor and predictand tend to be non-linear rather than linear. The development of machine learning algorithms provides effective measures to quantify the complicated relationship by constructing non-linear models. Artificial neural networks (ANNs) have the capabilities of simulating complex hydrological characteristics to an arbitrary degree of accuracy ^[34][35][34,35]. This makes ANN becomes an attractive measure in the downscaling researches, which have been applied to some typical regions, e.g., the Northern High Plains [36], California’s Central Valley [34], and the Lower Peninsula of Michigan [26]. Similarly, some tree-based machine learning algorithms (e.g., random forest (RF) and gradient boosting decision tree (GBDT)) become popular in regression tasks with the advantages of simplicity and effectiveness. The RF algorithm has been utilized to downscale GRACE observations and obtained satisfactory results in some areas ^[37][38][37,38]. As a kind of ensemble machine learning algorithm, GBDT performs well in constructing non-linear regression models, which is often employed to forecast ET [39] and urban flood ^[40][41][40,41], but rarely in GWLA. Furthermore, the multi-stage machine learning algorithm may have more powerful expressive performance than a single algorithm in downscaling GRACE products. For example, Seyoum et al. [42] designed a two-layer boosted regression trees (BRT) model by utilizing GRACE data and hydrological variables in a glacial aquifer system of the United States, which can predict groundwater level anomalies (GWLA) with a high spatial resolution.

2. Study Area and Data

2.1. Study Area

The North China Plain, located in the eastern coastal region of China, lies between latitude 35° N–41° N and longitude 113° E–120° E (Figure 1). It is one of the three great plains in China, covering an area of approximately 140,000 km^{2 [43]}[48]. The NCP is a central agricultural area in China, which produces about one-fourth of the country’s total grain yield [10]. The main crops include winter wheat and summer maize, and the NCP supplies more than 50% of the wheat and approximately 33% of the maize production in China ^[44][49].

Figure 1. Location and digital elevation model map of the NCP. The blue and magenta dots represent the training and testing wells, respectively.

The NCP belongs to a continental monsoon climate with an annual average temperature between 8 and 15 °C ^[33][45]. The annual precipitation, most of which occurs during the growth period of summer maize, ranges from 500 to 600 mm, and annual evaporation is 900–1400 mm ^[46][47]. The NCP contains a shallow unconfined aquifer (40–60 m) and three confined aquifers of different depths (120–170 m, 250–360 m, and 400–600 m) ^[48].

The NCP belongs to a continental monsoon climate with an annual average temperature between 8 and 15 °C [33,50]. The annual precipitation, most of which occurs during the growth period of summer maize, ranges from 500 to 600 mm, and annual evaporation is 900–1400 mm [51,52]. The NCP contains a shallow unconfined aquifer (40–60 m) and three confined aquifers of different depths (120–170 m, 250–360 m, and 400–600 m) [53].

2.2. Data

The fusion model, proposed in this study, is designed based on the water balance principle and machine learning methods (i.e., multiple linear regression and GBDT). Several datasets (terrestrial water storage anomalies, precipitation, runoff, evapotranspiration, soil moisture, snow water equivalent, and groundwater level) are chosen to owe to their close relationship with groundwater storage changes, as shown in Figure 2.

Figure 2. Summary of variable information (resolutions, units, processing, and sources) employed in the fusion model.

Specifically, some variables should be resampled from 0.25° to 1° for matching the spatial availability of GRACE products, and the study period covers from January 2005 to December 2014, with a total of 120 months. The schematic diagram of the water balance principle is shown in Figure 3.

Figure 3. Sketch map of the water cycle.

2.2.1. GRACE TWSA

The GRACE gravity satellites, jointly developed by NASA (National Aeronautics and Space Administration) and DLR (German Aerospace Center), were launched in 2002 and successfully completed its missions in 2017 [13]. They were designed to track global mass changes or gravity variations using the K-band ranging system and low-low satellite tracking satellite mode ^[13][49][13,54]. In this study, GRACE observations are provided by the Jet Propulsion Laboratory (JPL), Center for Space Research (CSR), and GeoForschungsZentrum Potsdam (GFZ), respectively. The gridded-gain factors are utilized to reduce the leakage error ^[50][55], which are available at ^[51][56]. Some discrepancies exist among these three solutions due to different processing strategies and tuning parameters ^[52][57]. Therefore, we utilize the ensemble average of different solutions as the representative TWS estimates in the following discussion.

2.2.2. TRMM Precipitation

The Tropical Rainfall Measurement Mission (TRMM) is a joint project of NASA and Japan Aerospace Exploration Agency (JAXA), aiming to analyze the impact of rainfall data on weather and climate [ ^[53]58]. The monthly precipitation products (2003–2015) used in this study are the TRMM 3b43 with the spatial resolution of 0.25° × 0.25°. The dataset can be obtained from the Goddard Earth Sciences Data and Information Services Center (GES DISC) ^[54][59]. Previous studies [16] have demonstrated that TRMM data match well with gauged stations compared to other remotes sensing products, thus utilized in this study.

2.2.3. GLDAS Data

The Global Land Data Assimilation System (GLDAS) is developed by the Goddard Space Flight Center (GSFC). The primary goal of the GLDAS is to ingest satellite- and ground-based observational data products, using advanced land surface modeling and data assimilation techniques, in order to generate optimal fields of land surface states and fluxes ^[55][60]. Up to now, there have been four land surface models (LSM), namely Mosaic ^[56][61], Community Land Model (CLM) ^[57][62], Noah ^[58][63], and Variable Infiltration Capacity (VIC) ^[59][64], respectively. The Noah model is selected to provide some water-budget variables, including runoff (R), soil moisture (SM), and snow water equivalent (SWE). The runoff includes surface runoff and underground runoff, and soil moisture is the sum of four soil water layers. The datasets from the Noah model include two kinds of temporal resolutions (3-h and monthly scale) and spatial resolutions (0.25° × 0.25° and 1° × 1°). Monthly datasets are chosen in this study, with the resolutions of 0.25° and 1°, which are available at ^[60][65].

2.2.4. GLEAM Product

GLEAM (Global Land Evaporation Amsterdam Model) is a set of algorithms dedicated to estimating global evapotranspiration by combining satellite observations and the Priestley and Taylor equation ^[61][66]. It has been continuously revised and updated since 2011, and the third version of the model was released in 2017 ^[62][67]. The latest version GLEAM v3.3 contains two kinds of data sets (v3.3a and v3.3b), differing in their forcing and temporal coverage ^[63][68] . In this study, the GLEAM v3.3a is employed to provide monthly estimates in evapotranspiration, spanning the 36 years from 1980 to 2018.

2.2.5. Groundwater Level

Monthly groundwater level, collected from the Haihe River Basin Water Resources Bulletin ^[64][69], is used to validate the accuracy of the downscaled TWSA and predicted water level. Groundwater monitoring wells are unevenly distributed across the NCP, and there are abnormal jumps, data gaps, and outliers in some wells. Therefore, the collected data require pre-processing as follows: (1) ignore the wells with more missing months and obvious errors; (2) aggregate the groundwater level to pixels values on 0.25° cells by using the simple average of groundwater observations within the pixel; (3) remove the mean value to obtain the groundwater level anomalies data. A total of 18 wells are selected for this study (Figure 1), and 12 wells are used to train the downscaled model, while the remaining 6 wells are used to test the performance of models.