## Abstract

In this study, a linked simulation–optimization approach is proposed for optimal dewatering system design for an excavation site. In the simulation part of the proposed approach, an analytical groundwater flow model is used to simulate the flow process in the subsurface. This model is then integrated with two different optimization models where differential evolution (DE) and harmony search (HS) optimization algorithms are used. The objectives of both models are to determine the best dewatering system configuration by considering locations and pumping rates as decision variables. During the search process, the constraints related to drawdown and well locations are included in the optimization by using the penalty function approach. The performance of DE- and HS-based solution approaches is evaluated on a hypothetical excavation site for different well numbers. Identified results indicated that the DE-based solution approach provides slightly better results than those obtained by using the HS-based solution approach.

## INTRODUCTION

As a result of continuous development, intensive urbanization and the growth of the population, a large decline has been observed in the number of suitable areas for construction projects. Therefore, regions that are less suitable in terms of soil properties are being chosen as construction sites. One of the most encountered problems in these areas is the high groundwater levels which not only cause slope stability problems but also prevent the ensuring of suitable working conditions. Therefore, high groundwater levels around construction sites need to be lowered before starting excavation work. This task is usually conducted by means of dewatering systems consisting of pumping or dewatering wells and a collector line. The cost and efficiency of the dewatering systems mostly depend on the numbers, depths and pumping rates of the dewatering wells. Therefore, the following questions regarding dewatering system installation can be asked before starting excavation work: *How many dewatering wells should be used for lowering the groundwater table below the bottom of the excavation site? Where should these dewatering wells be located over the field? What amount of groundwater should be pumped from each dewatering well?* These questions are usually answered by the site engineers depending on their previous experiences. However, this kind of system configuration may not be optimal in terms of the installation and operation costs. Furthermore, such a dewatering system may not result in the best solution in terms of ensuring safe working conditions. Therefore, optimum configuration of a dewatering system is determined by means of a linked simulation–optimization model.

Linked simulation–optimization models have been widely used in groundwater management for many years. These models consist of separate optimization and simulation parts such that while the optimization parts aim to generate the best management alternative from the set of solutions, the simulation parts determine the response of the groundwater system for the generated alternative (Ayvaz 2016). Note that available simulation–optimization models in the literature differ in terms of the considered simulation and optimization approaches. The first difference is observed in the simulation parts, which are the most important parts of these models since they transform the generated pumping plan to hydraulic head or drawdown fields over the aquifer system. This task is achieved by solving the governing partial differential equations (PDEs) depending on the given initial and boundary conditions. Different analytical and numerical solution approaches are proposed in the literature to solve these PDEs in the linked solution approaches. The second difference of these models is observed in the optimization parts, which are also important since they are used to determine the best solution from the set of alternatives. In early studies, groundwater management problems were usually solved by using deterministic optimization approaches. Although these approaches can determine a global optimum solution in reasonable solution times, they may not guarantee to find the global optimum for cases with strong nonlinear or non-convex solution spaces (Ayvaz 2009). Therefore, use of stochastic-based heuristic optimization approaches is preferred nowadays in order to solve these kinds of problems. There exist a huge number of heuristic optimization approaches used in groundwater management studies. Review of these approaches and state-of-the-art reports can be found in Gorelick (1983), Yeh (1992), and Das & Datta (2001).

