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

Abstract

Groundwater over-pumping in estuary cities leads to a series of groundwater environmental problems that seriously restricts economic development. On the basis of field 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 five 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 affected the layered intrusion area, as different permeability, tide and barrage further affected 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² under the same pumping quantity.

1. 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].

Seawater intrusion is a density-dependent flow problem [4,28,29]. Several methods have been developed to address it, including analytical solution, experimental and simulation methods [18,21,23,24,28,29,30,31,32]. Analytical solution is always used to locate a steady fresh–saltwater interface location under many assumed conditions [19,33,34,35]. Laboratory experiment is limited to scale effect, and requires higher testing operations and monitoring means [23,28,30,36]. The numerical method is a powerful tool for modeling and solving the problem of seawater intrusion because it has some advantages that can depict the migration law of the freshwater and saltwater interface under many complicated hydrogeological conditions and artificial factors [4,8,17,29,37,38,39]. The SUTRA [40], FEFLOW [1], VDFST-CFP [41], SEAWAT [42,43] and other types of numerical software were developed and are widely used in the simulations of seawater intrusion.

Complex hydrological models are controlled by a high number of parameters [44]. Sensitivity analysis can determine the most influential parameters of the model [45]; it includes local sensitivity and global sensitivity. Specific sensitivity analysis methods include Fourier amplitude sensitivity [46], the Morris method [47], the Sobol method [48], etc. The Sobol method, through variance decomposition of the model results, quantitatively obtains the sensitivity of each parameter, and has been widely used in hydrological models in recent years [44,49,50].

Groundwater overexploitation is one of the key factors in inducing seawater intrusion, as has been documented in many studies documents over the years [18,24,25,26,51,52]. Most of the previous researchers consider seawater intrusion distribution in single layer aquifers [1,18,20,23]. A coastal stratum deposit is affected by the alluvial–diluvial effects of coastal rivers and presents the stratification phenomenon. Some scholars have studied the law of seawater intrusion in the multi-layered aquifer system [7,19,21,42,53,54,55,56,57]. A few studies have considered seawater intrusion under the influence of pumping in multi-layered aquifers. For example, Guo et al. [21] studied saltwater intrusion in multi-layered aquifers in the laboratory, but only pumping in a confined aquifer. Chang et al. [42] studied the interface change between brine and freshwater under pumping in multi-layered aquifers, but studies regard the complex ground-well pumping as a two-dimensional cross-section. At present, the study of layered saltwater intrusion caused by well group pumping in multi-layered aquifers in the field is not enough.

In this study, a three-dimensional numerical model was built in the estuary of the Daqing River in Liaodong Bay, China, on the basis of field investigation and long-term monitoring data analysis. The seawater intrusion distribution affected by pumping in the multi-layered aquifers was calculated using SEAWAT-2000. Then, the extent and area of seawater intrusion in the multi-layered aquifers was estimated. The reasons for seawater intrusion in the multi-layered aquifers were analyzed based on the method of Sobol sensitivity analysis. Finally, the change trend of seawater intrusion in the multi-layered aquifers was predicted for the coming two decades, under the same pumping quantity and different pumping scenarios, considering the pumping-limit in the centralized water sources. The study results can provide references for the sustainable utilization of groundwater in the local area.

1.1. Study Area

The study area is located in the eastern coastal area of Liaodong Bay, which is in Yingkou City, Liaoning province. It is downstream from the hydrological station of the Daqing river basin, as shown in Figure 1. Geographically, the area range is between latitude 40°20′ and 40°26′ and longitude ranging from 122°10′–122°35′. The total area is about 170 km², with a length of about 24 km from east to west, and a length between 2 km and 13 km from north to south. The coastline length of in the study area is about 6 km. The climate of the study area is continental monsoon type, with four distinct seasons and the same season of rain and heat [58]. The average annual precipitation was 600–800 mm from 1956 to 2016 [58,59]. The average annual evaporation was high, ranging from 1000 mm to 1200 mm because the average amount of annual sunshine hours can be between 2600 h and 2880 h. The maximum evaporation depth was 4 m. The average temperature ranged from 9 °C to 10 °C for many years, and the frost-free period is 180–210 days in every year.

1.2. 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.

1.3. Hydrogeological Conditions

1.3.1. 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². 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³/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³/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.

1.3.2. 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 water level data of the pumping wells and monitoring wells in 2015.

