## Abstract

Sea-level rise due to global warming and over-exploitation of groundwater resources in coastal areas will induce seawater intrusion into inland groundwater which is leading to the migration of the transition zone, and is affecting the security of regional social economy and water resources. Taking the blue economic zone of Shandong Peninsula as the research background, selecting the coastal area of Longkou and Zhaoyuan as the study area, this paper firstly depicts the groundwater flow field pattern in the research area, and then carries out the quantitative analysis for the migration patterns of the transition zone under changing groundwater levels. The results demonstrate that the width of sea-land transitional zone in the research area is about 1.5–4.5 km. The groundwater level is higher, the migration speed of the sea-land transitional zone is smaller. The results are of significance to study the migration of the sea-land transitional zone in the blue economic zone of Shandong peninsula.

## INTRODUCTION

With the growth of population in coastal areas, demand for freshwater has increased. To meet the ever-increasing needs of freshwater, many coastal areas have exploited groundwater on a large scale and this has resulted in a sharp decline of the inland groundwater level (Shi & Jiao 2014; Tam *et al.* 2014). On the other hand, average global sea levels are rising continually due to global warming. Under these two stressful conditions, groundwater over-exploitation and sea-level rise, the natural dynamic balance between the seawater and the inland groundwater has been broken in many coastal regions (Werner *et al.* 2013), causing the freshwater–seawater interface to migrate inland (Xue *et al.* 1992, 1993; Wu *et al.* 1996). As a result, seawater intrusion (SWI) has occurred in many countries and regions today (Guo & Huang 2003; Anders *et al.* 2014; Cary *et al.* 2015). Moreover, SWI has a huge negative impact on local social development, which can pollute the inland fresh groundwater resources and destroy the groundwater environment (Zeng *et al.* 2016).

Since SWI has become a common issue around the world, many scholars have conducted in-depth research on SWI, mainly including the concept, influencing factors, monitoring methods, prevention measures and so on (Guo & Huang 2003; Salem & Osman 2017). Especially, the shape, migration mechanism and migration law of the freshwater-seawater interface are core issues (Guo & Huang 2003). In the early stage, researchers saw the freshwater–seawater interface as a sharp interface and proposed many formulas to determine the interface according to the phreatic water level under the condition of hydrostatic equilibrium (Tang *et al.* 2008; Zhou 2008). However, since seawater and freshwater are easily miscible in coastal aquifers, the actual freshwater–seawater interface is a transition zone, the shape and the thickness of which are mainly affected by geological structure, hydrodynamic characteristics, tidal fluctuation and groundwater exploitation. Because of the complexity of the seawater groundwater transition zone, the transition zone is commonly depicted by two partial differential equations. One equation is used to describe the movement of the mixed fluid of seawater and freshwater under the changing density, and the other equation is applied to describe the migration of the solute in the mixed liquid (Xue *et al.* 1992, 1993). These two equations effectively realize the coupling of the variable density, solute concentration and groundwater levels.

The transition zone can only be calculated correctly by numerical methods (Guo & Huang 2003). Pinder & Cooper (1970) were the first scholars to study the transitional zone model of SWI in coastal areas, and presented a numerical solution by the finite element method for Henry's problem. Several works have been published to understand the extent of transition zone using finite element methods which include a two-dimensional finite element model (Segol *et al.* 1975), a finite difference method (Gupta & Yapa 1982) and a three-dimensional finite element model (Huyakorn *et al.* 1987; Xue *et al.* 1992, 1993).

Nowadays, with the rapid development of computer science, much numerical simulation software, which can make SWI study more convenient (Praveena *et al.* 2011), has been used to assess SWI. Numerical simulation software applied by scholars to successfully study SWI mainly include visual MODFLOW (Hu *et al.* 2015) and FEFLOW (Lu *et al.* 2010). Furthermore, to guarantee the accuracy of simulation results, deficiencies of related data, boundary conditions and stratum structure, which are the most common problems in building the numerical simulation model to study SWI (Post 2005; Werner *et al.* 2013), have been improved by many different effective methods (Siarkos & Latinopoulos 2016).

