A Study on the Cause of Layered Seawater Intrusion in the Daqing River Estuary of Liaodong Bay, China

: Groundwater over-pumping in estuary cities leads to a series of groundwater environmental problems that seriously restricts economic development. On the basis of ﬁeld investigation and long-term monitoring data analysis, a three-dimensional numerical model was built in the estuary of the Daqing River in Liaodong Bay, China. The Quaternary overburden can be generalized into ﬁve layers according to particle composition and parameters in the vertical direction. There are many scattered irrigation wells pumping in the second layer, and three water source areas mainly pumping groundwater in the fourth layer. Long-term over-pumping in multi-layered aquifers causes onshore layered seawater intrusion. The laws of layered intrusion under the layered pumping were calculated and analyzed with SEAWAT-2000, and the sensitivity was analyzed with the Sobol method. Results showed that the intrusion area had an obvious layered law. Layered pumping directly a ﬀ ected the layered intrusion area, as di ﬀ erent permeability, tide and barrage further a ﬀ ected it. The prediction study showed that the cone of depression recovered after the pumping-limit of water source areas, and the intrusion area started to retreat in the fourth layer. At that time, the pumping quantity of irrigation wells became the main reason for the increase of the intrusion area. If the water source areas are used to bear part of the irrigation demand, so as to reduce the pressure of pumping in the second layer, the overall intrusion area can be reduced by about 0.23 km 2 under the same pumping quantity.


