Multi-well capture zones in strip-shaped aquifers

Capture zone equations for a multi-well system in strip-shaped confined and unconfined aquifers with and without regional flow are presented. The aquifer is limited by two parallel boundaries that are either no flow (barrier) or inflow (variable head) so that aquifers with four possible boundary configurations are formed. The wellfield includes any number of extraction or injection wells or a combination of both types. The flow field in the strip-shaped aquifer was converted to its equivalent extensive aquifer using conformal mapping and image well methods. To delineate the capture envelope, the potential, streamline and stagnation point equations were derived using velocity potential theory. The solution permits rapid determination of the effect of number, position and extraction/injection rate of wells, boundary type and direction, and rate of regional flow on the size, shape and pattern of well capture zones. The derived equations are readily extended to water quality and quantity management simulations, as shown by embedding the equations within two optimization schemes, viz., Particle Swarm Optimization (PSO) and Genetic Algorithm (GA), to automatically determine the most efficient wellfield designs for pump-and-treat remediation, contaminant plume containment and pumping policy projects.


Introduction
Demands on water resources are increasing due to climate-induced changes and population growth [1][2][3], emphasizing the important role of managing the quantity and quality of water resources [4]. In pumped groundwater systems, knowledge of capture zones is a fundamental element of groundwater management.
The well capture zone within an aquifer is defined as the region from which water is withdrawn by one or more pumping wells [5]. After pumping is initiated, the capture zone grows with time and reaches its maximum size at steady state, which defines the capture envelope [6]. Capture zones are important for aquifer management, groundwater remediation projects [7][8][9][10], surface-groundwater interactions, well head protection [11], water rights, and in delineating impacts on transboundary aquifers. For over-exploited aquifers (where annual extraction exceeds recharge), identification of capture zones underpins optimal pumping plans to recover and sustain the depleted storage [12]. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Capture zone determination for confined aquifers is long established in groundwater engineering [13][14][15][16][17][18][19]. For instance, a two dimensional (2-D) pump-and-treat system was presented by Javandel and Tsang [18] including type curves to determine the number and pumping rate of wells to contain a contamination plume. Three-dimensional (3-D) capture zones are available for both horizontal drains and vertical wells in homogeneous, anisotropic aquifers [20]. Zlotnik [21] studied the effects of anisotropy on capture zones for partially penetrating wells. Analytical and semi-analytical expressions for multiple wells and extraction rates and locations based on complex potential theory and superposition were presented by Christ and Goltz [22]. Their solution was extended by Fienen et al. [23] to the computation of stagnation point locations in a multi-well system. Lu et al. [24] determined stagnation points in a flow field of a single pumping well located at the center of a recharge area.
The above-mentioned studies focused on simple aquifer geometries and did not account for different boundary configurations or unconfined flow. Efforts to extend analytical results include the case of a single well between two parallel streams [5], and multi-well systems in wedge-shaped [6] and peninsula-shaped confined and unconfined aquifers [25]. Capture zones for multi-well systems in rectangular bounded aquifers are also available [26], both for confined and unconfined aquifers. Such analytical solutions are used in optimization schemes. For example, Mantoglou [27] developed analytical models of saltwater intrusion in coastal aquifers of finite size pumped by several wells that were embedded within an optimization algorithm to find the optimal pumping strategy.
The approach most commonly used to obtain analytical results is the image well method [28] in conjunction with complex discharge potential theory. For complex boundary conditions, the number of imaginary wells can be large, and the method can become unwieldy. In this situation, however, conformal mapping can be used, since the number of imaginary wells for individual real single well is then limited to three or less for any aquifer boundary configuration.
Conformal mapping is commonly employed in areas of physics and engineering where Laplace's equation applies [29][30][31][32][33]. In groundwater flow problems, the modeling of seepage flow was solved by Fukuo and Kaihotsu [34], followed by other applications, e.g., Anderson [35][36][37][38]. Fitts [39] obtained an exact solution for two-dimensional flow to a well in an anisotropic domain. Chahar [40] used conformal mapping to provide an analytical solution to the problem of seepage from a soil channel with a curvilinear bottom. Klammler et al. [41] used the conformal mapping method to obtain the solution for flow fields near drain-and-gate reactive barriers. Using the same method, Klammler et al. [42] also developed an approximate analytical solution to the permeable reactive barrier (PRB) capture and release behaviour considering doubly-symmetric funnel-and-gate as well as drain-and-gate configurations. Similarly, Lu et al. [43] derived analytical solutions for pumping in a fully bounded rectangular aquifer. This approach was used by Strack [44] to obtain analytical results for flow to a well between two parallel rivers or rivers that are at a fixed angle. These results are generalized here.
Our purpose is to obtain a general solution for the capture zone of a multi-well system in bounded strip-shaped aquifers with or without uniform regional flow. The well system includes any number of arbitrarily located extraction or injection wells. The strip-shaped aquifer is located between parallel inflow (variable head) boundaries (surface water bodies like streams, rivers, canals, lakes, and seas) and no-flow (barrier) boundaries (bedrock, fault blocks, glacial tills or impervious rocks). Conformal mapping is applied to determine the complex discharge potential and streamline equations for determining the capture envelope. The method permits efficient and rapid modeling of capture envelopes, as shown by examples where the effects of number, position and extraction rate of wells and regional flow direction and rate on the capture envelope are investigated. We demonstrate the versatility of the approach by its application to groundwater quantity and quality management problems. Specifically, the capture zone model is embedded within a Particle Swarm Optimization (PSO) [45] and Genetic Algorithm (GA) [46] schemes to optimize design of pump-and-treat remediation, contaminant plume containment and pumping projects. Fig 1 shows a schematic plan view of different aquifer (confined or unconfined) configurations. Each case shows an aquifer with a fully penetrating well bounded by two parallel boundaries of infinite extent forming a strip-shaped aquifer. The aquifer is isotropic and homogeneous with uniform thickness. Steady state, 2-D flow is considered. In Fig 1A, the boundaries are fully penetrating streams (variable head boundary) having no hydraulic resistance with the aquifer. An extraction or injection well is located at (x w ,y w ). Note that the solution is for a multi-well system with arbitrary extraction/injection rates, and steady regional flow at an arbitrary angle (β) relative to positive x-axis. On physical grounds, the flow direction must be consistent with the hydraulic properties of the boundaries in that if the direction of regional flow is perpendicular to streams (β = π/2 or 3π/2), they are treated as constant head boundaries. In Fig 1B and 1C, one boundary is an inflow stream (variable head boundary) and another one is an impervious (no-flow) boundary. The aquifer depicted in Fig 1C is different from that of Fig 1B due to the position of wells relative to boundaries and vice versa. As a result, they should be treated as two separate well-aquifer systems. There are two no-flow boundaries in Fig 1D. In this case, on physical grounds, the regional flow direction must be parallel to the impermeable boundaries.