The main purpose of this paper is to study the condition of SWI in the coastal area of Longkou and Zhaoyuan in Shandong province, China, which is seriously affected by SWI, based on the seawater groundwater transition zone model, the theory of groundwater movement and numerical simulation. According to the simulation results, this paper summarizes and analyzes the migration patterns of the seawater groundwater transition zone in the study area, and provides effective suggestions for the prevention and control of SWI to this area.

## METHODS

### Sea-land transitional zone model

As seawater exists in the underground aquifers below sea level, there is a contact zone between the freshwater flowing into the ocean and the saltwater. Besides, the saltwater and freshwater are two miscible fluids, and under the effect of hydrodynamic dispersion, the contact zone will exist in the form of a transitional zone (Guo & Huang 2003; Ding *et al.* 2004). Ding *et al.* (2004) proposed the flow process of the freshwater and saltwater in the seawater groundwater transition zone, which is shown in Figure 1.

The sea-land transitional zone model is described by two partial differential equations, one is the flow equation to describe the mixture liquid whose density changes, the other is the advection-dispersion equation to describe salt transportation in saltwater. The distribution range and concentration value of the seawater groundwater transition zone are obtained through two mathematical equations which relate to the density, concentration and water level of the fluid (Xue *et al.* 1992; Wu *et al.* 1996).

### The theory of groundwater movement

This theory mainly includes two equations, the differential equation of groundwater movement and the governing equation of the solute transport model.

*et al.*2009), as described in Equation (1): where is the hydraulic conductivity in the direction of

*x*,

*y*and

*z*;

*h*is the hydraulic head,

*W*is the source/sink term, and is the specific storage of aquifer.

*et al.*(2009) describes the migration and the direction of modeling factor (k) in the three-dimensional groundwater flow system as described in Equation (2): where are the distances along the

*x*and

*y*axes, respectively; is the dispersion coefficient tensor of hydrodynamic force; is the effective porosity; is the dissolved concentration of modeling factor (k); is seepage velocity; is the volume flux of the source sink term in aquifer per unit volume; is the concentration of the substance

*k*in the source and sink term; is the reaction terms.

## RESULTS AND DISCUSSION

### Study area

The study area of Longkou and Zhaoyuan is in the northwest of Jiaodong peninsula in Shandong Province, the southern coast of Bohai Bay, and across the sea to the northeast are cities such as Tianjing and Dalian. The total area is 424 km^{2} and the study area has a maximum horizontal distance of 23 km. The geographical location of the study area is shown in Figure 2.

The average year's annual precipitation is 592.7 mm, which is mainly concentrated in July and August in the monsoon period. The average annual surface evaporation of this area is 1,479 mm. The aquifer of the research area is a multi-structure which can mainly be divided into three layers. Generally, the aquifer consists of sandy soil, silty and silty clay from top to bottom. These data were collected from the Shandong Institute of Geological Survey. Nearly 70% of the urban water supply for research area comes from the groundwater, and the over-exploitation of groundwater resources has caused serious SWI.

### Establishment of the conceptual hydrogeological model

The construction of the conceptual hydrogeological model of this area is mainly based on three aspects of the regional aquifer system, the boundary condition and source and sink terms.

The regional aquifer system: the distribution of underground aquifers are continuous, and they are of a layered distribution in a vertical direction. Underground aquifers in the research area can be generalized into three layers, and the lithology of aquifers is by sandy soil, silt and silty clay, primarily from the surface down.

The boundary condition: vertically, the upper part of the aquifer is the transfer boundary and its lower part is the no flow boundary; horizontally, the study area relates to Bohai Sea on the north side, west side and northwest side, it is respectively bounded by the Huangshui River and the Jie River in the northeast and southwest corners which are generalized as Dirichlet type, and the southern boundary is an inland hilly area in Longkou which is generalized as Neumann type.