Introduction
Many large cities in the world are located in the coastal areas and river estuaries. Due to the utilization of groundwater resources for the economic development demands in coastal areas, and sea level rising caused by global warming, the research on the seawater intrusion in estuaries and coastal cities is a hot issue [1][2][3][4][5][6][7]. Seawater intrusion is an important environmental problem that leads to groundwater salinization, freshwater resources reduction, soil salinization and industrial machinery corrosion, endangers human health and seriously restrict the economic development of coastal cities [8][9][10][11][12]. The World Health Organization (WHO) indicated that freshwater is directly nonpotable when it is mixed with about 1% volume of seawater (250 mg/L chloride) [13]. There are two kinds of inducement factors for seawater intrusion, including natural factors (sea level rising, drought and tide) and human factors (groundwater pumping and mariculture) [14][15][16][17][18][19][20][21][22][23]. Compared with the minor equilibrium damage caused by sea level rise, seawater intrusion caused by a large imbalance of hydraulic gradient due to over-pumping groundwater is more important [1,18,[24][25][26][27].

Geological Conditions
Yingkou City is located in the south wing of Yingkou to the Kuandian uplift, an anticline of Liaodong platform, which belongs to the northern China platform [59]. The mountains were formed by the Yanshan movement, and are distributed in the area from east to west. The terrain gradually changes from high to low from east to west. Geomorphology regularly changes under the influence of structure, lithology and geotectonic movement from east to west, namely, low mountains, high hills, low hills and coastal plains [59]. The study area is a mountain valley plain formed by the alluviation and diluviation of the Daqing River, which originates in the eastern mountain of Yingkou. The Liaoning hydrogeology brigade carried out a geological survey of the middle-lower valley plain of Daqing River in 1973. According to the survey, the conceptual model of the aquifer medium in the study area was generalized (Figure 2), and the characteristics of the aquifer in the study area were analyzed.

Geological Conditions
Yingkou City is located in the south wing of Yingkou to the Kuandian uplift, an anticline of Liaodong platform, which belongs to the northern China platform [59]. The mountains were formed by the Yanshan movement, and are distributed in the area from east to west. The terrain gradually changes from high to low from east to west. Geomorphology regularly changes under the influence of structure, lithology and geotectonic movement from east to west, namely, low mountains, high hills, low hills and coastal plains [59]. The study area is a mountain valley plain formed by the alluviation and diluviation of the Daqing River, which originates in the eastern mountain of Yingkou. The Liaoning hydrogeology brigade carried out a geological survey of the middle-lower valley plain of Daqing River in 1973. According to the survey, the conceptual model of the aquifer medium in the study area was generalized (Figure 2), and the characteristics of the aquifer in the study area were analyzed. In the study area, the thickness of Quaternary overburden was greatly affected by the fluctuation of the underlying bedrock surface, ranging from 20 to 65 m [59]. As a whole, the overburden is gradually thinner from the central zone of the valley to the hills area on both sides, and the overburden on the hill in the north area is relatively thick due to the alluvial-proluvial deposit action. River discharge has changed with the seasonal climate in history; therefore, the plain aquifer formed by the river's alluvial-proluvial deposit has great differences that can be divided into five layers from top to bottom.

Characteristic of Daqing Basin
The Daqing River, which originated from the eastern mountain of Yingkou City, is the main surface water in the study area. The width of the riverbed is about 20-300 m. The basin area is about 1468 km 2 . In order to observe the water level, flow and velocity of Daqing River, the Wangbaoshan hydrological station was built in the middle-lower reaches of Daqing River in 1959. According to Wangbaoshan observation data, the average maximal and minimal values of the river flow are 57 and 0.316 m 3 /s, respectively. The flow significantly changes with the seasons. In order to adjust the flow variation with the seasons of Daqing River, the Shimen reservoir was built in the upstream of Daqing River. The flow rate is detected to be about 0.9 m 3 /s at Wangbaoshan hydrological station. In addition, a barrage was built in the lower reaches of Daqing River. It has the functions of storing water, irrigation and blocking tide ( Figure 1). In the upper reaches of Daqing River, the riverbank terrace is narrow and the area is small. The valley plain in the lower reaches of Wangbaoshan gradually opens. The observation data of water level in the hydrological station are complete. Therefore, the reach of the Wangbaoshan hydrological station is considered to be a constant head boundary that can reflect the influence of the upstream river on the downstream area.

Hydrogeological Characteristics of Aquifer System
The thickness of the first layer ranges from 5 to 10 m. The permeability of this layer is relatively low. The hydraulic conductivity of the silt coast is about 0.4 m/day. The aquifer system is unconfined In the study area, the thickness of Quaternary overburden was greatly affected by the fluctuation of the underlying bedrock surface, ranging from 20 to 65 m [59]. As a whole, the overburden is gradually thinner from the central zone of the valley to the hills area on both sides, and the overburden on the hill in the north area is relatively thick due to the alluvial-proluvial deposit action. River discharge has changed with the seasonal climate in history; therefore, the plain aquifer formed by the river's alluvial-proluvial deposit has great differences that can be divided into five layers from top to bottom.

Characteristic of Daqing Basin
The Daqing River, which originated from the eastern mountain of Yingkou City, is the main surface water in the study area. The width of the riverbed is about 20-300 m. The basin area is about 1468 km 2 . In order to observe the water level, flow and velocity of Daqing River, the Wangbaoshan hydrological station was built in the middle-lower reaches of Daqing River in 1959. According to Wangbaoshan observation data, the average maximal and minimal values of the river flow are 57 and 0.316 m 3 /s, respectively. The flow significantly changes with the seasons. In order to adjust the flow variation with the seasons of Daqing River, the Shimen reservoir was built in the upstream of Daqing River. The flow rate is detected to be about 0.9 m 3 /s at Wangbaoshan hydrological station. In addition, a barrage was built in the lower reaches of Daqing River. It has the functions of storing water, irrigation and blocking tide ( Figure 1). In the upper reaches of Daqing River, the riverbank terrace is narrow and the area is small. The valley plain in the lower reaches of Wangbaoshan gradually opens. The observation data of water level in the hydrological station are complete. Therefore, the reach of the Wangbaoshan hydrological station is considered to be a constant head boundary that can reflect the influence of the upstream river on the downstream area.

Hydrogeological Characteristics of Aquifer System
The thickness of the first layer ranges from 5 to 10 m. The permeability of this layer is relatively low. The hydraulic conductivity of the silt coast is about 0.4 m/day. The aquifer system is unconfined aquifer ranges from the second layer to the fifth layer. The second layer is mainly composed of sand gravel with thickness of about 10-25 m. Hydraulic conductivity is between 10 and 40 m/day, which gradually decreases from the valley to both sides. The third layer is mainly composed of sandy loam with thickness of about 0-8 m and hydraulic conductivity of about 5 m/day. The fourth layer is composed of gravel pebble, with thickness of about 0-15 m. The layer gradually becomes thinner from the riverbed to the hills, with hydraulic conductivity between 5 and 70 m/day. The fifth layer is mainly composed of clay containing some gravel. Hydraulic conductivity is about 4 m/day. The second and fourth layers are mainly water-supply aquifers with larger permeability. There is no pumping-well distribution in the third and fifth layers because their hydraulic conductivity was lower than that of the other layers.
In the study area, the recharge sources of aquifer mainly include precipitation infiltration, lateral recharge of river and lateral groundwater flow recharge of the alluvial-diluvial fan under natural conditions. Groundwater flow discharges into the sea from east to west along the terrain. Water resource demand increases along with the development of industry and agriculture; the river basin has a high number of agriculture wells in the second layer and four water source areas in the fourth layer (Appendix A). In addition, there are three long-term water level monitoring wells in the study area. Kriging interpolation was used to obtain Figure 3a according to the water level data from monitoring wells and pumping wells in 2015 (Appendix A, Table A1). The groundwater overexploitation has become the most important outflow of groundwater since the pumping wells were built. Groundwater is mined in large quantities, the groundwater level is lower than the river level and the groundwater is replenished by the river for a long time. The supply amount of the river is limited, and a cone of depression is formed near the water source areas and irrigation area ( Figure 3). aquifer ranges from the second layer to the fifth layer. The second layer is mainly composed of sand gravel with thickness of about 10-25 m. Hydraulic conductivity is between 10 and 40 m/day, which gradually decreases from the valley to both sides. The third layer is mainly composed of sandy loam with thickness of about 0-8 m and hydraulic conductivity of about 5 m/day. The fourth layer is composed of gravel pebble, with thickness of about 0-15 m. The layer gradually becomes thinner from the riverbed to the hills, with hydraulic conductivity between 5 and 70 m/day. The fifth layer is mainly composed of clay containing some gravel. Hydraulic conductivity is about 4 m/day. The second and fourth layers are mainly water-supply aquifers with larger permeability. There is no pumping-well distribution in the third and fifth layers because their hydraulic conductivity was lower than that of the other layers.
In the study area, the recharge sources of aquifer mainly include precipitation infiltration, lateral recharge of river and lateral groundwater flow recharge of the alluvial-diluvial fan under natural conditions. Groundwater flow discharges into the sea from east to west along the terrain. Water resource demand increases along with the development of industry and agriculture; the river basin has a high number of agriculture wells in the second layer and four water source areas in the fourth layer (Appendix A). In addition, there are three long-term water level monitoring wells in the study area. Kriging interpolation was used to obtain Figure 3a according to the water level data from monitoring wells and pumping wells in 2015 (Appendix A, Table A1). The groundwater overexploitation has become the most important outflow of groundwater since the pumping wells were built. Groundwater is mined in large quantities, the groundwater level is lower than the river level and the groundwater is replenished by the river for a long time. The supply amount of the river is limited, and a cone of depression is formed near the water source areas and irrigation area ( Figure 3).

Overview of Water Resource Utilization
The original exploitation had no detailed record in water source areas, with single well production of possibly less than 2000 m 3 /day. The structure of irrigation planting changed from natural planting to greenhouses planting, and the pumping increase with the pumping mode changed from seasonal to continuous pumping. Groundwater exploitation causes a series of environmental groundwater problems, and seawater intrusion is the most prominent. The first survey on seawater intrusion was carried out in 1991, and seawater intrusion occurred in the lower reaches of the branch of the Daqing River. The exploitation of groundwater was not controlled in the study area, in reverse, while groundwater withdrawal increases with the demand of water resources. For the sustainable utilization of groundwater, the exploitation quantity in all the water source sites decreased annually in 2012-2015, and the total quantity of water source areas was less than 9 million m 3 /year in 2015. However, the exploitation quantity of irrigation wells was not controlled. The total pumping quantity of all wells in the area in 1991-2015 is listed in Table 1.

Overview of Water Resource Utilization
The original exploitation had no detailed record in water source areas, with single well production of possibly less than 2000 m 3 /day. The structure of irrigation planting changed from natural planting to greenhouses planting, and the pumping increase with the pumping mode changed from seasonal to continuous pumping. Groundwater exploitation causes a series of environmental groundwater problems, and seawater intrusion is the most prominent. The first survey on seawater intrusion was carried out in 1991, and seawater intrusion occurred in the lower reaches of the branch of the Daqing River. The exploitation of groundwater was not controlled in the study area, in reverse, while groundwater withdrawal increases with the demand of water resources. For the sustainable utilization of groundwater, the exploitation quantity in all the water source sites decreased annually in 2012-2015, and the total quantity of water source areas was less than 9 million m 3 /year in 2015. However, the exploitation quantity of irrigation wells was not controlled. The total pumping quantity of all wells in the area in 1991-2015 is listed in Table 1.

Hydrogeological Conceptual Model
Taking the Beijing 54 coordinate system as the simulation coordinate system, according to the hydrogeological characteristics of the study area, the simulation area was determined. The southern and northern sides are bounded by low mountains and hills. The eastern side is bounded by the reach where Wangbaoshan hydrological station is located, and the western side extends 2 km from the coast to Liaodong Bay ( Figure 1). The boundary conditions of the model were set as follows: The northern and southern sides of the model are constant flow boundaries that accept various lateral recharges. The eastern boundary is a constant head boundary based on the monitoring data of the hydrological station. The western boundary extending 2 km into sea can be defined as transient head and has constant chloridion concentration with H(t) and 20,000 mg/L (the monitoring data in this area are mainly chloride concentration), respectively. The fluctuating water head under the tidal action of the upper part of the silt coast is the third kind of boundary condition. The barrage can be regard as a constant head boundary. The regional reach is treated as a river in the model (water head and riverbed elevation are constant). The top boundary is the flow boundary, and the precipitation is the annual average variation measured by Yingkou meteorological station in 1991-2015. The buried depth of groundwater in the area was larger than the critical depth of groundwater loss (4 m). Therefore, evaporation is ignored. The bottom boundary was regarded as a no-flow boundary. The sources and sinks in the area mainly include groundwater-pumping wells. The exploitation amount is shown in Table 1. The initial conditions of the model in Appendix B ( Figure A1 and Table A3).
The stratigraphic deposition in the area is in accordance with the sedimentary law of the river valley coastal plain, and the overburden was approximately divided into five layers according to the actual stratigraphic structure (Section 1.2), with different permeability. In addition, due to the different parameters from the riverbed to the hills on both sides and from the hydrological station to coast, each layer was divided into eight parameter zones ( Figure 4). The aquifer is homogeneous and anisotropic in the same parameter zone. The initial hydraulic conductivity was assigned according to the results obtained from the pumping tests in the survey report [59]. The horizontal value of hydraulic conductivity was assumed to be 10 times its vertical value, and the rainfall infiltration coefficient and effective porosity are 0.22 and 0.24, respectively [60,61].

Numerical Model
In this paper, the mathematical model for variable density groundwater flow was applied to simulate the seawater intrusion process, and the distribution of seawater intrusion in the multilayered aquifers was analyzed. SEAWAT-2000 is numerical simulation software that combines MT3D with MODFLOW-2000 (https://www.aquaveo.com/), considering the effect of density on groundwater flow.

Mathematical Model
The governing equation for the variable density groundwater flow [62] in the study area can be expressed as Equation (1): where hf is the equivalent fresh water head (L); ρf is the density of fresh water (M/L 3 ); qs is unit volume flow of the source (sink) (1/T); ρ is the density of flow (M/L 3 ); ρs is the density of the source (sink) (M/L 3 ); θ is the effective porosity of porous medium; Sf is the unit water storage coefficient of equivalent freshwater (1/L); Kf is the hydraulic conductivity of equivalent freshwater (L/T). Equation (1) and the corresponding conditions of determining the solution constitute a mathematical model of groundwater flow. The conditions of determining the solution can be expressed as follows: Initial condition (Equation (2)):

Numerical Model
In this paper, the mathematical model for variable density groundwater flow was applied to simulate the seawater intrusion process, and the distribution of seawater intrusion in the multi-layered aquifers was analyzed. SEAWAT-2000 is numerical simulation software that combines MT3D with MODFLOW-2000 (https://www.aquaveo.com/), considering the effect of density on groundwater flow.

Mathematical Model
The governing equation for the variable density groundwater flow [62] in the study area can be expressed as Equation (1): where h f is the equivalent fresh water head (L); ρ f is the density of fresh water (M/L 3 ); q s is unit volume flow of the source (sink) (1/T); ρ is the density of flow (M/L 3 ); ρ s is the density of the source (sink) (M/L 3 ); θ is the effective porosity of porous medium; S f is the unit water storage coefficient of equivalent freshwater (1/L); K f is the hydraulic conductivity of equivalent freshwater (L/T). Equation (1) and the corresponding conditions of determining the solution constitute a mathematical model of groundwater flow. The conditions of determining the solution can be expressed as follows: Initial condition (Equation (2)): Sustainability 2020, 12, 2842 8 of 22 Dirichlet boundary (Equation (3)): Neumann boundary (Equation (4)): The third kind of boundary (Equation (5)): H 0 (x, y, z) is the initial head value of each layer, H 1 (x, y, z, t) is the head boundary, q(x, y, z, t) is the unit discharge on the flow boundary, K 1 is the hydraulic conductivity of the silt layer, m 1 is the thickness of the coast silt, H n is the head outside the boundary, H is the head inside the boundary.
The equation of solute transport [63] is expressed as in Equation (6): where C is the dissolved concentration of substance (M/L 3 ); D is the hydrodynamic dispersion coefficient (L 2 /T); C s is concentration of substance in source or sink term (M/L 3 ). Equation (6) and the conditions for determining the solution constitute a mathematical model of solute transport. The conditions of determining the solution can be expressed as follows: Initial concentration condition (Equation (7)): The first kind of concentration boundary (Equation (8)): where C 0 (x, y, z) is initial concentration distribution, and C 1 (x, y, z, t) is the measured concentration for the first concentration boundary.

Model Domain Discretization
The model was discretized into rectangular grids with 100 m × 100 m in the horizontal direction, and divided into five layers according to the corresponding aquifer characteristics in the vertical direction. There are 90,175 cells within the valid boundary. The simulation period is from November 1991 to November 2015, two years were taken as a stress period, and the output time step in the stress period was two months.

Model Calibration
Initial aquifer parameters (Appendix B, Table A4) were assigned according to the results of the surveyed report in 1973 [61]. It is necessary to use the observation data for parameter identification due to there being a deviation of simulation by using the initial parameters directly. There are three head observation wells and three concentration observation points (Appendix A, Table A2) in the study area ( Figure 1). The water level data of the pumping wells were measured in November 2015. the hydraulic conductivity, specific yield and specific storage were adjusted to match the calculated head and concentration with the observation values so as to meet the distribution law of flow and concentration field over time. The results of the identified model parameters are listed in Table 2; meanwhile, the horizontal hydraulic conductivity is 4.9 times the vertical value after identification. According to the existing simulation cases, the longitudinal dispersivity and transversal dispersivity are 500 m and 50 m, respectively [42], and the diffusion coefficient is 1 × 10 −4 m 2 /day.
In order to verify whether the identified parameters can truly reflect the hydrogeological characteristics of the area, the model validation period was selected as November 2015 to November 2016. During the validation period, the precipitation was assigned according to the monthly precipitation, and other parameters remained unchanged. One month was taken as a stress period and three days as an output time step. The calculated values of head were compared with the observed values of the head at observation wells (W1, W2, W3). The simulated concentrations were compared with the observed ones at wells W4, W5 and W6, as shown in Figure 5. Figure 5 shows that the variation trend is the same between the observed data the calculated data, and the average relative error values between the calculated values and the observed values are less than 5%. According to the head comparison curve between the observation well data and the calculated data, the parameters obtained by the model identification are basically in line with the actual situation. The model can be further used to study the seawater intrusion laws in the study area.

Model Uncertainty Analysis
The aquifer division in the model was determined according to the structure and distribution of the actual aquifer. Since the thicknesses of the aquifer on both sides of the river valley are reduced from top layer to bottom layer, in order to facilitate the simulation calculation, it was assumed that the distribution range of each layer is the same, and the part of each layer beyond the distribution range of the actual aquifer was given as a smaller value of hydraulic conductivity. The error caused

Model Uncertainty Analysis
The aquifer division in the model was determined according to the structure and distribution of the actual aquifer. Since the thicknesses of the aquifer on both sides of the river valley are reduced from top layer to bottom layer, in order to facilitate the simulation calculation, it was assumed that the distribution range of each layer is the same, and the part of each layer beyond the distribution range of the actual aquifer was given as a smaller value of hydraulic conductivity. The error caused by this method can be ignored through multiple sets of simulation comparison.
The uncertainty of the precipitation value was analyzed (Appendix C). The error was analyzed by comparing the calculation results of the annual average and monthly precipitation at the same output time. The result (Appendix C, Figure A2) shows that the error between the two cases is small and can be ignored. Therefore, the average value of precipitation can be assigned in the model to study the laws of seawater intrusion.
Sensitivity analysis is conducted for the extended distance to the sea. The intertidal coast is flat, and the silt layer extends 2 km to the sea in the model. The tidal head boundary was set on the silt layer. For sensitivity analysis, it was assumed that the extension distance of the silt layer to the sea is 1 km and 3 km, respectively. The simulation results show that the extended distance has a tiny effect on the concentration results under the effect of the tide.

Sobol Global Sensitivity Analysis Method
In the process of seawater intrusion, the external parameters such as local precipitation, pumping quantity, sea level and other related parameters such as permeability parameters and dispersion have an important impact on seawater intrusion. Laboratory experiments have shown that the differences in permeability parameters of different layers will cause different seawater intrusion in the multi-layered aquifers [21]. However, the reason for the differences in seawater intrusion between every two layers is not only the permeability parameters in the field. As for the role of the factors in layered intrusion, Sobol global sensitivity analysis [48] is carried out in the first layer, second layer and fourth layer. The water exchanges between Daqing River and Bohai Sea in the first layer, which is influenced by average sea level. The barrage has the effect of blocking the tide and can be used as a barrier between saltwater and freshwater. Therefore, the freshwater level at the barrage is very important. In addition, this layer is directly supplied by atmospheric precipitation, which is closely related to the first layer. Therefore, the sensitivity analysis selected parameters 1 (the water level at barrage), 2 (the mean sea level) and 3 (precipitation). A large number of irrigation wells in the second layer receive the saltwater supply directly from the sea and the first layer. Therefore, parameters 2, 3 and 4 (Q1-single irrigation well pumping quantity) were selected for sensitivity analysis. The fourth and second layers are separated by the third layer, which is a relatively weak permeable layer. The main feature of this layer is that there are three water sources for large-scale pumping groundwater. Therefore, the sensitivity analysis was carried out by selecting parameters 5 (Q2-single well pumping quantity of Yong'an water source area), 6 (Q3-single well pumping quantity of Yinzhuhuafang water source area) and 7 (Q4-single well pumping quantity of Gaizhou 2 and 3 water source area).
Before the sensitivity analysis, the annual precipitation and the maximal pumping quantity of the single well were taken as the mean values, and the parameter distribution ranged from 0.5 to 1.5 times the mean values. In addition, the average sea level and the water level at the barrage were given as a fluctuation range. Mean values and the distribution range of the specific setting are listed in Table 3. Then, the Monte Carlo method is used to sample the distribution range to obtain the combination of parameter sensitivity analysis.

Results
After the model calibration and uncertainty analysis, the real hydrogeological characteristics of the area can be represented, and the simulation results can reflect the law of seawater intrusion in the area. The quantity of pumping reached the maximal value in the area during 2006-2011. The maximal quantity of pumping in the water source areas and irrigation area is about 55 million m 3 /year. The rainfall supply, lateral supply and river flow are about 17 million m 3 /year, 4 million m 3 /year and 25 million m 3 /year, respectively. The total freshwater flow (including river flow, rainfall supply and groundwater) is about 46 million m 3 /year in the global study area, which is less than the total quantity of pumping and belongs to over-pumping. The pumping-limit scenario was implemented to control the quantity of pumping at each water source area in 2012. From 2012 to 2013, the quantity of pumping was reduced to half of the original quantity at Yong'an and Tuandian water source areas. The quantity of single well pumping was less than 1000 m 3 /day at each water source area after the pressure mining was again carried out in 2014. However, the pumping of irrigation wells is not controlled. The area with concentration greater than 250 mg/L is regarded as the seawater intrusion area. The seawater intrusion area of the study area is defined as A 0 , and the corresponding seawater intrusion area is defined as A t at time t under no groundwater-pumping conditions. Therefore, the seawater intrusion area at corresponding any time node is ∆A t = A t − A 0 . The changes of seawater intrusion in the multi-layered aquifers were calculated by simulation after the pumping-limit ( Table 4). The seawater intrusion in the multi-layered aquifers meets the law: The intrusion area is from the first layer to the fourth layer, and the fifth layer is approximately equal to the fourth layer. While the intrusion area increases after the pumping-limit scenario, the intrusion rate decreases in the multi-layered aquifers, except the first layer, which is affected by the tide. The simulation results in November 2015 were used to analyze the laws of seawater intrusion in the area (Figures 6 and 7). Due to over-pumping in the area for many years, a depression cone is formed near the Yong'an water source area in the south of Daqing River, where the pumping wells are relatively concentrated (Figure 6). The lowest head in the center of the depression cone was historically about -16 m. Figure 7 shows the concentration distribution of each layer in the study area. Due to the differences in density, hydrogeological parameters and pumping quantity, the concentration isoclines show differences in different layers. The reasons for the differences in layers mainly include the three following aspects:  (1) The seawater intrusion was influenced by tide and regulation of barrage. Figure 7a shows that the saltwater intrusion is mainly controlled by the tide and barrage in the first layer. In the Daqing River channel near the estuary, the high tide level seawater moves upward along the channel, and there is a high concentration of seawater intrusion front in the channel. In the vicinity of the barrage, there is an intrusion line inclined to the sea in a small range, which is controlled by the high water level of the barrage.  (1) The seawater intrusion was influenced by tide and regulation of barrage. Figure 7a shows that the saltwater intrusion is mainly controlled by the tide and barrage in the first layer. In the Daqing River channel near the estuary, the high tide level seawater moves upward along the channel, and there is a high concentration of seawater intrusion front in the channel. In the vicinity of the barrage, there is an intrusion line inclined to the sea in a small range, which is controlled by the high  (1) The seawater intrusion was influenced by tide and regulation of barrage. Figure 7a shows that the saltwater intrusion is mainly controlled by the tide and barrage in the first layer. In the Daqing River channel near the estuary, the high tide level seawater moves upward along the channel, and there is a high concentration of seawater intrusion front in the channel. In the vicinity of the barrage, there is an intrusion line inclined to the sea in a small range, which is controlled by the high water level of the barrage.
(2) The seawater intrusion was influenced by different quantities of pumping in the different layers. The pumping wells are mainly distributed in the second and fourth layers with high permeability. The wells in the second layer are pumped for irrigation. The groundwater level decreases due to the large number of irrigation wells being pumped. The high sea level and the saltwater entering the river channel under the action of the tide together recharge the groundwater. The intrusion area increases with the increase of exploitation (Figure 7b). The fourth layer is dominated by the pumping wells in the water source areas. A depression cone was formed due to the long period of pumping in the wells in the water source area. The direction of local hydraulic gradient changes from pointing to the sea to pointing to the center of the depression cone. Saltwater intrudes into this layer with the largest intrusion area (Figure 7d).
(3) Seawater intrusion was influenced by different permeability parameters in the different layers. Because the hydraulic conductivity of the first layer is less than that of the second layer, when a large number of irrigation wells pumping groundwater to receive a high level of saltwater supply the river channel, the saltwater body enters the first layer and then infiltrates. The second layer also receives the saltwater supply directly from the sea, so the intrusion area of the first layer is smaller than that of the second layer. Guo et al. [21] showed that the clay lens in the aquifer is bypassed by the saltwater body in the process of seawater intrusion, and is slowly diffused under the effect of gravity. Therefore, the third layer is affected by the second and fourth layers, and the intrusion area is between the second and fourth layers. The fifth layer is mainly affected by the saltwater intrusion in the fourth layer, and there is almost no difference with the concentration distribution of the fourth layer.
Therefore, the main reason for the seawater intrusion in the multi-layered aquifers in this area is that the amount of layered pumping is different. The intrusion area is the largest in the fourth layer, and the intrusion area increases with the increase of irrigation pumping in the second layer. The difference in permeability parameters of each layer resulted in the larger intrusion area in the second layer. Under the action of gravity, the intrusion area in the fifth layer was approximately equal to the fourth layer, and the intrusion area in third layer was between the second and fourth layers. In addition, the tide and the barrage work together to form the unique intrusion phenomenon in the first layer. This section may be divided by subheadings.

Reason Analysis of Seawater Intrusion in Multi-Layered Aquifers
Sobol global sensitivity analysis [48] was carried out in the first, second and fourth layers. In the process of sensitivity analysis, six observation points O1-O6 ( Figure 1) were evenly selected in the intruded area, and the chloridion concentration values at different observation points were taken as the model output results. When the Sobol method was used for sensitivity analysis, the seawater intrusion model in the study area was operated for one year. According to the chloridion concentration of six observation points, the full-order sensitivity coefficient S Ti (Figures 8-10) of each parameter at different locations could be calculated [50]. Figure 8 is the calculation results of the full order sensitivity coefficient at each point in the first layer. The water level of the barrage further affected the concentration near the river channel by affecting the downstream river level. The O4 is most sensitive where closest to the river channel, and the sensitivity of the other points decreases with an increase in the distance from the river. The average sea level has an impact on each point, and the O3 surrounded by the river is less sensitive to the sea level.
process of sensitivity analysis, six observation points O1-O6 (Figure 1) were evenly selected in the intruded area, and the chloridion concentration values at different observation points were taken as the model output results. When the Sobol method was used for sensitivity analysis, the seawater intrusion model in the study area was operated for one year. According to the chloridion concentration of six observation points, the full-order sensitivity coefficient STi (Figures 8-10) of each parameter at different locations could be calculated [50].     Figure 9. The full-order sensitivity coefficient of the specific points in the second layer. Figure 9 is the calculation results of the full order sensitivity coefficient at each point in the second layer. Concentration is sensitive to the sea level and quantity of pumping at each point, and O4 is the most sensitive to the sea level in the lower part of the riverbed. The number of small-scale irrigation wells in the north plain of the river in the area is small, and the O4-O6 on the north are less sensitive to the pumping quantity of the irrigation wells. The valley plain on the south of the river is relatively flat and wide, the number of the irrigation wells is larger and there is river between the two banks, which, to some extent, interferes with the hydraulic connection. Therefore, the O1-O3 on the south are sensitive to the response of pumping quantity. Figure 10 is the calculation results of the full-order sensitivity coefficient at each point in the fourth layer. The concentration of each point is sensitive to the changes of the pumping quantity of the Yong'an water source, but almost insensitive to the other water sources. This may be because other water sources are located in the upstream far from the observation point with fewer pumping wells, and a watershed formed by pumping between the Yong'an water source area and others.   Figure 8 is the calculation results of the full order sensitivity coefficient at each point in the first layer. The water level of the barrage further affected the concentration near the river channel by affecting the downstream river level. The O4 is most sensitive where closest to the river channel, and the sensitivity of the other points decreases with an increase in the distance from the river. The average sea level has an impact on each point, and the O3 surrounded by the river is less sensitive to the sea level. Figure 9 is the calculation results of the full order sensitivity coefficient at each point in the second layer. Concentration is sensitive to the sea level and quantity of pumping at each point, and O4 is the most sensitive to the sea level in the lower part of the riverbed. The number of small-scale irrigation wells in the north plain of the river in the area is small, and the O4-O6 on the north are less sensitive to the pumping quantity of the irrigation wells. The valley plain on the south of the river is relatively flat and wide, the number of the irrigation wells is larger and there is river between the two banks, which, to some extent, interferes with the hydraulic connection. Therefore, the O1-O3 on the south are sensitive to the response of pumping quantity. In general, the sensitivity analysis results showed that the difference of the pumping quantity in different layers is the most important reason for the difference of the intrusion area in the second and fourth layers. The intrusion phenomenon in the first layer was formed by the influences of tide and barrage.

Prediction of Seawater Intrusion in the Multi-Layered Aquifers
Pumping was limited after 2015; since then, the total pumping quantity at the water source areas has not been more than 9 million m 3 /a. The total irrigation pumping quantity is about 25 million m 3 /a without pumping-limit. The quantity of pumping is less than the total amount of freshwater resources flowing through the area, and the quantity of pumping in the second layer is greater than in the fourth layer. According to the results of sensitivity analysis, the pumping quantity is the key factor affecting the intrusion area in the second and fourth layers. Therefore, in predictions for the next two decades, the intrusion area for different scenarios was mainly discussed under the condition of the same quantity of pumping. The specific scenario was designed as follows: In Scenario I, the quantity of pumping is 25 million m 3 /a in the second layer, which undertakes all irrigation tasks, that is, 9 million m 3 /a in the fourth layer. In Scenario II, the pumping quantity in both the second and fourth layers is 17 million m 3 /a, and part of the pumping quantity at the water source areas is used as irrigation. The specific simulation settings are as follows: The prediction simulation period was 20 years, from November 2015 to November 2035, the hydrogeological parameters remained unchanged, the precipitation took the annual average value and the evaporation was ignored. The groundwater flow and chloride concentration field in November 2015 were selected as the initial flow and concentration field, and the pumping quantity was set according to Scenarios I and II, respectively.
The simulated seawater intrusion areas in the multi-layered aquifers are listed in Table 5. The seawater intrusion area for Scenario I is larger in first to third layers, compared with that of Scenario II. However, the intrusion area was approximately equal in the fourth and fifth layers for both scenarios. That is to say, the pumping quantity at water source areas is less than that of irrigation wells after a pumping-limit, and the large-scale pumping of irrigation wells is the main reason for the increase of the intrusion area. However, the saltwater began to retreat after the recovery of the groundwater head surface at water source areas in the fourth layer. If the pumping wells in the water source areas bear part of the irrigation demand within their pumping capacity, so as to reduce the pressure of pumping in the second layer, the overall seawater intrusion area can be reduced by about 0.23 km 2 under the same pumping quantity. Therefore, in order to meet the production demand of the area and ensure that the intrusion area can be annually reduced, the pumping wells at water source areas can be used to undertake part of the irrigation demand.

Conclusions
On the basis of field investigation and long-term monitoring data analysis, a three-dimensional numerical model was built in the estuary of Daqing River in Yingkou City, Liaoning Bay, China. Based on the observation data, the model was calibrated so that it could reflect the hydrogeological characteristics of the study area. The laws of seawater intrusion in the layered aquifer under the layered pumping were calculated and analyzed by SEAWAT-2000, and the sensitivity was analyzed with the Sobol method, and the change of seawater intrusion in the multi-layered aquifers in the future was predicted. According to the study results, the following conclusions were drawn: (1) The seawater intrusion area in different layers is different, with an obvious layered law. There was a direct relationship between the layered intrusion area and the quantity of layered pumping. Due to the long-term over-pumping of groundwater, the value of the intrusion area in the fourth layer was the largest, while the intrusion area increased with the pumping quantity in the second layer.
(2) Difference in permeability and other factors further affects the layered intrusion. Under the effects of tide, barrage and different permeability, the seawater intrusion phenomenon is special in first layer and the intrusion area is less than that of the second layer. The saltwater in the fourth layer infiltrates into the fifth layer under the action of gravity, and its intrusion area is approximately equal to the fourth layer.
(3) The depression cone recovered after a pumping-limit in the three centralized water source areas; then, the intrusion area started to retreat in the fourth layer. Simultaneously, the pumping quantity of the irrigation wells became the main reason for the increases of the intrusion area. If the wells at water source sites are used to bear part of the irrigation demand, so as to reduce the pressure of pumping in the second layer, the overall intrusion area can be reduced by about 0.23km 2 under the same pumping quantity. Therefore, in order to meet the production demand of the area and ensure that the intrusion area annually decreases, the wells in the water source site are used to undertake part of the irrigation demand.
Beyond that, the model has some limitations. The barrage was treated as a constant head, but inflow from the upper reaches is sometimes not enough, and the water level at the barrage is lower than that at the constant head. Therefore, the predicted extent of seawater intrusion downstream from the barrage was less than that of the actual intrusion.

Appendix A
Water resources demand increase along with the development of industry and agriculture. The study area has a high number of irrigation wells and four water source areas. The four water source areas are Tuandian (built in 1965; total number of wells, 19; four wells in the study area), Yinzhuhuafang (built in 1970; eight wells), Gaizhou 2 and 3 (built in 1970; 13 wells), Yong'an (built in 1936; 18 wells) from upstream to downstream (Figure 1). All irrigation wells were arranged in the farming area, about 1.5 wells/km 2 , with a total of 90 wells. Local water level measurement and statistics were carried out in November 2015 (Table A1). Statistics included the part of the pumping wells. In addition, there are three water level monitoring wells and three concentration monitoring wells. Details of the monitoring wells are listed in Table A2.

Appendix B
The initial steady-state flow field was calculated according to the pumping of the water source areas and irrigation wells in 1991 ( Figure A1). The observed concentration values (Table A3) in 1991 are the initial concentrations in the study area. The initial parameters of the model are listed in Table A4.

Appendix B
The initial steady-state flow field was calculated according to the pumping of the water source areas and irrigation wells in 1991 ( Figure A1). The observed concentration values (Table A3) in 1991 are the initial concentrations in the study area. The initial parameters of the model are listed in Table  A4.   Figure A1. The flow field in 1991.

Appendix C
The annual average precipitation was directly assigned in the model because of the long simulation period. In order to analyze the impact of this on the results, monthly average and annual average precipitation were, separately, assigned in the model from 2011 to 2015. We compared the calculated results of the same output time at four different points ( Figure A2).