The literature currently includes several studies regarding the solution of dewatering optimization problems. Aguado *et al.* (1974) proposed a method for designing a dewatering system for a site. Their simulation model consists of the finite difference solution of the governing PDE, which is embedded as a constraint in a linear-programming-based optimization model. Aguado & Remson (1980) developed a model to solve the dewatering optimization problem for maintaining low water levels in an excavation site. They considered a fixed-charge problem with installation and pumping costs to be minimized. Galeati & Gambolati (1988) developed a solution model which integrates a three-dimensional finite element model with an optimization approach where mixed integer linear programming is used. They evaluated their model by considering different objective functions and constraints. Önder & Değirmenci (2000) developed a management model which integrates the simulation and optimization approaches for dewatering an excavation site. Their model simulates the flow process by means of a response matrix which is developed by using the outputs of the MODFLOW model. Their optimization model solves the problem by using linear and mixed integer linear programming approaches. Tokgöz *et al.* (2002) also used linear and mixed integer linear programming approaches to determine the optimal dewatering system design for a trench-type excavation. They considered both steady-state and transient model simulations using a response matrix. Demirbaş *et al.* (2008) developed an approach for dewatering an excavation area by using a combined simulation–optimization model. Their simulation model is developed for an unconfined aquifer by means of the response matrix approach and integrated with a genetic algorithm (GA) optimization approach. Ayvaz *et al.* (2011) developed a solution approach which integrates the heuristic harmony search (HS) optimization algorithm with an analytical groundwater flow model. They evaluated the model performance on a construction site for different well numbers. Jiang *et al.* (2013) developed a simulation–optimization model to solve a dewatering optimization problem in Shengli No.1 open-pit coal mine, China. Their simulation model was developed by using the MODFLOW model for steady-state flow conditions. The developed simulation model is then integrated to an optimization model where a Pareto dominance-based real-coded GA is used. They considered the safety of the mine, slope and seepage discharge of dewatering wells in the objective function. Shi *et al.* (2017) developed a solution approach which considered the dynamics process of construction dewatering. Their approach calculates the total water yield and dewatering time so that the initial groundwater level reduces to the target level. Shourian & Davoudi (2017) developed a simulation–optimization model to solve the dewatering optimization problem in urban areas. Their model consists of the MODFLOW model as the groundwater flow simulation tool and firefly optimization algorithm (FOA) as the optimization tool. They tested their model performance in an urban area in the southeastern part of Iran. The model results indicated that the cost-effective design noticeably outperformed the consulting engineers' proposals.

It should be noted that most of the studies given above considered a predefined number of fixed well locations. However, since the locations of the wells are also unknown just as their numbers and pumping rates, it is required to determine the optimum system characteristics including numbers, locations, and pumping rates of the dewatering wells by means of the simulation–optimization models. Therefore, the main objective of this study is to solve the construction site dewatering optimization problem by using differential evolution (DE) and HS optimization approaches. The objective of both approaches is to minimize the total pumping rate of the dewatering wells by considering the locations and pumping rates as the decision variables. In order to determine the groundwater system behavior for the generated system configuration, an analytical groundwater flow model is used as the simulation model. During the optimization process in both DE and HS, it is required to maintain some physical and managerial constraints, which are included in the optimization model by means of the penalty function approach. The performance of the proposed DE- and HS-based solution approaches are evaluated on a hypothetical excavation site for different numbers of dewatering wells and the identified results are compared in detail.

## SIMULATION MODEL

*h*is the hydraulic head value at the piezometer location after pumping ,

*H*is the hydraulic head value before pumping ,

*K*is the hydraulic conductivity of the aquifer ,

*R*is the radius of well influence , and

*r*is the measured radial distance from piezometer to well . After calculation of hydraulic head using Equation (1), the drawdown value at the piezometer location can be calculated as follows:

As can be seen from Equation (2), the drawdown value is the difference between the initial hydraulic head and the hydraulic head during pumping. It is important to note that this drawdown value is the result of a single pumping event on the aquifer system. In the case of more pumping wells in the aquifer system, by depending on the principle of flow superposition, the drawdown at any point can be calculated by aggregating all the drawdowns caused by the independent operations of each well (Kresic 1997; Usul 2008). Figure 2 illustrates this process for two pumping wells in an unconfined aquifer where the resulting drawdown is the sum of the two individual drawdowns . By using this principle in the simulation model, the resulting groundwater table elevations for a dewatering system configuration can be calculated by using the optimization model.

