Optimal location of cutoff walls for seawater intrusion

In this paper, the simulation–optimization method is used to study the optimal location of cutoff walls for seawater intrusion. The optimization model is based on minimizing the chlorine concentration of two water sources after 50 years. In order to reduce the computational complexity, a Kriging surrogate model simulation is coupled with the optimization model. Finally, a hypothetical case is used to evaluate the accuracy of the surrogate model and the performance of the optimization model. The results show that the outputs of the Kriging surrogate model and the variable density groundwater simulation model for the same cutoff wall design fit well, and the average relative error of the two outputs is only 2.2% which proves that it is feasible to apply the Kriging surrogate model to this problem. By solving the optimization model, the location of the cutoff wall which minimizes the sum of chlorine concentration of the two water sources after 50 years is obtained. This provides a stable and reliable method for the site selection of cutoff walls for future projects intended to prevent and control seawater intrusion.


Introduction
Overexploitation of groundwater lowers the level of groundwater in an aquifer, and the main reason for seawater intrusion is that the groundwater level in the aquifer is lower than that in the adjacent seawater level (Alfarrah and Walraevens 2018). Seawater intrusion will cause shortages of freshwater resources in coastal land, thus affecting the normal water use of local residents, threatening human health, reducing the inherent irrigation capacity for agriculture, and limiting industrial production (Huang et al. 2014). In addition to hindering the economic development of coastal areas, seawater intrusion can also cause harm to the local terrestrial and aquatic ecology in the form of soil salinization, degradation of shrub communities, and deterioration of water quality from mineral water sources (Arslan and Demir 2013;Qi et al. 2019).
There are many countermeasures to prevent seawater intrusion, such as strengthening groundwater resource management, artificial recharge of groundwater, building impervious groundwater curtains, saving irrigation water, and so on (Guo et al. 2013). Among them, the construction of cutoff walls (impervious groundwater curtains) is an engineering measure with good benefits, which is widely used in the field of seawater intrusion . At present, there are many engineering examples of the prevention of seawater intrusion via the construction of underground cutoff walls, both domestically and internationally. In 1988, U.S. Army engineers used one-dimensional and two-dimensional numerical models to design temporary submerged dams in the lower Mississippi River. The final underground dam successfully protected fresh water supplies in New Orleans and controlled salt water intrusion further upstream (Mcanally and Pritchard 1997). In 1997, Nagata and Kawasaki summarized in detail the construction methods of cutoff walls used to increase groundwater storage and control seawater intrusion. Their report mentioned that many cutoff walls have been built in southern Japan, and more cutoff walls are planned in future. These successful cutoff wall projects in Japan demonstrate the superiority of cutoff walls over other schemes, as they not only increase groundwater storage, but also control seawater intrusion (Mundzer 2001). In China, an underground cement curtain with a length of 6 km and an average depth of 26.7 m was built near the estuary of Huangshui River in Longkou City, Shandong Province, which successfully prevented seawater intrusion and formed a groundwater reservoir (Gao et al. 2004). The underground cutoff wall project of Mawan in Qingdao City in Jiaodong Peninsula has effectively solved the deterioration of groundwater quality caused by seawater intrusion (Sun et al. 2006).
Experts in many nations have studied the problem of seawater intrusion with regards to cutoff walls. In one such study, Luyun et al. observed the dynamic changes of residual saline water before and after the cutoff wall was setup through experiments in 2009. They then setup a numerical model. The behavior of the progressive saltwater intrusion wedge and back residual salt water after the construction of cutoff wall was predicted by SEAWAT software. In this study, the influence of cutoff wall height on seawater intrusion prevention is discussed through experimentation and numerical analysis (Luyun et al. 2009). In 2012, Allow used a three-dimensional finite difference model to represent seawater intrusion by the SEAWAT numerical code of Visu-alModflow software. This simulation demonstrated that the cutoff wall is a good solution to seawater intrusion problem. Compared with injection well method, the installation of a cutoff wall is suggested as being more economic, despite its initial high cost, since it would not have to be repaired over time (Allow 2012). In 2016, Strack et al. provided an approach for reducing saltwater intrusion in coastal aquifers by artificially reducing the hydraulic conductivity in the upper part of selected areas by using a precipitate. By presenting equations for the location of the tip of the saltwater wedge and verifying them through a sand-tank experiment, this excavation-free method of barrier installation proved to be feasible. Further development and experimentation, however, is needed to evaluate the practicality of this method (Strack et al. 2016). In 2017, Abdoulhalik et al. proposed a mixed physical barrier (MPB), which combines a cutoff wall and a semi-permeable underground dam, as a new barrier system to control seawater intrusion. They used the numerical program SEAWAT to evaluate the consistency of experimental data and to test the sensitivity of blocking performance within various key parameters. The results showed that an MPB can be more effective at controlling seawater intrusion at a lower cost than traditional physical barriers . Recently, Li et al. have conducted the comparative studies using sandbox test as a physical model and the FEFLOW-based groundwater numerical model to better understand the dynamic processes of saltwater intrusion with and without subsurface barriers. They also studied the influence of the subsurface barriers with different permeability to the migration distance and diffusion of saltwater front by testing in the laboratory and using the two-dimensional numerical simulation model. The result indicated that the subsurface barriers with low permeability (K = 3.7 × 10 -8 m/s) could prevent the saltwater intrusion effectively (Li et al. 2018a, b).
In the design of cutoff wall, the predecessors mostly studied the material or permeability of cutoff wall, but seldom studied the location of cutoff wall. Regarding the site selection of the cutoff wall, most of the previous studies performed qualitative analysis by moving the position of cutoff wall, but did not provide quantitative calculation methods for selecting the optimal position of the cutoff wall. In 2018, Wu Yajie and others used OpenGeoSys software to simulate saltwater intrusion under the influence of an underground cutoff wall. They revealed the saltwater intrusion rule for intrusion under the influence of underground cutoff wall and analyzed the influence of location of the cutoff wall with regards to its effectiveness in preventing saltwater intrusion by comparing six different positions of cutoff wall .
Therefore, on the basis of previous studies, this paper considers a hypothetical problem about the site selection of a cutoff wall in the face of seawater intrusion. It demonstrates how to choose the position of a 220 m cutoff wall with the same height and direction such that it can minimize the amount of chlorine concentration at two water supply sources over a 50-year period, on the premise that the chlorine concentration at each water supply source after 50 years does not exceed the concentration value after 50 years without the cutoff wall. By establishing a threedimensional variable density numerical simulation model of seawater intrusion, the simulation-optimization method is used to study the optimal scheme of cutoff wall location in the hypothetical problem. In order to reduce the computational load caused by invoking the simulation model, this paper establishes a Kriging surrogate model with reliable accuracy to replace the simulation model and couple it with the optimization model. This study provides a stable and reliable method for site selection of cutoff walls in future prevention and control of seawater intrusion.