Capture zone
Using conformal mapping, the conceptual model of Fig 1 (z 1 −plane) is first mapped onto the upper half-plane in the z 2 −plane (Fig 2) using z 2 = exp(πz 1 /d) (see S1 File for details). When the boundaries are a combination of barrier-inflow and inflow-barrier, the z 1 −plane is transformed into the second auxiliary plane, the z 3 −plane, using the transformation z 3 ¼ z 2 1=2 . Transforming the zones in Fig 1A-1D  Transformation from the z 1 −plane to the z 3 −plane is given by: The complex potential for flow towards a single well at position z in the z 2 −plane with the regional uniform flow in an aquifer bounded by two parallel streams is [44]: where, Q w is the pumping or injection rate (L 3 T -1 ), q 0 is the uniform flow rate in the direction β (radians) relative to the positive x axis, ϕ 0 is the arbitrary potential in the absence of extraction and regional flow. Writing Eq (2) in the z 3 −plane gives: The parameter c is +1 for extraction wells and -1 for injection wells, respectively. Generalizing Eq (3) for the four boundary configurations of Fig 3 results in: where parameters J 1 to J 4 take values as given in Table 1. Because of the linearity of Laplace's equation, Eq (3) may be extended for a multi-well system using the principle of superposition: where N is the number of wells. The substitution of z 3 for z 1 gives: Using the following dimensionless quantities: Eq (6) becomes: with f 1Dj , f 2Dj , g 1Dj and g 2Dj as given in Appendix A. The real part of Eq (8) gives the dimensionless potential, ϕ D :