WELL	East Longitude (◦)	North Latitude (◦)	Head (m)	WELL	East Longitude (◦)	North Latitude (◦)	Head (m)
Yong’an-1	122.325	40.401	5	Yinzhuhuafang-1	122.398	40.364	4.9
Yong’an-2	122.322	40.398	−6.94	Yinzhuhuafang-2	122.396	40.367	4.94
Yong’an-3	122.327	40.394	−10.2	Yinzhuhuafang-3	122.394	40.372	5.13
Yong’an-4	122.325	40.393	-15	Yinzhuhuafang-4	122.392	40.368	5.1
Yong’an-5	122.326	40.388	−13.6	Yinzhuhuafang-5	122.383	40.373	4.89
Yong’an-6	122.315	40.393	−13	Gaizhou 2&3-1	122.356	40.383	4.98
Yong’an-7	122.325	40.384	−13.4	Gaizhou 2&3-2	122.357	40.380	4.98
Yong’an-8	122.318	40.389	−13.21	Gaizhou 2&3-3	122.359	40.377	5.04
Yong’an-9	122.328	40.296	−11	Gaizhou 2&3-4	122.365	40.374	5.03
Tuandian-1	122.456	40.426	22.49	Gaizhou 2&3-5	122.352	40.373	4.89
Tuandian-2	122.452	40.423	22.51	W1	122.257	40.358	8.75
Tuandian-3	122.452	40.422	22.53	W2	122.393	40.357	6.7
Tuandian-4	122.451	40.420	22.49	W3	122.290	40.389	1.63

). 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).

1.3.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³/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³/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. Groundwater exploitation in the study area.

Time	Single Well Pump-Out (m3/day)
	Field Irrigation	Tuandian Water Source	Gaizhou 2 and 3 Water Source	Yinzhuhuafang Water Source	Yongan Water Source
1991–1995	400	800	1500	1500	2000
1996–2000	400	800	1500	2000	2000
2001–2005	500	1000	1800	2000	2300
2006–2010	600	1300	1800	2400	3000
2011–2015	600–700	400–1300	800–1800	200–2400	500–3000

.

2. Materials and Methods

2.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 observed concentration values in 1991.

Point Number	East Longitude (◦)	North Latitude (◦)	Depth (m)		Concentration (mg/L)
1	122.258	40.412	25	270
2	122.278	40.435	5	382
3	122.322	40.451	10	263
4	122.307	40.441	12	205
5	122.292	40.459	4	1184
6	122.293	40.421	8	107
7	122.274	40.416	30	165
8	122.284	40.468	6	18701

).

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].

2.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.

2.2.1. Mathematical Model

The governing equation for the variable density groundwater flow [62] in the study area can be expressed as Equation (1):

$$\begin{array}{l}\frac{\partial }{\partial x}\times \left[\rho K_{fx}\times \left(\frac{\partial h_f}{\partial x}+\frac{\rho -{\rho }_f}{{\rho }_f}\times \frac{\partial z}{\partial x}\right)\right]+\frac{\partial }{\partial y}\times \left[\rho K_{fy}\times \left(\frac{\partial h_f}{\partial y}+\frac{\rho -{\rho }_f}{{\rho }_f}\frac{\partial z}{\partial y}\right)\right]\\ +\frac{\partial }{\partial z}\times \left[\rho K_{fz}\times \left(\frac{\partial h_f}{\partial z}+\frac{\rho -{\rho }_f}{{\rho }_f}\frac{\partial Hz}{\partial x}\right)\right]=\rho S_f\times \frac{\partial h_f}{\partial t}+\theta \times \frac{\partial \rho }{\partial c}\frac{\partial c}{\partial t}-{\rho }_s\times q_s\end{array}$$

(1)

where h_f is the equivalent fresh water head (L); ρ_f is the density of fresh water (M/L³); q_s is unit volume flow of the source (sink) (1/T); ρ is the density of flow (M/L³); ρ_s is the density of the source (sink) (M/L³); θ 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)):

$$H\left(x,y,z,0\right)=H_0(x,y,z)$$

(2)

Dirichlet boundary (Equation (3)):

$$H(x,y,z,t)|{}_{{\Gamma }_1}=H_1(x,y,z,t)$$

(3)

Neumann boundary (Equation (4)):

$$K\frac{\partial H(x,y,z,t)}{\partial n}|{}_{{\Gamma }_2}=q(x,y,z,t)$$

(4)

The third kind of boundary (Equation (5)):

$$K\times \frac{\partial H}{\partial n}|{}_{{\Gamma }_3}=\frac{K_1}{m_1}\times (H_n-H)$$

(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₁ is the hydraulic conductivity of the silt layer, m₁ 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):

$$\frac{\partial C}{\partial t}=\nabla \left(D\times \nabla C\right)-\nabla \left(\overrightarrow{v}C\right)-\frac{q_s}{\theta }{C}_S$$

(6)

where C is the dissolved concentration of substance (M/L³); D is the hydrodynamic dispersion coefficient (L²/T); C_s is concentration of substance in source or sink term (M/L³).

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)):

$$C(x,y,z,0)=C_0(x,y,z)$$

(7)

The first kind of concentration boundary (Equation (8)):

$$C(x,y,z,t)|{}_{{\Gamma }_1}=C_1(x,y,z,t)$$