Study area
Our study area is a simplified hypothetical example based on the seawater intrusion area in Longkou, China in order to better present the method of selecting site for cutoff walls. In Longkou area of China, the thickness of aquifers is small and the shortage of fresh water resources is severe. The longterm overexploitation of groundwater causes the groundwater level to drop sharply, which destroys the original hydraulic balance between the salty water body and the fresh water body. At the same time, the marine sediments near the coastline of Longkou area are fine sand, coarse sand and fine gravel with good water permeability, and there is a good hydraulic connection between the sea water and the coastal aquifer. Therefore, the seawater intrusion in Longkou area has been serious for a long time. Local decision makers consider the construction of cutoff walls as one of the options available to control the seawater intrusion. In this case, we chose Longkou as the prototype to design a hypothetical area to study the optimal location of cutoff wall.
In the water flow model, the northern part of the study area is bounded by the foot of a mountain, and is generalized to the known flow boundary. Both the eastern and western sides are bounded by watersheds and can be generalized as known flow boundaries. The southern part is the boundary between seawater and groundwater. Considering the influence of sea level rise on seawater intrusion, the seawater head rises by an average of 0.01 m annually, so it can be generalized to a known head boundary. In the solute transport model, the northern boundary is generalized to the known hydrodynamic dispersion flux boundary, the eastern and western boundaries are generalized to the known hydrodynamic dispersion flux boundary, and the southern boundary is the known concentration boundary. The upper boundary is the phreatic surface with water exchange and no solute transport, while the lower boundary is the aquifer impervious floor without water exchange and solute transport. There is no surface water in the area. Under the control of landform and topography, the direction of groundwater flow is generally from north to south.
The movement factors of groundwater flow in the study area change with time. The velocity of groundwater flow cannot be neglected in both vertical and horizontal directions. Therefore, groundwater flow can be generalized as three-dimensional unsteady flow. The stratum lithology in this area is coarse sand, mainly loose sediments. From north to south, sandstone particles have the characteristics of changing from coarse to fine, and the permeability coefficient also changes from large to small. The internal structure of the aquifer can be generalized as a heterogeneous isotropic phreatic aquifer. According to the internal structure of the aquifer, the study can be divided into two parameter zones, the north and south (Fig. 1). The hydrogeological parameters of each zone are as follows: Table 1. Atmospheric precipitation infiltration recharge and lateral runoff recharge are the main recharge modes of regional aquifers, and the main discharge modes are phreatic evaporation and artificial exploitation. There are six pumping wells and two water supply sources in the study area. During the process of groundwater exploitation, the trend of seawater intrusion appears. A three-dimensional numerical model of variable density seawater intrusion is established to predict the results of seawater intrusion in the study area after 50 years, which can be used as the basis for the design of cutoff wall in the study area.

