Impacts of Regional Groundwater Flow and River Fluctuation on Floodplain Wetlands in the Middle Reach of the Yellow River

: Floodplain wetlands are of great importance in the entire river and ﬂoodplain ecosystems. Understanding the hydrological processes of ﬂoodplain wetlands is fundamental to study the changes in wetlands caused by climate change and human activities. In this study, ﬂoodplain wetlands along the middle reach of the Yellow River were selected as a study area. The hydrological processes and the interactions between the river and the underlying aquifer were investigated by combining remote sensing, hydraulic monitoring, and numerical modeling. Wetland areas from 2014 to 2019 were extracted from Landsat 8 remote sensing images, and their correlation with the river runo ﬀ was analyzed. The results indicate that the river ﬂow had a limited impact on the wetland size and so did groundwater levels, due to the strong reliance of wetland vegetation on water levels. Based on hydrological and hydrogeological conditions, a surface water–groundwater coupled numerical model was established. The comparison and correlation analysis between the monitored groundwater head and the simulated river stage also show that river ﬂow did not play a ﬁrst-order role in controlling the groundwater levels of wetlands in the study area. The simulation results also suggest that it is the regional groundwater ﬂow that mainly sustains shallow groundwater of ﬂoodplain wetlands in the study area. The ﬂoodplain wetland of the study area was dynamic zones between the regional groundwater and river, the contrasting pattern of hydrological regimes on both banks of the Yellow River was due to a combination of regional groundwater ﬂow and topography.