## OPTIMIZATION MODEL

As indicated previously, the optimum dewatering system configuration is determined by using both DE- and HS-based optimization models. The main reason for selecting DE and HS as the optimization approaches is their strong capability of finding the optimum solutions. Furthermore, they are not required to perform complex coding/decoding processes during programming since their computational structures are very simple compared with other models such as GA. In the following sections, after describing the main computational steps of DE and HS, how to solve a dewatering system optimization problem using these approaches is explained.

### DE optimization algorithm

DE, first proposed by Storn & Price (1997), is a population-based evolutionary optimization algorithm. Basically, it has a similar computational scheme to GA where the given solutions are evolved by means of mutation, crossover, and selection operators. The key difference between DE and GA is that all the individual solutions in DE are subjected to genetic evolution whereas the same task is performed based on a probability in GA. After this evolution process, all the evolved solutions are directly transferred to the next generations if their corresponding objective function values are improved. The basic computational steps of DE can be summarized as a pseudo-code as follows (Gökçe & Ayvaz 2014a):

Randomly initialize all agents (e.g. candidate solutions) in the population ( being the number of individuals in the population).

Repeat the following until a termination criterion is met:

• For each agent in the population:

• Randomly select three distinct solutions

*a*,*b*and*c*from the population.• Pick a random index ( being the dimension of the problem).

• Compute the agent's potentially new position as follows:

• For each

*k*, pick a uniformly distributed random number, .• If ( is the crossover rate) or then set ( is the differential weight) otherwise set .

• If replace the agent in the population with the improved candidate solution, that is replace with in the population.

Pick the agent from the population that has the highest fitness or lowest cost and return it as best found candidate solution.

### HS optimization algorithm

HS, first proposed by Geem *et al.* (2001), is a heuristic optimization algorithm which is created by inspiration from the musical improvisation process. In this analogy, each musician in an orchestra corresponds to a decision variable, and similarly, the notes in the musicians' memories correspond to the potential values of the decision variables. Note that the primary objective of the musicians is to find a fantastic harmony based on the following musical operations (Geem *et al.* 2001): (i) playing a note from the harmony in memory; (ii) playing a note which is close to another one stored in memory; (iii) playing a note randomly from the possible note range. These musical operations can be adapted to engineering optimization problems as follows: (i) new variable values are selected from harmony in memory; (ii) new variable values are further replaced with the ones which are close to the current ones; (iii) new variable values are randomly selected from the possible range. The combination of these three operations allows searching for a global optimum solution. The following computational scheme describes the required solution steps for solving an optimization problem using HS (Gökçe & Ayvaz 2014b):

*Step 1:*Generate random solution vectors as many as the harmony memory size , then, store them in harmony memory .*Step 2:*Generate a new solution vector . For each element :• with probability of (harmony memory considering rate), pick the stored value from such that where is the uniform random number;

• with probability of , pick a random value within the allowed range.

*Step 3:*Perform additional adjustment if the value in Step 2 came from :• with probability of (pitch adjusting rate), change the value of by a small amount such that where is the fret width, which can be defined as the amount of the maximum change in pitch adjusting process;

• with probability of , do nothing.

*Step 4:*If value of is better than the worst vector in , replace with .*Step 5:*Repeat from Step 2 to Step 4 until termination.

### Problem formulation

*x*and

*y*directions, respectively.

As can be seen from the optimization formulation given above, the objective of the DE- and HS-based optimization models is to minimize the sum of the pumping rates of the dewatering wells. This kind of an objective function is considered rather than the dewatering configuration cost since the cost of the dewatering system is already related to the numbers and the pumping rates of the dewatering wells. For such a problem, the decision variables of the optimization model are , , and for all . The problem is solved for different well numbers and the identified results are evaluated and compared together.

## NUMERICAL APPLICATION