Establishment of mathematical model
Considering the transition zone, the problem of seawater intrusion needs to be described by two equations (Zeng et al. 2018). One is the governing equation of water flow which reflects the change of concentration (density) of solution and the change of head and velocity of the liquid (mixture of seawater and freshwater). The other is the governing equation of water quality which reflects salt transport in groundwater. The flow control equation and the water quality control equation, together with the corresponding definite solution conditions, constitute the mathematical flow model and the mathematical solute transport model of a seawater intrusion aquifer. Because the concentration term is contained in the flow control equation, and the actual flow velocity in the water quality control equation needs to be calculated by  water head, the above two equations need to be solved by the motion equation.
The change of groundwater density in the transition zone leads to the constant change of the actual head value under the same pressure. For this reason, the concept of equivalent fresh water head is used to establish the flow equation: h = p ρ 0 g + z , where 0 is freshwater density. At this time, the mathematical model of flow is as follows: where K is the hydraulic conductivity of the aquifer (m/d); c is the solution concentration (mg/l); n is porosity; η = ε c s is the density coupling coefficient, c s is the maximum density, ε = ρ s −ρ 0 ρ 0 is the density difference, ρ s is the corresponding liquid concentration; μ s is the water storage rate (m −1 ); q is the volume (or sink) flow of a porous medium with a unit volume (d −1 ); L 1 is the known flow boundaries on the northern side; L 2 , L 3 are the water barriers on both the eastern and western sides; L 4 is the sea boundary; h 1 (x, y, z, t), h 0 (x, y, z) and q 0 (x, y, z) are known functions. The mathematical models of solute transport in seawater intrusion are as follows: where the direction of the X axis is consistent with the mean velocity; D xx , D yy , D zz are diffusion coefficient tensors (m 2 /s); u x , u y , u z are the actual average velocities of groundwater in the respective X, Y and Z directions (m/s); c * is the concentration of salt in the liquid for extraction or injection (mg/l); c 1 (x, y, z) and c 0 (x, y, z) are known functions.
The equations of motion for solving coupling problems are: is the hydraulic conductivity under fresh water conditions (m/d).
The mathematical water flow model, solute transport model and motion equation listed above constitute a complete mathematical model of seawater intrusion in the study area.

SEAWAT
The mathematical model of seawater intrusion mentioned above is solved by the existing program SEAWAT. SEAWAT can simulate the three-dimensional, variable density and unstable groundwater flow system in porous media (Cobaner et al. 2012). The source code of SEAWAT is developed by combining MODFLOW and MT3DMS into a single program, which can solve the coupled equation of water flow and solute transport. By replacing the flow volume in the matrix equation with the flow mass and introducing the appropriate density term, the modified MODFLOW can solve the variable density flow equation. Flow density is usually determined as a function of solute concentration, while the effect of temperature on flow density is not considered temporarily. Using the program in MT3DMS, SEAWAT can simulate the spatiotemporal variation of salt concentration. SEAWAT uses explicit or implicit methods to couple groundwater flow equations with solute transport equations. In the explicit method, the flow equation is solved for each time step, and the advection velocity can be used to solve the solute transport equation. This method will be repeated until the end of each stress period and simulation period. When the implicit method is used for coupling, the equations of water flow and solute transport will be solved many times in the same time step, until the maximum difference between the two successive iterations of water flow density is less than the user specified threshold.