Source and sink terms: The main source of recharge to the aquifer are rainfall and runoff from the hilly regions, while the sink from the aquifer mainly includes the evaporation and exploitation of groundwater. According to hydrogeological data, this paper assigns recharge rates for calculation units throughout the study area.

### Calculation and discussion of the groundwater flow field

*et al.*2009): where is the permeability coefficient in the direction of the

*x*axis and is the permeability coefficient in the direction of the

*y*axis, m/d; is the specific yield of the phreatic aquifer;

*Z*is the bed level value of the aquifer, m;

*H*is the phreatic water level, m;

*H*

_{1}is the water level of the first-type boundary point, m;

*q*is the discharge per unit width of the secondary-type boundary, m

^{3}/d/m;

*D*is the range of the calculation area; and are the first-type and second-type boundary.

The groundwater flow field in the study area is calculated by a MODFLOW module in Groundwater Modeling System (GMS) software, the meshes generation of the research area is carried out under the 3DGRID module in GMS software as shown in Figure 3.

Since three layers can be regarded as the homogeneous layer based on the analysis of the hydrogeological condition in the study area, initial values of the hydrogeological parameter (specific yield, permeability coefficient and effective porosity) of each aquifer in the study area are as shown in Table 1.

Aquifer . | Lithology . | Horizontal permeability coefficient (m/d) . | Horizontal permeability coefficient/coefficient of vertical permeability . | Specific yield . | Effective porosity . |
---|---|---|---|---|---|

First layer | Sandy soil | 32 | 4 | 0.09 | 0.28 |

Second layer | Silt | 5 | 3 | 0.059 | 0.15 |

Third layer | Silty clay | 1.3 | 3 | 0.04 | 0.05 |

Aquifer . | Lithology . | Horizontal permeability coefficient (m/d) . | Horizontal permeability coefficient/coefficient of vertical permeability . | Specific yield . | Effective porosity . |
---|---|---|---|---|---|

First layer | Sandy soil | 32 | 4 | 0.09 | 0.28 |

Second layer | Silt | 5 | 3 | 0.059 | 0.15 |

Third layer | Silty clay | 1.3 | 3 | 0.04 | 0.05 |

Taking the period July 2009–June 2010 as the calibration period, and July 2010–June 2011 as the validation period, it was mainly verified whether the dynamic change of the groundwater level in the simulated study area is consistent with the measured dynamic change of the groundwater level.

First, we calibrated the model according to the data of the calibration period, and two groundwater level observation points were randomly chosen to validate the model. The calculated value of 16 groundwater observation points in the calculation area were fitted with the measured values, it turns out that the fitting degree is better and the proportion of the absolute error that is less than 0.5 m is 93.75%. The comparison between the observed value and the calculated value of the water level in two typical groundwater level observation points during the verification period are as shown in Figures 4 and 5.

As can be seen from Figures 4 and 5, the maximum difference value between observation values and calculated values of the water level in all groundwater level observation points is 0.48 m, and the fitting degree between the calculated values and the observation values is acceptable. Besides the groundwater flow field, the equilibrium quantity of the source sink term and the hydrogeological parameter are verified separately, and the fitting degree of each term with the actual situation is good. As a result, it is concluded that the calculations of the established numerical model of the groundwater flow field can meet the accuracy requirements, which reflects groundwater movement correctly and can be applied to simulate the groundwater solute transport in the study area.

Based on the groundwater flow field simulation, we found that groundwater levels that are effected by the rainfall in a certain degree show a periodic fluctuation during the year, which would be high in the flood season and are likely to become low in the non-flood season, and annual changes of groundwater levels in these three layers are very little, the maximum of which is no more than 2 m. Moreover, the hydraulic gradient in three aquifers is small and the direction of the groundwater flow in the study area is from southeast to northwest.

### Numerical simulation and discussion of the SWI condition