In this section, the applicability of the proposed DE- and HS-based solution approaches is evaluated on a hypothetical excavation site dewatering problem given in Figure 3. This is a shopping mall construction site which has a deep foundation with a depth of . Hydrogeological studies around this site indicate that the groundwater table is located below the ground surface. The aquifer system is classified as an unconfined aquifer with a uniform thickness of . The hydraulic conductivity and the radius of the well influence are determined as and , respectively. In order to measure the variation of the groundwater table elevations, this site is surrounded by nine piezometers, all of them fully penetrating the entire aquifer thickness. In order to provide safe working conditions inside the excavation site and to prevent slope stability problems, groundwater table elevations must be lowered below the bottom of the excavation site. This task should be conducted by drilling new dewatering wells equipped with submersible pumps with a maximum pumping capacity of . Since these pumps will be operated during the entire construction process, it is essential to minimize their pumping rates in order to solve the dewatering problem with minimum operation cost. Therefore, the objective here is to determine the numbers, locations, and pumping rates of these submersible pumps by minimizing Equation (3) subject to the constraints in Equations (4)–(8). The related optimization parameters of DE are: , and . Similarly, the parameters for HS are taken as: , , and . Values of these parameters are selected depending on previous experiences and recommendations in the literature. In order to guarantee convergence for each solution, the maximum generation number in DE is taken as 10,000, which corresponds to 200,000 simulation model runs. Thus, the maximum improvisation number of HS is taken as 200,000 in order to obtain comparable results. For both DE and HS, the bounds of the decision variables and the values of the penalty coefficients are taken as follows: . Using these parameter sets, the problem is solved for to 5 dewatering wells for both DE and HS. Convergence plots for each solution are given in Figure 4.

As can be seen from Figure 4, as the number of dewatering wells increases, the complexity of the problem also increases, which is why more simulation model executions are required for convergence. However, it is not required to conduct 200,000 model executions since the most complex models of DE and HS require approximately 85,000 and 108,000 model executions, respectively. For these convergence profiles, Table 1 lists the identified results.

Number of wells . | Objective function value . | Total penalty value . | Identified pumping rates . | |||||
---|---|---|---|---|---|---|---|---|

. | . | . | . | . | ||||

DE | 1 | 9,973.78 | 9,933.78 | 40.00 | – | – | – | – |

2 | 3,538.91 | 3,458.91 | 40.00 | 40.00 | – | – | – | |

3 | 118.18 | – | 40.00 | 38.38 | 39.80 | – | – | |

4 | 113.09 | – | 37.05 | 24.63 | 11.40 | 40.00 | – | |

5 | 108.78 | – | 40.00 | 5.89 | 21.59 | 13.94 | 27.36 | |

HS | 1 | 9,973.79 | 9,933.79 | 40.00 | – | – | – | – |

2 | 3,539.27 | 3,483.02 | 40.00 | 40.00 | – | – | – | |

3 | 118.37 | – | 39.10 | 39.64 | 39.63 | – | – | |

4 | 113.46 | – | 32.29 | 32.76 | 20.77 | 27.64 | – | |

5 | 111.91 | – | 26.39 | 5.29 | 24.06 | 22.43 | 33.74 |

Number of wells . | Objective function value . | Total penalty value . | Identified pumping rates . | |||||
---|---|---|---|---|---|---|---|---|

. | . | . | . | . | ||||

DE | 1 | 9,973.78 | 9,933.78 | 40.00 | – | – | – | – |

2 | 3,538.91 | 3,458.91 | 40.00 | 40.00 | – | – | – | |

3 | 118.18 | – | 40.00 | 38.38 | 39.80 | – | – | |

4 | 113.09 | – | 37.05 | 24.63 | 11.40 | 40.00 | – | |

5 | 108.78 | – | 40.00 | 5.89 | 21.59 | 13.94 | 27.36 | |

HS | 1 | 9,973.79 | 9,933.79 | 40.00 | – | – | – | – |