PLOS ONE
while the imaginary part gives the dimensionless stream function, ψ D : Plotting Eqs (9) and (10) results in a flow net that illustrates the capture curves of a multiwell system in a strip-shaped aquifer in dimensionless form. The flow net is conveniently plotted using the contour syntax in MATLAB. The stream function has a branch cut for each well that may be removed with the unwrap syntax in MATLAB.
In unconfined aquifers, the saturated thickness varies, so the potential equals to: Eqs (8-10) can be rewritten for unconfined aquifer by replacing the product Kb with Kh 0 (discharge per unit width) and using the following dimensionless terms:

Stagnation points
The capture envelopes that pass through the stagnation points (where the flow velocity is zero) in the flow field separate the regions of flow. To find the position of stagnation points, the derivative of the complex potential (Eq 5) is taken with respect to z 3 and the result set equal to zero: The roots of Eq (13) are the locations of stagnation points in the z 3 −plane, which are easily transformed to the z 1 −plane. Eq (13) was solved using MATLAB. The values of the stream function at the stagnation points are calculated (Eq 10) and the capture envelopes are drawn by the stream line passing through the stagnation points.

Drawdown
The steady state drawdown (s) in a confined aquifer from a fully penetrating extraction well is [47]: where h 0 and h are hydraulic heads at distance r 0 and r from the well, respectively. Since the

PLOS ONE
velocity potential is ϕ = Kh and the dimensionless head is h D = h/d, the dimensionless drawdown (s D ) can be written in terms of dimensionless potential as: where potentials ϕ 0D and ϕ D are calculated using Eq (9). The steady state drawdown (s) in an unconfined aquifer at a distance (r) from a fully penetrated extraction well with radius (r w ) is [47]: Using Eq (11) the dimensionless head in an unconfined is written as: From Eqs (16) and (17), the dimensionless drawdown in an unconfined aquifer is: Eqs (15) and (18) can be used to demonstrate the application of proposed model to water quantity management (see Water quantity management projects Section below).

Results and Discussion
In the following sections, Eqs (9), (10) and (13) are solved for strip-shaped aquifers with different boundary configurations and flow conditions, after which we plot the capture envelopes of the well(s). First, the effects of various boundary configurations (Fig 1) on the shape and properties of the capture zones are investigated, then the effect of extraction rates and regional flow direction and rate are explored. Due to numerous parameters (number, type, location and extraction/injection rates of wells, and the rate and direction of regional flow) that influence the flow field, many capture envelopes can be simulated. Here, representative examples of capture zones in each parallel boundary configuration are presented, along with the water quantity management and optimized remediation schemes for a pump-and-treat method. The dimensionless values of parameters (i.e., x wD , y wD , Q wD ) can be calculated using Eq (6). https://doi.org/10.1371/journal.pone.0229767.t002