Based on the calculation of the groundwater flow field in the study area, the numerical simulation of the migration in the inland groundwater for in the saltwater is carried out because is the most stable and major element in saltwater, which can directly reflect the migration status of the saltwater. In the coastal area of Longkou and Zhaoyuan, the concentration of in the groundwater is more than 200 mg/L and can be regarded as the SWI area (Xue *et al.* 1993; Wu *et al.* 1996). Because of the over-exploitation of the groundwater in the research area, the migration of the sea-land transitional zone is based on the over-exploitation of groundwater.

There are around 12 km of coastline in the northwestern part of the study area, which is regarded as the simulation boundary to study the migration process into inland groundwater based on the dynamic change of , and the initial value of concentration acting on the coastline is 3 g/L. The calculated span of simulation migration is 7,200 days, and every 300 days is a stress calculation period, which is divided into 24 stress calculation periods.

#### Discussion of the concentration changes of in aquifers of the horizontal direction

After the migration of in inland groundwater for 1,000, 3,000 and 6,000 days, the distribution of the concentration (mg/L) in the horizontal direction in the first layer and second layer is as shown in Figures 6 and 7.

As can be seen from Figures 6 and 7, the process of intrusion from saltwater into the inland groundwater in the study area is obvious during the first 1,000 days of the simulation period, and after that, migration of has a slight increase. Furthermore, the concentration distribution of *Cl*^{–} has no change in two layers at the simulation point of 6,000 days. The concentration distribution diagram of after its migration in aquifers for 6,000 days can be used as the basis for dividing the seawater groundwater transition zone. What is more, from saltwater migrating into groundwater is more obvious in the first layer than that in the second layer which is almost not affected by SWI. The serious SWI area is long and narrow and is also small when it is compared with the total area of the study area. So, we can know that SWI in the study area would be steady now if there is no new change and SWI is more serious in some small areas in the first layer than that in the second layer.

#### Discussion of the concentration changes of in aquifers in the vertical direction

To study the condition of SWI in the vertical direction by concentration changes of , selecting four profiles in the model, the vertical concentration distribution diagram of in each profile after the migration for 6,000 days was presented. The location of the selected profile in the model is as shown in Figure 8.

These four profiles, named as the I-I profile, J-J profile, K-K profile and L-L profile respectively, are scattered and representative and the locations are shown as in Figure 8. The vertical concentration distribution of these four profiles at the simulation point of 6,000 days is as shown in Figure 9.

According to Figures 8 and 9, when the movement of has a steady statement, at the simulation point of 6,000 days, it shows that the migration in the groundwater inland of in the saltwater mainly occurs in the sandy soil (the first layer), the silt layer (the second layer) has been less affected by , and has minimal effects on the silty clay layer (the third layer), which means the maximum depth of the inland groundwater polluted by the saltwater intrusion is the thickness of two layers in the study area.

#### Discussion of concentration changes of over the time at different distance from the coastline

To study SWI at different distances from the shoreline, choose six points such as A, B, C, D, E and F which are different distances from the coastline, and draw the concentration curves of at six points based on simulation results to analyze the changing concentration of at each point over time. The points A, B, C, D, E and F are respectively 3, 4.5, 8.4, 12.5, 14.8 and 16 km from the coastline, and the location of each point is as shown in Figure 10.

The concentration changes of over time are as shown in Figure 11 at each point, namely A, B, C, D, E and F.

As we can see from Figure 11, the concentration variation curves of over time at points A, B, C, D, E, F are obviously different. Taking the simulation point of 6,000 days as the steady movement statement of in the aquifer, the concentration value at the final stability is very different between these six points, such as the final concentration of at point A is about 250 mg/L, but the final concentration of at point F is only about 50 mg/L. After the migration of from the saltwater into inland groundwater for 500 days, the concentration of at point B that is away from the coastline of 4.5 km reached 200 mg/L. Considering the distribution map of the concentration in the horizontal direction in the first aquifer after the migration of in inland groundwater for 3,000 days and the concentration curve of over the time at six points, the area that the concentration of is always more than 200 mg/L exists in the coastal areas, and this kind of area has been polluted by the SWI with the width of 1.5–4.5 km, that is the width of the sea-land transitional zone in research area. Xue *et al.* (1992) from Nanjing University carried out research on the SWI condition of this area through a large number of field observations in around 1991, and they considered that the width of seawater groundwater transition zone in this area was about 1.5–2.5 km. However, after the research, this paper considers that the width of the transition zone in the coastal area of Longkou and Zhaoyuan is about 1.5–4.5 km (wider than the width of sea-land transitional zone obtained by Xue Yuqun) which illustrates that after around 30 years, saltwater in this area continues to migrate into the inland, and the saltwater intrusion has not been effectively curbed if the over-exploitation of the groundwater is continuous.