Introduction
Wetlands play a vital role in sustaining global ecosystems for their significant contribution to biodiversity support, water quality improvement, flood abatement and carbon sequestrations [1,2]. Over the last decade, wetland coverage has continued to shrink [3,4], leading to widespread concern about changes in the wetland area, vegetation, and water quantity [5][6][7][8]. Understanding hydrological processes is of great importance to study mechanisms driving wetland changes, which is essential for protecting such fragile ecosystems. Floodplain wetlands, as one type of wetlands, occupy nearly one third of the global wetland area and provide far more ecosystem services than the relatively small area [4]. The floodplain wetlands have been regarded as one of the most concerned wetland types for its unique and highly complex hydrological process.
The high complexity of the hydrological process for floodplain wetlands is mainly due to complicated and heterogeneous hydraulic connectivity with adjacent water bodies such as rivers and underlying aquifers [9,10]. Previous studies have proposed the flood pulse concept and pointed out that periodic river flooding is a fundamental force driving the changes in floodplain wetland systems [11,12]. Since then, the majority of published papers have focused on the hydraulic connectivity between rivers and wetlands, including the river's response to flood pulse and its important influence on water as well as the nutrient cycle of floodplain wetland [13,14]. For example, extensive studies have revealed that rivers are the main factor in promoting the exchange of sediments and nutrients, and a major source of water for floodplain wetland [15][16][17][18][19].
Groundwater also plays a significant role in maintaining the health of floodplain wetlands [20]. Shallow aquifer is of great importance for its direct interaction with rivers, and has been found to contribute significant water to wetlands by storing water from rivers in flood seasons [21][22][23][24][25][26][27]. Moreover, evidence accumulated from monitoring data, as well as numerical models, suggest that regional groundwater flow may also exert strong control upon the hydrological process of floodplain wetlands, particularly when the wetlands are underlain by highly permeable aquifers, such as chalk [28,29], or located in groundwater discharge zones [30]. Thus, regional groundwater flow plays a significant role in wetland nutrient status and vegetation distribution [28][29][30][31][32][33][34][35]. Despite the existing studies, how regional groundwater flow affects floodplain wetlands is far from clear. Research is warranted to assess how regional groundwater flow affects floodplain wetlands where the underlying shallow aquifer and the adjacent river interact frequently.
Due to the complexity of hydrological processes in floodplain wetlands, investigating the interactions of river, wetland, and groundwater through field observations alone is often difficult. Therefore, extensive wetland studies have employed numerical simulations instead. Simulations of hydrological processes within fully integrated or coupled groundwater-surface water models have been applied to explore hydrological regimes in wetlands [29,[36][37][38]. However, the studies that simply relied on models are less capable of relating model results to varied patterns of wetlands. In this regard, remote sensing images were introduced by studies to analyze the hydrological regimes in floodplains for their direct display of changes in wetlands [6,7]. Therefore, a more comprehensive understanding of the site's hydrological functioning may be obtained through the combined application of multiple methods.
In this paper, we attempted to examine how regional groundwater flow and river stage changes affect floodplain wetlands. Wetlands along the middle reach of the Yellow River, typical of floodplain wetlands, were selected as the study area. Multiple methods, including remote sensing, hydrological observation and numerical modeling were used to perform this study. Wetland area was first extracted from remote sensing images between 2014 and 2019, then changes in the wetland area with river fluctuation were discussed to understand the connectivity between the wetlands and the adjacent river. Moreover, the measured groundwater heads and river stages at three transects were compared and their correlations were analyzed to infer the connectivity between the shallow groundwater underlying the wetlands and the river. Afterward, a coupled river-groundwater numerical model was established to examine the hydrological processes in relation to the wetlands at three transects.

Study Area
Wetlands along the middle reach of the Yellow River were chosen as our study area ( Figure 1). These wetlands are protected and identified as nature reserves, which cover an area of about 1442 km 2 [39]. The river that runs through the reserves is about 384.5 km in length, with the Longmen and Tongguan gauging stations as the upstream and downstream ends, respectively. Approximately 40% of the wetlands are located on the right bank of the river and belong to the Shaanxi Yellow River Wetland Provincial Natural Reserve (SWPNR), whereas the others are situated on the left and regulated by the Shanxi Yuncheng Wetland Provincial Natural Reserve (YWPNR). These two reserves serve as important habitats for migratory birds in the central and western regions of China. These wetlands consist of mainly permanent river wetlands and floodplain wetlands (riverine wetlands cover 74.78% of the total area). However, these wetlands are not connected to the river directly due to constructed dikes.
Water 2020, 12, x FOR PEER REVIEW 3 of 18 regulated by the Shanxi Yuncheng Wetland Provincial Natural Reserve (YWPNR). These two reserves serve as important habitats for migratory birds in the central and western regions of China. These wetlands consist of mainly permanent river wetlands and floodplain wetlands (riverine wetlands cover 74.78% of the total area). However, these wetlands are not connected to the river directly due to constructed dikes. Figure 1. Location of the study area: (A) shows the geographic location of the study area within the Yellow River (marked in the red square); (B) is the close-up of the study area; and (C) shows the three groundwater monitoring transects that this study focused on.

Hydrology and Meteorology
According to historical records, the mean annual runoff of the Yellow River at the Longmen gauging station is 1060 m 3 /s, with the maximum runoff at 21,000 m 3 /s and the minimum runoff at 53.2 m 3 /s. There have been 25 floods (≥10,000 m 3 /s) since 1950, most of which occurred in July or September [40,41]. The river width increased from several hundred meters to several kilometers during these floods. However, the wetlands were not widely inundated by floods due to the constructed dikes.
The mean annual precipitation for the SWPNR on the right bank is around 545 mm, and approximately 51.6% of the precipitation occurs during summer. In comparison, the mean annual precipitation for the YWPNR is around 635 mm. The mean annual potential evapotranspiration of the entire study area is estimated at 1200 mm.

Hydrology and Meteorology
According to historical records, the mean annual runoff of the Yellow River at the Longmen gauging station is 1060 m 3 /s, with the maximum runoff at 21,000 m 3 /s and the minimum runoff at 53.2 m 3 /s. There have been 25 floods (≥10,000 m 3 /s) since 1950, most of which occurred in July or September [40,41]. The river width increased from several hundred meters to several kilometers during these floods. However, the wetlands were not widely inundated by floods due to the constructed dikes.
The mean annual precipitation for the SWPNR on the right bank is around 545 mm, and approximately 51.6% of the precipitation occurs during summer. In comparison, the mean annual precipitation for the YWPNR is around 635 mm. The mean annual potential evapotranspiration of the entire study area is estimated at 1200 mm.

Geology and Hydrogeology
Quaternary strata cover the entire study area near the land surface. The Quaternary stratum is mainly composed of fine sand, a mixture of fine and medium sand to a mixture of fine sand and clay. On the right bank of the Yellow River, the Quaternary stratum is underlain by the Cambrian-Ordovician carbonate rocks [42]. On the left bank of the Yellow River, the Quaternary stratum overlies the Cambrian-Ordovician carbonate rocks in the north and Archaean metamorphic rocks in the south [43,44].
Groundwater in the study area can be divided into two types: Quaternary pore water that is mainly distributed in the loose deposits near the surface, and karstic water within the underlying Cambrian-Ordovician carbonate rocks. The Quaternary aquifer has a thickness from 60 m to 150 m for the right bank [42,45] and a thickness from 300 to 500 m for the left bank of Yellow River [43,46,47], respectively. Regional groundwater flow in the Quaternary aquifer is directed from northwest to southeast on the right bank and from northeast to southwest on the left, respectively, primarily driven by the topography. The Quaternary aquifer is recharged mainly by precipitation infiltration, lateral flow or upward leakage from the underlying aquifers, and episodic recharge from the Yellow River during flooding. Groundwater discharge occurs in the form of evaporation, springs, and lateral outflow to the Yellow River.
Karstic water exists within structural and corrosive fissures of the Ordovician carbonate rocks with a thickness of up to 1000 m [42,48]. In the northwest of the study area, the karstic stratum is outcropped at the ground surface due to high elevation and receives precipitation directly as groundwater recharge through sinkholes. In other parts of the study area, the karstic aquifer is mostly covered by the Quaternary sediments. For the right bank of the Yellow River, the regional groundwater flow in the karstic aquifer is from northwest to southeast. The discharge of the karstic groundwater is thought to be springs (upward flow) in the Shaanxi Yellow River wetlands (e.g., the Virgin Spring or the Fen Spring).

Remote Sensing and In-Situ Monitoring
In this study, 12 periods (a total of 24 images) of Landsat 8 images from 2014 to 2019 were collected from the websites of the United States Geological Survey (USGS) (http://glovis.usgs.gov/). Two periods of each year were selected to analyze the change of wetlands area: February or March representing dry seasons, and July representing wet seasons. Pre-processing steps (including Geometric correction and radiance calibration) were conducted in ENVI 5.3 to acquire the atmospherically corrected data. After three-band red-green-blue images were stacked in each image, two sensed images in each period were mosaicked to obtain specific images of the study area. Then, each period's wetland area was extracted from the images by classification combining the field investigation and maximum likelihood classifier. In the field investigation, 94 samples were collected, which were used as training data to perform the classification. The overall classification accuracy was greater than 90%.
In order to obtain groundwater head underneath the wetlands, three transects (A, B, C) were set up with a total of six monitoring wells. Each transect has two monitoring wells, which were separately located on the different bank of the river. The groundwater head was monitored daily and the monitoring time started on the 16 July 2018 and finished on the 09 August 2019. The data obtained from these wells were used in numerical modeling as well as data analysis. It should be noted that the monitored data in B2 and C1 only lasted until the 14 March 2018, when these two wells were damaged. In addition, the measured hydrological data, including river stages and flow rates at Longmen and Tongguan, were collected for the following analysis and simulation.

Numerical Modeling
In this study, we established a physically-based numerical model that coupled surface water and groundwater flow to examine river-aquifer interactions and regional groundwater flow. The numerical model was built with an integrated hydrological modelling system MIKE by DHI (Danish Hydraulic Institute, 2014). The numerical code MIKE SHE in the modeling system was used to describe the groundwater flow, and another code, Mike 11, was adapted to describe the river flow. Then, the two models were integrated with the coupling module within MIKE SHE.

Model Setup
The modelled area covered the entire wetlands as shown in Figure 2. A digital elevation map with a spatial resolution of 90 m was obtained from SRTM (Shuttle Radar Topography Mission) to describe the surface topography.
Water 2020, 12, x FOR PEER REVIEW 5 of 18 In this study, we established a physically-based numerical model that coupled surface water and groundwater flow to examine river-aquifer interactions and regional groundwater flow. The numerical model was built with an integrated hydrological modelling system MIKE by DHI (Danish Hydraulic Institute, 2014). The numerical code MIKE SHE in the modeling system was used to describe the groundwater flow, and another code, Mike 11, was adapted to describe the river flow. Then, the two models were integrated with the coupling module within MIKE SHE.

Model Setup
The modelled area covered the entire wetlands as shown in Figure 2. A digital elevation map with a spatial resolution of 90 m was obtained from SRTM (Shuttle Radar Topography Mission) to describe the surface topography. The subsurface component was characterized as a saturated zone, and the 3D finite-difference Darcy flow method was adapted to calculate the groundwater flow. The hydrogeologic units were determined according to the model conceptualization as described in Section 2.2. The northern and southern sides of the model were specified with inflow boundary conditions due to high elevation. The inflow rates were defined based on hydrogeological investigation reports [47][48][49][50], with a uniform value of 0.006 m 3 /s for the northern boundary and 0.002 m 3 /s for the southern boundary. The western and eastern boundaries were defined as inflow boundaries varying with time, whereby the inflow rates were calculated using 90 groundwater level data in 33 monitor wells in the study area and the surrounding region collected from the National Groundwater Monitoring Project of China Geological The subsurface component was characterized as a saturated zone, and the 3D finite-difference Darcy flow method was adapted to calculate the groundwater flow. The hydrogeologic units were determined according to the model conceptualization as described in Section 2.2. The northern and southern sides of the model were specified with inflow boundary conditions due to high elevation. The inflow rates were defined based on hydrogeological investigation reports [47][48][49][50], with a uniform value of 0.006 m 3 /s for the northern boundary and 0.002 m 3 /s for the southern boundary. The western and eastern boundaries were defined as inflow boundaries varying with time, whereby the inflow rates were calculated using 90 groundwater level data in 33 monitor wells in the study area and the surrounding region collected from the National Groundwater Monitoring Project of China Geological Survey. The western inflow boundary had a mean flow rate of 0.08 m 3 /s, a maximum value of 1.2 m 3 /s and a minimum of 0.003 m 3 /s. The eastern inflow boundary had a mean flow rate of 0.054 m 3 /s, a maximum value of 0.7 m 3 /s and a minimum of 0.001 m 3 /s. The model surface representing the maximum wetland was specified with the drainage boundary condition, which allows groundwater to discharge when water table depth is less than the drainage level from the land surface. The drainage level and the drainage time constant were determined during model calibration.
Daily rainfall data in Yuncheng (60 km away from the Yellow River) were obtained from NOAA(National Oceanic and Atmospheric Administration), with a maximum value of 53.6 mm/d, and a mean value of 1.2 mm/d. The extraction was uniformly defined as 160 mm/a and directly subtracted from the precipitation. The Net Rainfall Fraction in MIKE SHE, which represents the proportion of rainfall that infiltrated the aquifer after evapotranspiration, was defined as 0.14~0.32 according to previous work [46,47,51]. Other initial hydrogeological parameters (including the hydraulic conductivity, specific yield and specific storage) were also obtained from the hydrogeological reports [48,[50][51][52]. The initial hydraulic conductivity of the Quaternary aquifer ranged from 6 to 25 m/d and the specific yield was approximately between 0.02 and 0.28, the initial hydraulic conductivity of karstic aquifer was 5.6 m/d and the specific storage was 0.035. Hydraulic conductivities and the specific storages of the Quaternary aquifer and karstic aquifer were varied during the model calibration.
The surface component included the overland flow and river. The overland flow was uniformly distributed, with Manning's roughness coefficient defined as a calibration parameter. The Yellow River was digitized in MIKE 11, which extends from the Longmen gauging station at the upstream end to the Tongguan gauging station at the downstream end. The fully dynamic one-dimensional Saint-Venant equations were employed to describe the channel flow. The measured daily average flow rates at the Longmen gauging station were specified to the upstream end of the surface model as an inflow boundary condition. The measured daily average river stages at the Tongguan gauging station were specified to the downstream end of the model with a water level boundary condition. In addition, the river shape within study area and the distribution of river bed resistance were interpolated based on the measured data from 9 hydrologic cross-sections.
In this paper, all MIKE 11 branches were specified to be coupled to MIKE SHE. Considering the complex interaction between the river and the aquifer (the river could possibly be both recharge and discharge), the channel bed was specified to be in full contact with the saturated zone, so that the exchange between the river and aquifer was controlled by both the hydraulic conductivity of the aquifer and the river bed material.
The numerical model was discretized into 400 rows and 435 columns with an elemental size of 362 m × 362 m. The simulation period was from January 2014 to July 2019. Given that the exact initial condition was unknown, the period between January 2014 and June 2018 was used as the model warm-up period. The period between July 2018 and July 2019, when the groundwater heads were monitored, was then utilized for model calibration. The initial time step was 6 h, and the following time steps were adjusted automatically by MIKE SHE.

Model Calibration
In this paper, the model calibration target was the groundwater heads from monitoring wells, measured river stages at the Longmen gauging stations, and the measured flow rates at the Tongguan gauging stations. Therefore, two objective functions, including the Nash-Sutcliffe coefficient (E f ) [53] and the correlation coefficient (R) for the heads were implemented to optimize the numerical model. The E f and R are given as follows: where obs i,t and sim i,t represent the observed and simulated heads, respectively, at a location i and at time t. n is the number of observations at the location i. An efficiency of 1 (E f = 1) implies a perfect match between the modelled and observed heads, while an efficiency of less than 0 (E f < 0) suggests that the model is not reliable because its prediction is not as accurate as the mean of the observed heads. Moreover, an efficiency of 0 (E f = 0) indicates that the model yields prediction as accurate as the mean of the observed heads is. Model performance criteria were adapted from Henriksen et al.  Table 1 in this paper).

Remote Sensing and In-Situ Monitoring
The remote sensing images between 2014 and 2019 were processed to extract important features that reflect on the changes in the wetland area. In wet seasons, the wetland area was defined as the sum of the continuously inundated area and the vegetated area dominated by reeds ( Figure 3). As the growth of reeds is highly sensitive to water levels [55], the changes of the area covered by the reeds can also reflect on groundwater heads in the shallow aquifer. In dry seasons, the wetland area was extracted from the continuously inundated area only (Figure 4). This is because the reeds have withered by the dry season of a year and cannot be used to reflect on water depth anymore.
According to Figure 5, the annual runoff declined from 2014 to 2016, but the wetland area increased instead. The wetland area reached the maximum value in 2016, whereas the annual runoff at Longmen dropped to its lowest value. Moreover, the annual runoff at Longmen was quite stable in 2016-2017. On the contrary, the wetland area decreased, suggesting a decrease in the wetland water level. In 2017-2018, the annual runoff, as well as the runoff in the wet season rose significantly, but the wetland area only increased lightly. In 2019, the continuous increase in the annual runoff may result in an increase in the wetland area during dry season. However, the wetland area in the wet season fell dramatically. Overall, it can be concluded that the river runoff and the wetland area in wet season are weakly correlated. This weak correlation indicates that the regional groundwater flow may play a more critical role. that reflect on the changes in the wetland area. In wet seasons, the wetland area was defined as the sum of the continuously inundated area and the vegetated area dominated by reeds ( Figure 3). As the growth of reeds is highly sensitive to water levels [55], the changes of the area covered by the reeds can also reflect on groundwater heads in the shallow aquifer. In dry seasons, the wetland area was extracted from the continuously inundated area only (Figure 4). This is because the reeds have withered by the dry season of a year and cannot be used to reflect on water depth anymore.   According to Figure 5, the annual runoff declined from 2014 to 2016, but the wetland area increased instead. The wetland area reached the maximum value in 2016, whereas the annual runoff at Longmen dropped to its lowest value. Moreover, the annual runoff at Longmen was quite stable in 2016-2017. On the contrary, the wetland area decreased, suggesting a decrease in the wetland water level. In 2017-2018, the annual runoff, as well as the runoff in the wet season rose significantly, but the wetland area only increased lightly. In 2019, the continuous increase in the annual runoff may result in an increase in the wetland area during dry season. However, the wetland area in the wet season fell dramatically. Overall, it can be concluded that the river runoff and the wetland area in wet season are weakly correlated. This weak correlation indicates that the regional groundwater flow may play a more critical role. water level. In 2017-2018, the annual runoff, as well as the runoff in the wet season rose significantly, but the wetland area only increased lightly. In 2019, the continuous increase in the annual runoff may result in an increase in the wetland area during dry season. However, the wetland area in the wet season fell dramatically. Overall, it can be concluded that the river runoff and the wetland area in wet season are weakly correlated. This weak correlation indicates that the regional groundwater flow may play a more critical role. Observed groundwater heads were presented together with simulated heads in Section 4.2 below. Observed groundwater heads were presented together with simulated heads in Section 4.2 below.

Model Calibration
The numerical model was calibrated by comparing the simulated groundwater heads, river flow rates and river stages to the observed ones using R and E f . A sensitivity analysis was also implemented in the calibration process. The MIKE Zero automatic calibration procedure was used with individually varied parameters. The model performance was most affected by the river bed resistance and hydraulic conductivities of the Quaternary aquifer. These sensitive parameters were then listed in Table 2. The calibrated model reproduced the observed surface water and groundwater data reasonably well (Figures 6 and 7).
As demonstrated in Figure 6, the changes in the river flow rate and the river stage during the calibrated period was well captured by the model. According to the performance criteria (Table 1), both E f and R suggested that the model performed well in simulating the Yellow River. The R and E f for the river stage at the Longmen station were 0.96 and 0.87, respectively, indicating excellent model performance. The R for the flow rates at the Tongguan station represented excellent performance, and the E f suggested a very good performance. Despite the model's satisfied performance for the river, we found the model's ability to match peak flows in the flood season was not as good as that in the rest of the year. At the Longmen station, the river stages were slightly overestimated in 07.2018-09.2018, then underestimated during 07.2019. At the Tongguan station, the flow rates were slightly underestimated during the same period. The discrepancies between the simulated and observed channel flow values may be in part due to the river bed resistance varied with time. The river bed resistance is affected by many factors, amongst which the in-stream vegetation contributed most. Besides the channel morphology, bed material and adjacent conditions, the varied seasons is the most important factor for the growth and distribution of in-stream vegetation. The seasonal effects of vegetation growth on river bed resistance were not accounted for within the model; instead, a time constant resistance factor was applied throughout the MIKE 11. with time. The river bed resistance is affected by many factors, amongst which the in-stream vegetation contributed most. Besides the channel morphology, bed material and adjacent conditions, the varied seasons is the most important factor for the growth and distribution of in-stream vegetation. The seasonal effects of vegetation growth on river bed resistance were not accounted for within the model; instead, a time constant resistance factor was applied throughout the MIKE 11. The simulated groundwater heads also fitted reasonably well with the observed ones ( Figure 7). Overall, the simulated data reproduced the trend that groundwater heads decline along with the river downstream. For the , the values were between 0.44 and 0.89. Except for B1 and C1, the model performance in other the wells was fair to good. The R for the groundwater heads averaged 0.81, suggesting that the model performance was very good. The model parameters obtained for the best fit are shown in Table 2.
As illustrated in Figure 7, the model results were better for monitoring wells at transect A and C compared with those at transect B. The mean value of R and for transect A were 0.90 and 0.53, for transect C were 0.79 and 0.68, while the mean value of R and for transect B were 0.75 and 0.47. The locations of three monitoring transects may contribute to the difference in model performance. The transect B was sited in the middle of the simulated river reach, where the hydrological data (river The simulated groundwater heads also fitted reasonably well with the observed ones ( Figure 7). Overall, the simulated data reproduced the trend that groundwater heads decline along with the river downstream.
For the E f , the values were between 0.44 and 0.89. Except for B1 and C1, the model performance in other the wells was fair to good. The R for the groundwater heads averaged 0.81, suggesting that the model performance was very good. The model parameters obtained for the best fit are shown in Table 2.
As illustrated in Figure 7, the model results were better for monitoring wells at transect A and C compared with those at transect B. The mean value of R and E f for transect A were 0.90 and 0.53, for transect C were 0.79 and 0.68, while the mean value of R and E f for transect B were 0.75 and 0.47. The locations of three monitoring transects may contribute to the difference in model performance. The transect B was sited in the middle of the simulated river reach, where the hydrological data (river stage and flow rate) were calculated based on the measured data at gauging stations. By the contrast, transect A and C sited near the Longmen station and the Tongguan station, respectively. The measured hydrological data at these gauging stations restrained the interaction between the groundwater and the river at transects A and C, obtaining better fitted groundwater heads.
Generally, the localized and time-varying effect stems from the parameters of channel flow (river bed resistance, river stage and flow rate at cross-sections), which could therefore undermine the model's performance.
For transect A, the river stage rose from July to October, then came down in the next months. The maximum river stage in the monitored period was 375.82 m (on 23 September 2018), and the minimum was 372.93 m (on December 2018). The maximum head of Well A1 was 369.19 m (on 23 October 2018), the minimum value was 368.19 m (on 16 July 2018), with a difference of less than 1 m. This small difference suggests that the groundwater at Well A1 was rather stable within a year. The groundwater head at A2 varied in a similar fashion to that at A1 from July 2018 to March 2019. However, the groundwater head at A2 did not rise along with the river stage after April 2019, unlike A1.
For transect B, the groundwater head exhibited a similar temporal trend with the river stage, but the variation range was smaller: 3.17 m for the river stage, 0.99 m for Well B1, and 0.84 m for Well B2. Moreover, the river stage was lower than the groundwater head of Well B1 during 15-24 December 2018.
The river stage at transect C was more stable than transects A and B, with a seasonal fluctuation of less than 2 m. The highest groundwater heads were observed in winter (December) and the lowest in summer (August for C1 and June for C2), opposite to the river stage. It should be noted that C2 had the most stable groundwater head among all the monitoring wells, with a variation range of 0.63 m.  A linear correlation analysis was conducted with the Pearson correlation coefficient (Pearson's R) to investigate the relationship between the groundwater heads and the Yellow River (Figure 8). The results showed that the groundwater heads at transect B were most correlated with the river stage among three transects. The groundwater heads for all the monitoring wells were largely correlated with the corresponding river stage except A1. In addition, the groundwater heads of the left bank had a higher correlation coefficient than the right.

Connectivity between the River and the Wetland Groundwater
A linear correlation analysis was conducted with the Pearson correlation coefficient (Pearson's R) to investigate the relationship between the groundwater heads and the Yellow River (Figure 8). The results showed that the groundwater heads at transect B were most correlated with the river stage among three transects. The groundwater heads for all the monitoring wells were largely correlated with the corresponding river stage except A1. In addition, the groundwater heads of the left bank had a higher correlation coefficient than the right.
Water 2020, 12, x FOR PEER REVIEW 12 of 18 Figure 8. Observed groundwater heads, simulated river stages and their correlations at three transects. River stages on the groundwater monitoring transects were acquired from the numerical model.
The river stage was correlated with the groundwater head in the wetlands according to the linear correlation analysis results in Figure 8. Combined with a comparison of water levels, it could be concluded that there was supposed to be some recharge from the river to groundwater in the wetlands as the direction of the hydraulic gradient showed most of the time. What should be noted is that the correlations did not suggest that the groundwater head in the wetlands was controlled, or severely affected by river flow. The apparent differences in the temporal patterns of groundwater head and river stage indicated that the groundwater head in the wetland was probably not mainly controlled by river flow. The correlation analysis results also showed a difference in the river-wetland The river stage was correlated with the groundwater head in the wetlands according to the linear correlation analysis results in Figure 8. Combined with a comparison of water levels, it could be concluded that there was supposed to be some recharge from the river to groundwater in the wetlands as the direction of the hydraulic gradient showed most of the time. What should be noted is that the correlations did not suggest that the groundwater head in the wetlands was controlled, or severely affected by river flow. The apparent differences in the temporal patterns of groundwater head and river stage indicated that the groundwater head in the wetland was probably not mainly controlled by river flow. The correlation analysis results also showed a difference in the river-wetland connections for both banks, which implies that the hydrological processes may differ for wetlands on both banks.
The connectivity between wetland groundwater and river was also explored based on the simulation results. Based on the water balance results from the simulation, the proportional contribution of the river to the Quaternary aquifer within the wetland area was discussed. Because of the possible difference of river-wetland connection implied by correlation analysis, simulation results for both banks were calculated, respectively.
For the wetlands on the right bank, the rainfall and recharge from the regional groundwater flow (including the lateral inflow and the upward flow) took up most of the recharge of the target aquifer (21.8% and 67.7%, respectively). Water from the river only accounted for 10.4% of the total recharge. For the left bank, the river seepage took about 26% of the total recharge, followed by precipitation (about 24.1%). Recharge from the groundwater flow was approximately 49.9%. These data suggested that instead of the river, it is the regional groundwater flow that played as the most crucial source of water for the wetlands on both banks, sustaining the groundwater of the wetland.

Difference in the Hydrological Process on Both Banks
Results from the remote sensing data, observed data, and the numerical simulation suggested that the Yellow River had a limited impact on groundwater heads within the shallow aquifer and that the regional groundwater flow plays a more important role. These results also indicated the difference in the hydrological processes of wetlands on both banks. Groundwater on the left was more correlated with the river than on the right bank ( Figure 8). The river-aquifer exchange flux at three transects was plotted in Figure 9.
Water 2020, 12, x FOR PEER REVIEW 13 of 18 connections for both banks, which implies that the hydrological processes may differ for wetlands on both banks. The connectivity between wetland groundwater and river was also explored based on the simulation results. Based on the water balance results from the simulation, the proportional contribution of the river to the Quaternary aquifer within the wetland area was discussed. Because of the possible difference of river-wetland connection implied by correlation analysis, simulation results for both banks were calculated, respectively.
For the wetlands on the right bank, the rainfall and recharge from the regional groundwater flow (including the lateral inflow and the upward flow) took up most of the recharge of the target aquifer (21.8% and 67.7%, respectively). Water from the river only accounted for 10.4% of the total recharge. For the left bank, the river seepage took about 26% of the total recharge, followed by precipitation (about 24.1%). Recharge from the groundwater flow was approximately 49.9%. These data suggested that instead of the river, it is the regional groundwater flow that played as the most crucial source of water for the wetlands on both banks, sustaining the groundwater of the wetland.

Difference in the Hydrological Process on Both Banks
Results from the remote sensing data, observed data, and the numerical simulation suggested that the Yellow River had a limited impact on groundwater heads within the shallow aquifer and that the regional groundwater flow plays a more important role. These results also indicated the difference in the hydrological processes of wetlands on both banks. Groundwater on the left was more correlated with the river than on the right bank ( Figure 8). The river-aquifer exchange flux at three transects was plotted in Figure 9. Exchange flux between the river and the wetland at points during the simulation period (positive value represented that the groundwater is recharging the river, negative value indicated that the river supplies groundwater); inflow flux of both banks from the regional groundwater flow.
According to the simulation results for the right bank, the river recharge only accounted for 10.4% of the total recharge for the shallow aquifer. The primary discharge of the wetlands was evaporation and extraction (79.2%), whereas drainage to the river occupied the second place (around 16%), and the lateral outflow only accounted for about 4.8%. This indicated that the floodplain wetlands in the Figure 9. Exchange flux between the river and the wetland at points during the simulation period (positive value represented that the groundwater is recharging the river, negative value indicated that the river supplies groundwater); inflow flux of both banks from the regional groundwater flow.
According to the simulation results for the right bank, the river recharge only accounted for 10.4% of the total recharge for the shallow aquifer. The primary discharge of the wetlands was evaporation and extraction (79.2%), whereas drainage to the river occupied the second place (around 16%), and the lateral outflow only accounted for about 4.8%. This indicated that the floodplain wetlands in the right bank was a dynamic environment where drainage to the river and recharge from the river both occurred. The spatial and temporal variation in the exchange flux of the right bank (R1, R2 and R3 shown in Figure 9) also supported this conclusion. The exchange flux of R1 suggested that the groundwater was gaining water from the river at all times, while R2 and R3 implied that the groundwater was losing water to the river most of the year. A temporal variation could be observed in R2 and R3: the wetland groundwater heads were generally higher than the river stage and drained to river, but the hydraulic gradient was reversed during the wet season, and instead, the river stage rose rapidly and thus supplied water for the wetlands.
For the left bank, river seepage made up over 1/4 of the total recharge (about 26%). Major discharge occurred in the form of evaporation and extraction. The rest drained to the river (about 3.4%). The recharge from the river was far more than the amount drained to the river, implying that although the shallow groundwater in the wetland was mostly supported by rainfall and regional groundwater flow, the river also acted as one of the main sources for the left bank. Negative values of L1, L2 and L3 located on the left bank further demonstrated that the wetland was continuously gaining water from the river during the whole simulation period. The inflow flux from the regional groundwater flow into the wetlands for both banks was also shown in Figure 9, which was generally one order magnitude larger than the exchange flux between the groundwater and the river.
Generally, the floodplain wetland groundwater in the middle reach of Yellow River was the dynamic zones between the regional groundwater and the river. Moreover, as already speculated above, it was the regional groundwater flow that maintained the groundwater quantity of the wetlands. This reflected the larger-scale catchment processes and the position of the Yellow River as well as its floodplain wetland in the discharge zones of regional groundwater flow systems. Regional groundwater flow perhaps also resulted in spatial variations of the interaction between the wetland groundwater and the river. The underlying karstic aquifer in the right bank not only discharged to the river, but probably also generally interconnected with the wetland groundwater by numerous fractures and faults. Especially, the faults-induced karstic springs around Heyang (near groundwater monitoring transect B) continuously supplied water for the wetland, which probably caused the groundwater-river exchange flux difference between transect A and B (see R1 and R2 in Figure 8). By contrast, the wetland groundwater in the left bank has a much larger aquifer thickness than the right bank, so it did not get as much water support from the underlying karstic aquifer. As a result, the wetland on the left bank has lower groundwater heads and thus relies on river recharge more. As the regional groundwater flow direction is majorly controlled by topography in the study area, this contrasting pattern of hydrological regimes on both banks of the Yellow River is probably due to a combination of hydrogeological conditions and topography.

Summary and Conclusions
The hydrological and ecological response mechanism of the floodplain wetland systems to the changed flow regimes of surface water caused by human activity and climate change has long been the research focus worldwide in recent decades. The results are considered uncertain, mainly due to the complex hydrological process of floodplain wetlands, which has not yet been adequately addressed. This study takes regional groundwater flow into account, employing remote sensing images, field data and a coupled surface water-groundwater numerical model to investigate the hydrological process of wetlands distributed along the middle Yellow River (the Shaanxi Yellow River Wetland Provincial Natural Reserve and Shanxi Yuncheng Wetland Provincial Natural Reserve). The wetland areas were extracted from the remote sensing images of 12 periods from 2014-2019 (in the wet season and dry season, respectively), the changes were then compared with a river runoff to explore the possible relationship between water quantity and river flow variety. As a result, the river flow was found not to be the main factor controlling the water quantity of the wetland, and its impact is limited. The in situ monitoring data, including the groundwater heads in wetland monitored by six wells on three monitoring sections cross the river, and the measured hydrological data collected from two key gauging stations (Longmen and Tongguan) were used to establish and calibrate a coupled surface water-groundwater numerical model. The outcomes of the simulation, combined with hydrogeological conditions, were used to investigate the hydrological process in the wetland. The results indicated that the regional groundwater flow is likely to be the most important source to support groundwater quantity in floodplain wetlands in the study area. In addition, the river functions as a discharge for the right bank of Yellow River and a recharge source for the left bank. Thus, the wetland in the left bank relies more on river recharge than the right. These contrasting patterns of the hydrological process in the study area are probably due to a combination of hydrogeological conditions and topography.
The understandings of the hydrological process of wetlands along the middle reach of the Yellow River obtained from this study are critical for following researches in this area, and more monitoring work is required for further detailed research. This study also highlights the importance of taking regional groundwater flow into account when managing floodplain wetlands in groundwater discharge zones.