Model prediction
The study area is based on 20-m × 20-m squares as the basic unit which is evenly divided into 50-unit × 100-unit grids. The simulation time is 50 years, and the monitoring frequency is once a year. When the water quality of underground aquifers is deteriorated due to seawater intrusion, the most obvious change is the increase of chloride ion concentration in the water. Therefore, chloride ion is usually used as an indicator to determine whether seawater invades and the extent of invasion. Using groundwater chlorine concentration exceeding 250 mg/l as the mark of seawater intrusion, Fig. 2 shows the distribution of seawater intrusion after 50 years of simulation. Table 2 shows the chloride concentration values of the two water sources in each period.
According to the calculation results, it can be found that the No. 1 water supple resource in the study area will soon face the problem of seawater intrusion after 30 years.
The chlorine concentration in the water source area is 206.24 mg/l, and the chlorine concentration in the No. 2 water supple resource is 21.90 mg/l. Fifty years later, the No. 1 water supple resource in the study area has been completely damaged by seawater intrusion. The chlorine concentration at that water supply resource is 1063.35 mg/l, while the No. 2 water source area is at 98.34 mg/l. In order to verify the prevention effect of a cutoff wall on seawater intrusion, a cutoff wall crossing the East-West boundary is designed at 400 m away from the coastline, and plugged into the established seawater intrusion model for the calculation. The final result is shown in Fig. 3.
The results show that the cutoff wall effectively blocks the seawater within 400 m from the coastline by cutting off the seawater intrusion channel, and that the two water sources are not polluted. However, considering that the distance between the two watersheds is tens of kilometers under actual conditions, it is unrealistic to build such a long cutoff  wall, so the above design scheme for a cutoff wall would be difficult to achieve in practice. Therefore, the simulation-optimization method is used to discuss the design of a 220 m long cutoff wall in the study area. It requires that the sum of chlorine concentration at the two water sources should be minimized over a 50-year period, and that the chlorine concentration at both water sources at the 50-year mark should not be higher than that of natural seawater intrusion. The study of this problem provides a reference method for site selection of cutoff wall in actual situations.

Simulation-optimization method
In order to study the site selection for a cutoff wall using a simulation-optimization method, it is necessary to setup the optimization simulation such that it will maximize the achievement of the stated objectives within the given laws of groundwater movement. During the process of solving the optimization model, the simulation model needs to be mobilized many times, which creates a huge computational load. To help alleviate this problem, the Kriging method (Zhou et al. 2018) is introduced in order to construct a surrogate model to replace the seawater intrusion simulation model, which can then be coupled with the optimization model, greatly saving on calculation costs.