2 | 3,539.27 | 3,483.02 | 40.00 | 40.00 | – | – | – | |

3 | 118.37 | – | 39.10 | 39.64 | 39.63 | – | – | |

4 | 113.46 | – | 32.29 | 32.76 | 20.77 | 27.64 | – | |

5 | 111.91 | – | 26.39 | 5.29 | 24.06 | 22.43 | 33.74 |

As can be seen from Table 1, for and , the identified pumping rates for both DE and HS are , which is the value. For these solutions, the final objective function includes penalty values in terms of well location suitability and drawdown constraints. These results indicated the dewatering optimization problem cannot be solved by using one or two pumping wells if their maximum rates are fixed at . For to 5, the final objective function values have no penalties in both DE and HS, which indicates that the dewatering optimization problem can be solved without any violations in terms of Equations (7) and (8). When the objective function values are compared, it can be seen that DE provides slightly better results than those obtained in HS for all solutions. The final values decrease with inclusion of additional pumping wells. However, this improvement is not significant. For identified results, Table 2 shows the measured drawdowns at the piezometer locations for both DE and HS. As can be seen from Table 2, the drawdown constraint given in Equation (8) is not satisfied for and 2 for both DE and HS. Although all the piezometers satisfy the given drawdown constraint for , big drawdowns are obtained in DE and HS especially at P2, P5, and P9. Since the pumping cost is directly related to the magnitude of the drawdowns, use of this solution may not be feasible in terms of the operation cost of the dewatering system. Note that these big drawdown values are significantly reduced by including an additional pumping well when . For this solution, only P1 and P5 have the drawdown values of 25.66 and 26.65 m while the others are around 20 m in DE. For the same solution with HS, the piezometers of P7 and P9 have the drawdown values of 23.76 m and 20.86 m while the other piezometers are around 20 m. The reason for not obtaining the same drawdown pattern for DE and HS when is the difference between the identified well locations and pumping rates. When the final drawdown distributions are compared, it can be seen that HS provided more uniform drawdown distribution than DE although its final objective function value (113.46) is slightly greater than that of DE (113.09). Finally, when , all the drawdown values become more uniform except for the piezometer P5 for both DE and HS. After evaluating the model results in terms of the final pumping rates and measured drawdown values, the solution of may be selected as the feasible solution for both DE and HS in terms of the engineering point of view. However, by considering well drilling, equipment, and operation costs, it is required to conduct a detailed analysis, which is beyond the scope of this study. For , variations of the identified well locations during the search process and the final locations of the pumping wells can be seen in Figure 5.

Piezometers . | Number of wells . | |||||||||
---|---|---|---|---|---|---|---|---|---|---|

1 . | 2 . | 3 . | 4 . | 5 . | ||||||

DE . | HS . | DE . | HS . | DE . | HS . | DE . | HS . | DE . | HS . | |

P1 | 4.49 | 4.49 | 12.65 | 12.63 | 20.00 | 20.00 | 25.66 | 20.00 | 20.00 | 20.22 |

P2 | 9.31 | 9.31 | 14.59 | 14.57 | 32.19 | 32.93 | 20.75 | 20.02 | 20.00 | 20.19 |

P3 | 9.31 | 9.31 | 14.24 | 14.23 | 20.00 | 20.00 | 20.00 | 20.11 | 20.00 | 20.62 |

P4 | 11.47 | 11.47 | 20.00 | 20.00 | 21.47 | 21.52 | 20.00 | 20.04 | 20.31 | 20.01 |

P5 | 23.53 | 23.53 | 20.06 | 20.04 | 38.46 | 37.88 | 26.65 | 20.03 | 28.70 | 24.67 |

P6 | 9.88 | 9.88 | 14.86 | 14.85 | 20.00 | 20.03 | 20.00 | 20.22 | 20.00 | 20.05 |

P7 | 7.09 | 7.09 | 15.73 | 15.73 | 20.13 | 20.14 | 20.12 | 23.76 | 20.00 | 20.07 |