Effects of boundary configuration on the capture zones
Capture zone of well(s) in a strip-shaped aquifer with inflow-inflow boundaries. In this example, we consider a strip-shaped aquifer bounded with inflow-inflow boundaries ( Fig  1A) in which the direction β and the rate q 0D of the uniform regional flow are 0 rad and 0.001, respectively. The parameter J is defined in Table 1. Eqs (9), (10), and (13) are solved for 1, 2, 3 and 5 wells (the position and pumping rate of extraction wells are defined in Table 2). Fig 4A-4D illustrate, respectively for the four cases, their potential, streamlines, stagnation points and capture envelopes. The black diamond shows the well position and the green diamond illustrates the position of the stagnation point. The dashed lines indicate the velocity potential and the solid lines represent the streamline reaching the location of the extraction well. In Fig 4A, the well extraction rate is provided from regional flow and two inflow boundaries and the capture envelope extends symmetrically. Due to the direction of regional flow, the capture envelope trends toward the west (left) side. In Fig 4B, well (2) extracts water from the streams as well as the uniform regional flow. A sharp groundwater divide, effectively a barrier, is established between two wells so that the uniform regional flow provides no water to well (1). Streamlines run parallel along the divide. In Fig 4C, well (3) has the largest capture envelope because it has the maximum extraction rate. The three wells gain water from both boundaries. Two groundwater divides are also formed, which separates the capture zone of well (2) from those of the other wells. In Fig 4D, each well has its own capture envelope. Wells (1) and (2) gain water from the north and south side boundaries, respectively, whereas well (3) is recharged by the two boundaries and regional flow. Its capture zone shape shows the influence of the other wells. Well (4) is fully supplied by the northern boundary. A groundwater divide separates the capture zone of well (5) from those of the others while both boundaries supply water to it. In both Fig 4C and 4D, a stagnation point exists for wells (2) and (5) Figure B in S1 File.
Capture zone of well(s) in a strip-shaped aquifer with inflow-barrier boundary conditions. For this boundary configuration (Fig 1B), the aquifer thickness, hydraulic conductivity and the direction and regional flow rate are the same as the previous examples, as specified in Table 2 with J values as given in Table 1. Eqs (9), (10) and (13) are solved for 1, 2, 3 and 5 wells and the resulting capture curves are plotted in Fig 5. In Fig 5A, the case of 1 well, the extracted water is provided from both the regional flow and the inflow boundary. The streamlines are parallel to the no-flow boundary as expected. In Fig 5B, well (1) is mainly supplied by the stream, while well (2) gains water from regional flow and the stream. In Fig 5C, well (3) has the largest extraction rate, with water supplied by the stream and the regional flow. A small portion of regional flow reaches well (1). The capture zone of well (1) extends towards the stream as well as the west and east sides of the aquifer. The capture envelope of well (2) is smaller, due to its short distance from the north boundary and lower extraction rate. Fig 5D demonstrates the competition among wells in gaining water from the north boundary based on their pumping rate and distance from the boundaries. Only well (2) benefits from regional flow and its capture zone is limited by the no-flow boundary and manages to gain a small portion of water from the stream along three narrow zones: a) between the capture envelopes of wells (3) and (5), b) the capture envelopes of wells (1) and (3) and c) the capture zone of well (5) and the noflow boundary.
Capture zone of well(s) in a strip-shaped aquifer with barrier-inflow boundaries. Fig 6  shows results computed using Eqs (9), (10) and (13), solved for the configuration shown in Fig  1C, with parameters as given in Tables 1 and 2. The aquifer thickness, hydraulic conductivity and direction and rate of regional flow are the same as in the previous examples. Fig 6A depicts  that well (1) is supplied by the stream as well as regional flow. In Fig 6B, wells (1) and (2) gain water from the stream and well (2) captures the main portion of regional flow and allows a smaller portion to reach well (1). In Fig 6C, the stagnation points of wells (2) and (3) are positioned on the barrier boundary and form a small stagnation zone. Well (3) prevents regional flow to reach wells (1) and (2) and encompasses the capture envelope of well (1). As a result, wells (2) and (1) are only fed by the inflow boundary. Fig 6D demonstrates that well (1) is recharged by regional flow only and wells (2) and (5) by the stream. Water reaches well (4) via a narrow pathway connecting the well to the upstream regional flow and via two sections of the inflow boundary. Well (3) is fed by the stream and its capture envelope encompasses that of well (2).