The establishment of the surrogate model
Using a surrogate model is an effective way to couple simulation models and optimization models Chu and Lu 2015;Hou et al. 2017). It replaces the input-output response relationship of the simulation model, and can describe the input-output data of the simulation model with a higher accuracy with a smaller number of calculations. In the hypothetical case, the input is the design location of the cutoff wall, and the output is the sum of chlorine concentration of the two water sources after 50 years. In this paper, the Kriging method is used to build the surrogate model. It is an unbiased estimation method based on minimum estimated variance (Li et al. 2018a, b). The surrogate model established by the Kriging method can avoid the computational load caused by multiple invocation of the simulation model on the premise of guaranteeing high accuracy in describing the function of the simulation model. The basic form of Kriging surrogate model is (Kwon et al. 2014;Kersaudy et al. 2015): where y(x) is the pollutant concentration value required in the substitute model, and ŷ(x) is the estimate value of y(x) , which can be divided into two parts, f(x) T β is the linear regression part, and z(x) is the random is the base function of the known regression model, the quadratic function-type base function is selected in this paper. The undetermined parameter β = [β 1 , β 2 , ⋯ , β k ] is the coefficient corresponding to the Fig. 3 Results of seawater intrusion after cutoff wall design from 400 m location of coastline base function, which can be obtained by using the prepared training data. The random part of z(x) satisfies the following conditions: where R(x i , x j ) is a spatial correlation equation between x i and x j at any two sampling points, i.e., a correlation function. There are many forms of correlation function. The EXP model is adopted in this paper: where θ k is an undetermined parameter; and x i k is the K dimension coordinate of the i sample. According to the Kriging model, the prediction value of the response value, y(x) , at the prediction point, x, is: where r(x) is the correlation vector between point x and n sampling points of training samples (x 1 , x 2 , ⋯ , x n ) , r(x) = [R x, x 1 , R x, x 2 , ⋯ , R(x, x n )] ; y is the response value of pollutant concentration corresponding to n sampling points, which is a vector of n × 1 ; β is the undetermined component parameter of the linear regression, which can be obtained by the optimal linear unbiased estimation: β = (f T R −1 f) T f T R −1 y ; R is a n × n rank correlation matrix consisting of n sampling point correlation coefficients: The variance 2 estimate is determined by the following formula: Therefore, the establishment of the Kriging surrogate model is actually the solution of the above nonlinear unconstrained optimization problem (Hou et al. 2015). After the undetermined parameter θ k is obtained, the pollutant concentration response value can be obtained by establishing the Kriging model. The θ k can be obtained through an unconstrained optimization formula: Based on the construction principle of the Kriging surrogate model, this paper uses the software of MATLAB to compile and debug the surrogate model. The procedure of invoking the model is simple, and the workload is greatly reduced. The Latin hypercube sampling method (Davey 2008) was used to sample the coordinates of the center of the cutoff wall. Because the cutoff wall cannot effectively protect the water quality of the No.1 water source when it is built to the north of the No.1 water source, only locations for the cutoff wall between the No.1 water source and coastline are considered. Therefore, the value of the abscissa is [100,900], and the value of the ordinate is (0,700]. Forty sets of uniformly distributed coordinates were selected as input values of the samples, and the 220-m cutoff wall centered on the coordinates was substituted into the simulation model. The concentration sum of chlorine at both water supply sources 50 years later was calculated by the SEAWAT program as output values for each of the samples. The 40 sets of input and output data are substituted into the compiled MATLAB program, and the Kriging surrogate model is trained. The Latin hypercube method is used to extract 10 sets of data to test the Kriging surrogate model after training, so as to ensure the accuracy of the model. The average relative error and correlation coefficient are used to measure the degree of fit between the surrogate model and the simulation model. The test results are shown in Fig. 4. According to the analysis of the test results, the average relative error between the surrogate model and the simulation model is only 2.2% for the output of 10 groups of test samples, which shows that the Kriging surrogate model is feasible for the selection of cutoff wall location and can be used to replace the function of the seawater intrusion simulation model. In order to more directly demonstrate the superiority of establishing surrogate models, the calculation time of the two processes is recorded. It takes about 115 s to run a simulation model, so it takes about 128,340 s (35.6 h) to solve the optimal site selection problem using the simulation-optimization method. In contrast, it only takes 0.15 s to run the surrogate model, and the process of establishing the surrogate model runs 50 times, which takes about 5750 s (1.6 h). Therefore, replacing the simulation model with the Kriging method saves about 34 h. When the method is applied to a wider range of actual case, its advantages will be more obvious.