P8 | 4.31 | 4.31 | 15.35 | 15.39 | 20.00 | 20.01 | 20.35 | 20.00 | 20.00 | 20.07 |

P9 | 4.79 | 4.79 | 20.00 | 20.00 | 32.02 | 31.57 | 20.00 | 20.86 | 20.00 | 20.69 |

Piezometers . | Number of wells . | |||||||||
---|---|---|---|---|---|---|---|---|---|---|

1 . | 2 . | 3 . | 4 . | 5 . | ||||||

DE . | HS . | DE . | HS . | DE . | HS . | DE . | HS . | DE . | HS . | |

P1 | 4.49 | 4.49 | 12.65 | 12.63 | 20.00 | 20.00 | 25.66 | 20.00 | 20.00 | 20.22 |

P2 | 9.31 | 9.31 | 14.59 | 14.57 | 32.19 | 32.93 | 20.75 | 20.02 | 20.00 | 20.19 |

P3 | 9.31 | 9.31 | 14.24 | 14.23 | 20.00 | 20.00 | 20.00 | 20.11 | 20.00 | 20.62 |

P4 | 11.47 | 11.47 | 20.00 | 20.00 | 21.47 | 21.52 | 20.00 | 20.04 | 20.31 | 20.01 |

P5 | 23.53 | 23.53 | 20.06 | 20.04 | 38.46 | 37.88 | 26.65 | 20.03 | 28.70 | 24.67 |

P6 | 9.88 | 9.88 | 14.86 | 14.85 | 20.00 | 20.03 | 20.00 | 20.22 | 20.00 | 20.05 |

P7 | 7.09 | 7.09 | 15.73 | 15.73 | 20.13 | 20.14 | 20.12 | 23.76 | 20.00 | 20.07 |

P8 | 4.31 | 4.31 | 15.35 | 15.39 | 20.00 | 20.01 | 20.35 | 20.00 | 20.00 | 20.07 |

P9 | 4.79 | 4.79 | 20.00 | 20.00 | 32.02 | 31.57 | 20.00 | 20.86 | 20.00 | 20.69 |

As can be seen, the final locations of the pumping wells are different in the DE and HS solutions although the final objective function values are very close to each other. For both DE and HS, the pumping wells are located around the borders of the excavation site, which is an expected behavior. For this dewatering configuration, Figures 6 and 7 show the cross-sections of A-A, B-B, and C-C given in Figure 3 to compare the groundwater levels for both DE and HS. As can be seen, the minimum drawdown condition is satisfied in all the piezometers and safe working conditions are satisfied by means of the determined dewatering system.

## CONCLUSIONS

In this study, a linked simulation–optimization model is proposed to solve excavation site dewatering problems. In the simulation part of the proposed model, the analytical solution of the groundwater flow equation is used to calculate the drawdown values at the piezometer locations. This simulation model is then integrated with two different optimization models where DE and HS optimization approaches are used. The objective of both optimization approaches is to determine the locations and pumping rates of the dewatering wells by considering the sum of all pumping rates as the objective function. Two main constraints are considered in the proposed model and both of them are included to the optimization process by means of the penalty function approach. After solving the problem for different well numbers, identified results indicated that both DE- and HS-based solution approaches can efficiently determine the dewatering system without requiring any trial-and-error process. Since both approaches considered the simulation of the groundwater flow process by using an analytical solution of the groundwater flow equation, it can only be used for problems with idealized hydrogeological conditions. However, this problem can be handled by integrating a more sophisticated numerical groundwater flow model, such as MODFLOW, to both DE- and HS-based optimization models. Furthermore, the objective function of the optimization model can be re-organized by considering well drilling, operation, and maintenance costs to make the proposed model more useful for practical cases.

Technical Bulletin of the Turkish Chamber of Civil Engineers (Denizli), 67