Capture zone of well(s) in a strip-shaped aquifer with barrier-barrier boundaries.
The aquifer in this example corresponds to the case in previous section. The difference is that the in-flow boundary is replaced by barriers (Fig 1D). According to Table 1, for this boundary configuration the value of parameters J 1 , J 2 , J 3 and J 4 are all equal to +1. As above, the regional flow direction is from west to east (β = 0). Eqs (9), (10) and (13) are solved for 1, 2, 3 and 5 well systems-potential and streamline functions are plotted in Fig 7. As shown in the figure, the streamlines are parallel to the barrier boundary before merging to the wells and equipotential lines are perpendicular to the barrier boundaries, as expected. In Fig 7A, the extraction rate is provided by regional flow from the west side of the aquifer and also from the infinite boundary in the east. As a result, the west and east flowing fronts intersect and two stagnation points form on the boundaries. In Fig 7B, regional flow is fully captured by well (2) and well (1) gains water from the east infinite side. In Fig 7C, the capture envelopes of well (3) and well (2) extend to the west and east sides, respectively, while that of well (1) extends to both sides along two narrow channels parallel to the barrier boundary. In Fig 7D wells (1) and (2) capture water from the west side of the aquifer and wells (4) and (5) from the east side while well (3) captures water from both sides. Well (3) forms two stagnation points each located on one barrier boundary. Comparing this figure with Figs 4-6, one can see how boundary types, pumping rates, and well(s) position control the shape, size and pattern of capture envelopes. Wells interact differently with each other and with the boundaries. In a regional pumping scheme, these differences are important for water quantity management, allocation of water rights, and for pumping rate permits.

Effects of regional flow direction and rate on the capture zones
In the previous examples, we assumed that the direction of regional flow was from west to east (β = 0). To demonstrate the effect of regional flow direction on the capture envelopes, Figs 4D, 5D, 6D and 7D are replotted for β = π in Fig 8 while Fig 4B and 4D are replotted for β = π/2 in Figure C in S1 File. The difference due to the flow direction is profound, as readily seen in the shape, size and pattern of capture envelopes and the interaction of wells and boundaries. In Fig 8A, in contrast to its corresponding Fig 4D, well (5) fully captures and prevents regional flow to reach other wells. As the result, it is fed along a shorter segment of the north boundary; the capture zone of well (1) leans somewhat to the west and its stagnation point is displaced further south; well (2) is only fed by the southern stream; well (3) is fed along a longer section of the northern stream and a shorter segment of the southern stream while the capture envelope of well (4) slopes further to the east. In Fig 8B, contrary to Fig 5D, regional flow is mainly captured by well (5) although and a small portion of it is captured by well (2) through a narrow flow channel along the no-flow boundary. In Fig 8C, in comparison to Fig 6D, the pattern of capture envelopes changes considerably and regional flow is captured by wells (4) and (3) instead; well (1) reaches the stream and gain water from it; the capture envelope of well (4) encompasses that of well (5), which moves further to the east. In Fig 8D,   greatly change in size and shape and pattern; the capture envelope of well (1) becomes larger than that of well (3), which extends to the east by encompassing that of well (5); well (2) is only supplied by the southern boundary. The dimensionless regional flow is increased to 0.02 compared to that of the cases presented in Figs 4-8 (which was 0.001), with the new results given in Figures D-H in S1 File. All the capture envelopes in these figures extend in the opposite sense to the regional flow direction, i.e., the regional flow contributes to the extraction rate of all wells. It is clear that the size, shape, dimensions and pattern of well capture zones and the interaction of boundaries and wells are greatly changed compared to Figs 4-8.
Figs 4-8 and Figures C-H in S1 File demonstrate that the capture zone model developed above can help delineate the interaction between surface-subsurface flows in water resources management projects such as allocating water rights, issuing well drilling permits, pumping rate permits, etc., and for the sustainable development of water resources in a region.

Validation of the solution
To validate the developed capture zone solutions, the capture zones of 5 wells of Fig 8 were generated by MODFLOW 2000 [48] and MODPATH [49] and plotted in Fig 9. For the four plots shown in Fig 9, the well coordinates as well as their extraction/injection rates were set to those of Fig 8. The capture zones and the pattern of streamlines generated by the numerical models are comparable to the analytical model capture curves (i.e., Fig 8). The slight difference is mainly due to the approximation inherent in the numerical models. It is worth mentioning that the usage of the capture envelope equations developed in this paper is far easier, less time consuming, and more accurate than using the abovementioned numerical models, although numerical models are more flexible in dealing with complicated boundary conditions [10]. Figs 8 and 9 demonstrate the potential of developed analytical capture zone models to verify the accuracy of numerical models.