(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.

2.2.2. 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.

2.2.3. Model Calibration

Initial aquifer parameters (Appendix B,

Table A4. The initial parameters in the model.

Partition	Kh/m·d−1					Sy	Ss/m−1
	layer1	layer2	layer3	layer4	layer5	layer1	layer2	layer3	layer4	layer5
1	6	11	5	5	4	0.11	8 × 10−4	1 × 10−3	1 × 10−3	1 × 10−3
2	10	22	10	50	4	0.12	6 × 10−4	7 × 10−4	2 × 10−4	1 × 10−3
3	12	15	5	20	4	0.12	7 × 10−4	8 × 10−4	6 × 10−4	1 × 10−3
4	10	20	5	5	4	0.12	6 × 10−4	1 × 10−3	1 × 10−3	1 × 10−3
5	12	40	5	70	4	0.12	2×10−4	8 × 10−4	5 × 10−5	1 × 10−3
6	10	35	5	50	4	0.12	4 × 10−4	8 × 10−4	1.0 × 10−4	1 × 10−3
7	5	25	5	35	4	0.10	5 × 10−4	8 × 10−4	3 × 10−4	1 × 10−3
8	0.4	20	5	30	4	0.05	7 × 10−4	1 × 10−3	4 × 10−4	1 × 10−3

) 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. The detail of the monitoring wells.

WELL	East Longitude (◦)	North Latitude (◦)	Depth(m)	Screen Depth (m)	Observation Type
W1	122.257	40.358	25	15–25	Head
W2	122.393	40.357	40	25–40	Head
W3	122.290	40.389	20	10–20	Head
W4	122.259	40.393	28	10–28	Concentration
W5	122.292	40.401	45	18–45	Concentration
W6	122.324	40.446	25	10–25	Concentration

) 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. Model parameters.

Partition	Kh/m·day−1					Sy	Ss/m−1
	layer1	layer2	layer3	layer4	layer5	layer1	layer2	layer3	layer4	layer5
1	5.23	10.35	4.83	4.32	3.98	0.11	7.9 × 10−4	1.0 × 10−3	1.0 × 10−3	1.1 × 10−3
2	9.86	21.15	9.47	39.46	4.11	0.12	5.8 × 10−4	7.3 × 10−4	1.8 × 10−4	1.0 × 10−3
3	10.12	14.36	4.98	19.27	4.03	0.12	7.1 × 10−4	8.3 × 10−4	5.9 × 10−4	1.0 × 10−3
4	9.23	19.22	5.17	4.98	4.13	0.12	6.2 × 10−4	1.0 × 10−3	1.0 × 10−3	1.0 × 10−3
5	10.82	41.23	5.98	62.13	4.92	0.12	2.0 × 10−4	8.1 × 10−4	4.7 × 10−5	1.0 × 10−3
6	9.37	32.04	5.88	46.17	4.89	0.12	3.9 × 10−4	7.9 × 10−4	1.0 × 10−4	1.0 × 10−3
7	4.99	24.76	5.95	34.64	4.06	0.10	4.8 × 10−4	8.2 × 10−4	2.9 × 10−4	1.0 × 10−3
8	0.52	14.52	4.33	26.07	4.02	0.05	7.2 × 10−4	1.0 × 10−3	4.4 × 10−4	1.0 × 10−3

; 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⁻⁴ m²/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.

2.3. 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.

2.4. 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. The parameters for sensitivity analysis.

Parameter	Definition	Mean	Distribution Range
1	The water level at barrage (m)	2	(1.5,2.4)
2	mean sea level (m)	0.4	(−0.2,1.0)
3	precipitation (mm/a)	620	(310,930)
4	Q1 (m3/day)	700	(350,1050)
5	Q2 (m3/day)	3000	(1500,4500)
6	Q3 (m3/day)	2400	(1200,3600)
7	Q4 (m3/day)	1800	(900,2700)

. Then, the Monte Carlo method is used to sample the distribution range to obtain the combination of parameter sensitivity analysis.

3. Results and Discussion

3.1. 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³/year. The rainfall supply, lateral supply and river flow are about 17 million m³/year, 4 million m³/year and 25 million m³/year, respectively. The total freshwater flow (including river flow, rainfall supply and groundwater) is about 46 million m³/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³/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₀, 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 $\Delta 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 layered intrusion area changes year by year after pumping-limit.

Time	Layered Intrusion Area (km2)
	1	2	3	4	5
2011	16.865	17.785	18.353	18.566	18.566
2012	16.940	18.033	18.550	18.777	18.777
2013	17.170	18.205	18.765	18.956	18.956
2014	17.410	18.325	18.890	19.166	19.166
2015	17.557	18.514	19.037	19.281	19.281

). 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 (Figure 6 and Figure 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.

(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.

3.2. 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 (Figure 8, Figure 9 and Figure 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.

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.

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.

3.3. 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³/a. The total irrigation pumping quantity is about 25 million m³/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³/a in the second layer, which undertakes all irrigation tasks, that is, 9 million m³/a in the fourth layer. In Scenario II, the pumping quantity in both the second and fourth layers is 17 million m³/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 layered intrusion area after 20 years in different scenario.

Scenario	Layered Intrusion Area (km2)
	layer1	layer2	layer3	layer4	layer5
Scenario I	17.875	18.569	19.033	19.226	19.226
Scenario II	17.655	18.405	18.943	19.227	19.227

. 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² 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.

4. 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² 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.