Establishment and solution of optimization model
The objective of this study is to minimize the sum of chlorine concentration at two water sources after 50 years. The 0-1 integer programming model is used to describe the problem. The range of choosing the center coordinate of cutoff wall is mentioned above. There are 1116 grids in this range. When the grid is 0, the center position of a noncutoff wall is indicated. When the grid is 1, the center of a cutoff wall is represented. The optimization model generally consists of three parts: decision variables, objective functions, and constraints. The specific models of this case are as follows: where F indicates the minimum sum of chlorine concentration at two water sources after 50 years, C 1 and C 2 denote chlorine concentration at the No.1 and No.2 water sources after 50 years, respectively; j denotes grid number and x j is the decision variable. The first constraint condition represents the established surrogate model, which reflects that the optimization model conforms to the inherent law of the groundwater system in the study area. The second constraint means that the concentration of chlorine in the case with a cutoff wall is not higher than that in the case of natural seawater intrusion 50 years later. The third constraint condition indicates that there is only one central point of the cutoff wall, that is, only one cutoff wall is setup. The fourth constraint condition is the constraint of decision variables. Genetic algorithm (Sedki and Ouazar 2011;Sreekanth et al. 2017) is used to solve the optimization model, and the grid number is converted to the coordinate range.
By solving the optimization model, the optimal solution for the coordinates of the center point of the cutoff wall is finally set at: x ∈ [450, 470], y ∈ [670, 690] . Therefore, the optimal design location of cutoff wall in this case is obtained. Table 3 demonstrates the coordinate range of the complete cutoff wall and Fig. 5 shows the optimum design location of cutoff wall. More specifically, the simulation-optimization coupling model finds the best design position of cutoff wall in 1116 grids, and the central position of this cutoff wall is in a square grid of 20 m × 20 m ( x ∈ [450, 470], y ∈ [670, 690] ). Correspondingly, the horizontal coordinate range of the 220 m cutoff wall is [350,570]. In the longitudinal direction, the position of cutoff wall should be somewhere between 670-690. Decision makers can choose based on  Table 3 and Fig. 5.

Result analysis
The location of cutoff wall has an important influence on the seawater intrusion of the two water supply sources. Inappropriate locations may cause the seawater intrusion in the two water sources to become more severe. The simulation-optimization method can effectively find out the optimal position of cutoff wall that meets the requirements. The optimization goal of this hypothetical case is to minimize the sum of chloride ion concentration in the two water supply sources. According to the results of the four contrast schemes (Table 3), when the cutoff wall of the optimal scheme is shifted to the east or the coastal zone, although the chloride ion concentration of the No. 2 water source will decrease, the chloride ion concentration of the No. 1 water source will increase sharply, thus leading to an increase in their sum. On the other hand, when the cutoff wall shifts westward, the chloride ion concentration in the No. 2 water source will increase sharply, even exceeding the value without the cutoff wall, which does not meet the requirement of the second constraint.
In fact, decision makers can adjust the objective function and constraints based on local needs, and the position of the optimal cutoff wall obtained at this time will change accordingly.
The optimal scheme is run into the simulation model to check the results of seawater intrusion, as shown in Fig. 6. Table 4 shows the chloride concentration values of the two water sources in each period with the optimum design location of cutoff wall. Comparing the results of Table 4 with the results of natural seawater intrusion without designing a cutoff wall (Table 2), it is not difficult to see that neither of the two water supply sources was harmed by seawater intrusion within 30 years after the cutoff wall was installed according to the optimal results. In fact, 40 years later, the transitional zone of seawater intrusion had just begun to extend to the No.1 water supply resource. It can be seen that the optimal layout scheme protected both water supply sources from pollution for a long time, and ensured the minimum sum of chlorine concentration at the two water supply sources after 50 years. The calculation shows that the sum of  chloride ion concentration in the two water sources after 50 years can be reduced by 79.59% by using the optimal layout of cutoff wall compared with that before the construction of cutoff wall.

Conclusions
The construction of underground cutoff walls can effectively prevent seawater from invading terrestrial aquifers. Under ideal conditions (that is, when the cutoff wall completely cuts off the channel of invading seawater), seawater cannot invade the upstream area of the cutoff wall. The surrogate model established by the utilization of the Kriging method can obtain an input-output relationship similar to the variable density groundwater numerical simulation model with a relatively small amount of calculation. In the optimization of cutoff wall location, the application of the Kriging model instead of the simulation model to couple with the optimization model can greatly reduce the calculation load. In the case study, the average relative error between the Kriging surrogate model and the groundwater numerical simulation model is only 2.2%, which saves 34 h of calculation time. Based on the 0-1 integer programming, the optimization model takes the central position of the cutoff wall as the decision variable, which can minimize the sum of chloride ion concentrations at the two water sources after 50 years within the given constraints. In the case study, by solving the optimization model, it is found that the optimal center position of the 220-m cutoff wall is within the region of x ∈ [450,470], y ∈ [670,690] . The site selection scheme reduces the total chloride ion concentration of the two water sources by 79.59% after 50 years compared with that before the construction of the cutoff wall. This paper provides a stable and reliable quantitative analysis method for site selection of seawater intrusion cutoff walls.