Effects of the well distance from boundaries on the capture zones
The positions of the wells are kept constant relative to the south side stream as in Fig 4D, and the northern boundary is moved away so that its distance to the wells increases by a factor of 2 ( Fig 10A) or 4 (Fig 10B). Fig 10A illustrates that the capture zones are reshaped such that well (1) no longer is fed by the north side stream, but is instead fully fed by regional flow, while well (3) gains less water from the north side stream. Well (4), which was fully supplied with the north side stream, now gains water from both streams. Well (5), which was fed by both streams, is now recharged by the southern stream only. In Fig 10B, the distance between wells and the northern stream is increased, so that the capture zone of wells (1), (3) and (5) do not reach it, while well (1) is fed only by regional flow. Wells (3) and (5) are no longer supplied by the northern stream. The contribution of the north side stream to well (4) is also negligible. Practically speaking, for this case, the strip-shaped aquifer acts as an aquifer located within a semi-infinite spatial domain.

Groundwater remediation scheme design
In this section, use of dimensionless capture type curves in water quality management such as groundwater remediation projects is demonstrated. A variety of remediation technologies can be used to ensure contaminant plumes are either contained, removed and/or treated so they no longer pose a threat to water resources. In a typical in situ remediation project such as bioremediation, pump-and-treat, plume containment, etc., polluted groundwater is extracted using one or more wells, treated and possibly reinjected [50]. In such projects the main

PLOS ONE
elements of an efficient and cost-effective design are determination of the capture zone, the optimal number of wells, the optimal injection/extraction rates, and the layout of wells with respect to the plume [7-10, 18, 26, 50-52].
We consider a confined aquifer with the boundary configuration of type (a) in Fig 1. The problem setup includes known dimensions and hydraulic properties (aquifer thickness, length and width are 20 m, 1000 m and 500 m, respectively; K = 5 m/d and q 0 = 0.1 m 2 /d). A contaminant plume is located along the direction of regional groundwater flow (from the west to the east, β = 0). The aim is to contain the plume hydraulically and prevent its extension downgradient using two or three wells. A logical step is to position an injection well (to inject treated water) in the vicinity of the contaminant source and to place one or more extraction wells at the leading edge of the plume. A solution could be determined manually by using Eqs (9) and (10) to generate a set of capture zone curves for various well layouts and extraction/injection rates. This approach, although cumbersome, could achieve a satisfactory plume capture design.
Two automated parameter optimization schemes are used instead, i.e., Eqs (9) and (10) are embedded into the PSO algorithm (Appendix B) and GA (Appendix C). Fig 10 illustrates the optimal capture envelopes for two-and three-well remediation schemes, with the corresponding well positions and extraction/injection rates given in Table 3. As is evident in Fig 11, both systems involve a closed region that contains and encompasses the plume. Table 3 shows that both optimization schemes result in similar optimal solutions. This approach can easily be extended to include financial constraints associated with well installation and operation. In Fig  11, the three-well solution has a lower total injection/extraction rate, so it represents the optimal solution if the remediation scheme is dominated by pumping costs. Fig 10 was also generated by the abovementioned numerical models and presented in Figure I in S1 File. In this simulation, if the well positions and flow rates were unknown, finding the optimum plan for the plume containment would be difficult to determine manually.

Water quantity management projects
The availability and sustainability of groundwater in many major aquifers in the world are threatened due to excessive depletion. To avoid this, a typical management strategy is to allocate well pumping rates in such a way that the drawdown in the aquifer does not exceed a certain value (permissible drawdown). As an example, we consider a confined aquifer with the boundary configuration of type (a) in Fig 1. The problem setup includes known dimensions and hydraulic properties (aquifer thickness, length and width are 20 m, 1000 m and 500 m, respectively; K = 5 m/d and q 0 = 0.1 m 2 /d). Further, it is assumed that the aquifer is pumped by five wells with positions given in Table 4. The objective of this example is to calculate extraction rates such that the maximum drawdown is less than 4 m (Fig 12). Table 3. Extraction/injection rates and well positions in two remediation scenarios.