#### Discussion of the influence of inland groundwater level on the transition zone migration

Based on the current groundwater flow field, change the groundwater level of the study area, namely reducing by 1.5 m and raising by 1.5 m, and then carry out the numerical simulation on the migration of from the saltwater in the research area under different groundwater tables. Finally, analyze the influence of inland groundwater level on the migration of the seawater groundwater transition zone in the study area. After the groundwater level is raised by 1.5 m and reduced by 1.5 m, the concentration distribution of in the first layer after migration for 1,000, 3,000 and 6,000 days is as shown in Figure 12.

According to the concentration distribution diagram of in the first layer after migration for 1,000, 3,000 and 6,000 days at different groundwater levels, the influence of the groundwater level on the migration of is mainly shown within the first 1,000 days. While the groundwater level has been raised by 1.5 m, the width of the area where the concentration of is more than 200 mg/L is about 3 km after the migration of from the saltwater in inland groundwater for 1,000 days; under the condition of the current groundwater level, the width of the area where the concentration of is more than 200 mg/L is about 4 km; while the groundwater level has been reduced by 1.5 m, the width of the area where the concentration of is more than 200 mg/L is about 4.5 km after the migration of from the saltwater in inland groundwater for 1,000 days, and the area where the concentration of is more than 200 mg/L is becoming larger with the decrease of the groundwater level. According to the result, we can speculate that the maximum distance of SWI would increase 1 km when the groundwater level drop every 1.5 meters. As a result, SWI is affected by geological conditions, configuration of the strata, soil properties and other factors, the horizontal distance or vertical distance of SWI is not uniform in the coastal strata, the maximum distance of SWI is highly nonlinear with the groundwater level.

## CONCLUSIONS

Combining the sea-land transitional zone model and the theory of groundwater movement, this paper applies a numerical simulation model, GMS, to study the conditions of SWI in the study area by concentration changes of . According to the discussion of the concentration changes of in aquifers of horizontal direction and vertical direction, the discussion of concentration changes of over time at different distances from the coastline, as well as the discussion of the influence of inland groundwater level on the transition zone migration, we can mainly draw five conclusions about SWI in the study area:

- 1.
SWI in the study area mainly occurs in the sandy soil (the first layer), and there is almost no SWI in the other two layers.

- 2.
Based on conservative estimation, the maximum vertical depth of SWI is no more than the sum calculated by the first layer and the second layer.

- 3.
Areas of serious SWI in the study area are quite small compared to the total area, and SWI should be governed by taking some measures as quickly as possible.

- 4.
The horizontal width of SWI in the study area is between 1.5 and 4.5 km.

- 5.
With the decrease of the inland groundwater level, the horizontal and vertical depth of SWI would have a non-linear growth.

After verification of the simulation results by comparing with published research results about the study area, we can consider that the conclusions above are correct to a large extent based on the analysis of field data. Especially, the horizontal distribution of SWI and the relationship between groundwater table decline and SWI are proven to be credible and have high application values. Therefore, according to study results and discussion, it is suggested that relevant measures are taken to raise the groundwater level in this area to control SWI. For instance, we can build up the recharge engineering, make full use of rain-flood resources, and strictly control groundwater exploitation, and prohibit the exploitation of deep groundwater.

## ACKNOWLEDGEMENTS

This work was supported by a hydraulic research project of Shandong province (SDSLKY201402), the National Natural Science Foundation of China (41272325).