Pump-and-treat remediation scenarios
Type of well PSO GA x wD y wD Q wD x wD y wD Q wD

PLOS ONE
Based on Eq (7), the wells and the aquifer parameter values are converted to dimensionless values and Eq (15) is solved using the abovementioned optimization algorithms. The calculated optimum pumping rates by both optimization algorithms are presented in Table 4 with similar results. The capture envelope patterns of the 5 wells for the optimum solution are presented in Fig 11A and the corresponding dimensional drawdown contours calculated by Eq (15) are shown in Fig 11B. Note that the permissible 4-m drawdown is not exceeded. Table 5 presents the definition of variables in the equations.

Conclusions
Analytical models of the capture zone (stream function, velocity potential and stagnation point) for a multi-well system with any number of extraction/injection wells in confined and unconfined strip-shaped aquifers, with/without uniform regional flow are developed from the theory of complex velocity potentials and the Schwarz-Christoffel conformal mapping. The solution is provided for four strip boundary configurations: barrier-barrier, barrier-inflow, inflow-barrier and inflow-inflow. The solution, which is straightforward to compute, can be used to generate dimensionless capture type curves. The applicability of the developed model to evaluate surface-subsurface water interaction, well-head protection plans and water rights adjudication-which are key factors in the management of water resources-were demonstrated through numerical examples (Figs 4-8 and 12). These examples also displayed that the size and pattern of the capture envelopes are controlled by the extraction/injection rate, distance of wells from the boundaries, distances between wells, type and number of wells, the rate and Fig 11. Remediation of groundwater by the pump and treat method a) two-well and b) three-well system. "I" shows the injection wells and "E" the extraction wells. The colored area is the original contaminant plume. The thick black line surrounding the plume is the dividing streamline that defines the capture zone that separates the plume from the rest of the aquifer and β = 0.
https://doi.org/10.1371/journal.pone.0229767.g011 Table 4. Well positions and the calculated extraction rates of wells in the quantity management scenario.

Well number
Well coordinates PSO GA PLOS ONE direction of regional flow and the initial head along the boundaries. Further, the derived solution was embedded with two optimization schemes (PSO and GA) for the determination of optimal plume capture schemes for pump-and-treat remediation projects and for determination of the optimal pumping policy for sustainable water extraction schemes, with similar results obtained by both methods. The solution was validated by numerical models, and is thus suitable for quantifying the accuracy of such models for complex water injection/extraction scenarios. In general, the proposed models can be used by hydrogeologists, engineers and planners to design optimal layouts of groundwater remediation projects (e.g., in situ bioremediation or chemical oxidation schemes) that involve pump-and-treat and plume containment operations, to design the most cost-effective pumping schemes for sustainable development of water resources, to study the interaction of surface and groundwater resources and to verify numerical models.    (B-6)where c 1 and c 2 are acceleration constants (learning factors), r 1 and r 2 are random numbers uniformly distributed between 0 and 1, subscripts p and g respectively indicate the pbest and gbest values of the particles, ω is the inertia weight used to control the impact of the previous history of velocities on the current value, χ is the constriction coefficient, applied to restrain the velocity (v), number of iterations (t), well rates (Q w ) and the particle position (z w ).

5.
Repeat steps 3 and 4 until the fitness criterion or the maximum number of iterations is reached.
The above procedure was performed within MATLAB.
For water quantity management projects, the aim is to optimize extraction rates contingent upon a certain value of drawdown (permissible drawdown, s p ) where the objective function and constraints are: Maximize Q wjD ðB À 7Þ (B-7)subject to: Q jmin � Q wjD � Q jmax s max � s p ðB À 8Þ (B-8)where, Q wjD is the extraction rates based on Eq (8), Q jmin and Q jmax are userdefined minimum and maximum extraction rates for the j th well.
Supporting information S1 File. Mapping the conceptual model from the physical plane z 1 to the imaginary planes z 2 and z 3 , and effect of regional flow rate and direction on the shape and orientation of well capture zones. (DOCX)