The Method of Images Revisited: Approximate Solutions in Wedge‐Shaped Aquifers of Arbitrary Angle

This paper focuses on deriving new approximate analytical solutions in wedge‐shaped aquifers. The proposed methodology is applicable to any type of aquifer namely, leaky, confined and unconfined, under both steady state and transient flow conditions. By applying the method of images and separating the flow field into sections using physical arguments, approximate analytical expressions are obtained for the drawdown function, which in contrast to the conventional theory, are applicable to any arbitrary wedge angle. Moreover, the solutions fully observe the boundary conditions, while they preserve the continuity of the drawdown, which can be calculated directly at any point of the flow field. Nevertheless, comparison of the results of the new approximate analytical solutions to numerical ones, has been considered necessary to check their validity. MODFLOW, a well‐known numerical tool is used to calculate the numerical results. The discrepancies between the numerical results and those of the approximate analytical solution are negligible. The main advantages of the proposed methodology are the following: (a) The model needs only finite number of terms compared to conventional analytical and numerical solutions that involve infinite series, (b) The computational load is low, so it can be easily used in conjunction with meta‐heuristic algorithms to solve groundwater resources optimization problems, (c) Stream depletion rate can be calculated rather accurately and (d) The method is applicable to related flow problems.


Introduction
Aquifers bounded by two boundaries, either constant head such as streams and lakes or no flux, such as impermeable rocks, intersecting at angle smaller than 90°are called wedge-shaped (Mahdavi, 2021).Their study has attracted increasing research interest over the years.Most of the research papers have focused on analytical, approximate analytical and semi-analytical expressions for drawdown calculation and estimation of the aquifer's parameters (Singh, 2001;Yeh et al., 2008).Analytical procedures have been followed in several research papers to obtain expressions for the drawdown function (Zong et al., 2022).The well-known Hankel transform has been used by Chan et al. (1978), Yeh and Chang (2006), and Chuang and Yeh (2018) to obtain analytical solutions in wedge-shaped aquifers under steady state and transient flow conditions respectively.Chen et al. (2009) have applied the method of images to describe the aquifer's response to a constant-rate pumping well.Other method used to derive analytical solution to wedgeshaped aquifers is the Laplace transform (Lin et al., 2018).Recently, more complicated aquifer's shapes, such as triangle-shaped, annular wedge-shaped, trapezoidal-shaped, have been studied analytically (Asadi-Aghbolaghi & Seyyedian, 2010;Kacimov et al., 2017;Leray et al., 2019;Mahdavi, 2019Mahdavi, , 2022;;Mahdavi & Yazdani, 2021;Nagheli et al., 2020;Zlotnik et al., 2015).
Regarding wedge-shaped aquifers, in particular, dimensionless type curves of flux-time and drawdown-time are given for homogeneous aquifers by Sedghi et al. (2010Sedghi et al. ( , 2012) ) as well as for heterogeneous ones by Samani and Sedghi (2015), using integral transform methods.Moreover, estimation of hydraulic parameters and prediction of the discharge of a qanat in alluvial aquifers is achieved via the semi-analytical approach introduced by Sedghi and Zhan (2022).Approximate analytical solutions are simplified expressions aiming to describe complex problems with good accuracy (Moutsopoulos et al., 2022).Expressions obtained from approximate procedures have been combined very efficiently with meta-heuristic methods in groundwater optimization problems (Christelis et al., 2019;Karpouzos & Katsifarakis, 2021;Mallios et al., 2022;Rodriguez-Pretelin & Nowak, 2019).Their main asset is their low computational load, in comparison with numerical solutions.Low computational volume is very crucial when meta-heuristic methods are used, since the simulation model should be executed hundreds or thousands of times, in order not to compromise the optimization process.Further discussion about approximate analytical solutions will follow in Section 3.
In this framework, approximate analytical solutions for wedge-shaped aquifers are sought, based on the image well theory.Observance of both boundary conditions, either constant head or no flux, as well as continuity of the drawdown function were set as prerequisites.The concept of the proposed methodology is the division of the real flow field into two sections, where different fictitious wells are taken into account.The proper division of the real flow field in sections, the required number of image wells, together with their location and type (pumping or injection), are discussed in the following sections.The respective geometrical proof is presented in Appendix A.

Outline of the Method of Images
The basic concept of the method of images is that a boundary can be "removed" by adding a number of fictitious (or image) wells, symmetrical of the real ones with respect to it, resulting into an equivalent infinite flow field (Bear, 1979;Haitjema, 2006;Mahdavi, 2020;Nikoletos, 2020).The sign of the flow rate of each image well depends on the boundary condition and guarantees its observance (Katsifarakis et al., 2018;Samani & Zarei-Doudeji, 2012).From the mathematical point of view, it is a specific application of the Green's function and is applicable to problems described by the Poisson equation (Mohamed & Paleologos, 2018).Its use is extensive in many scientific fields such as groundwater hydraulics (Dewandel et al., 2022;Kuo et al., 1994), electrostatics (Nguyen & Mehrabian, 2021), magnetics (Curtis et al., 2015) and optics.The method of images has been widely used in groundwater flow simulation problems to calculate hydraulic head level drawdown (Atangana, 2014;Nikoletos & Katsifarakis, 2022;Penny et al., 2020), to describe interaction between ground and surface water (Anderson, 2003) and to optimize the management of aquifers (Katsifarakis, 2008) and especially coastal aquifers facing saltwater intrusion problems (Etsias et al., 2021;Mantoglou, 2003).It is worth mentioning that the method of images gives exact solutions in wedge-shaped aquifers bounded by two boundaries intersecting at angles of: 90°, 60°, 45°, 30°etc.(each angle verifying Equation 1) where θ, is the boundary intersection angle and N, the number of fictitious wells (Ferris et al., 1962).

Approximate Analytical Solutions
Due to the complexity of many flow fields, exact solutions cannot be found.In such cases, numerical methods are usually applied, as long as the required computational resources are available.Approximate analytical solutions could be a good alternative, if the computational volume is low, and the introduced error is acceptable.As solutions produced by numerical methods are inherently approximate, too, convergence of the results of approximate analytical and numerical methods point out the validity of the proposed solutions.

Basic Concept of the Proposed Solutions
The aim of the proposed solutions is to calculate the drawdown distribution in wedge-shaped flow fields, with satisfactory accuracy, while observing the boundary conditions.Following the approach developed by Nikoletos and Katsifarakis (2022), we divide the real flow field into two sections.In each section a number of fictitious wells is used in a way that observance of the boundary conditions is achieved.The proposed division of the flow field in two sections, does not disrupt continuity of the drawdown, but the flow velocity field is discontinuous, along the straight line that separates the field sections.

Water Resources Research
10.1029/2022WR034347 The main assets of the proposed solutions are: (a) The low compuational volume, which facilitates their combination with meta-heuristic optimization methods and (b) The full observance of the boundary conditions, which facilitates calculation of stream depletion rates.The accuracy of the results as well as the applicability range of the proposed approximate solutions are discussed in the following sections.Kuo et al. (1994) proposed a new approach for predicting drawdown in aquifer with irregularly shaped boundaries using the image well method.More recently, a novel methodology for estimation of stream filtration from a meandering stream, is introduced by Huang and Yeh (2015) by applying the image well theory.In both papers the decision variables were the flow rates of the image wells which are determined by solving a system of equations.

Comparison With Previous Studies
In this paper the decision variables are the number of fictitious wells and their locations, taking into account the physical interpretation of the method.The proposed methodology is more suitable to problems where the intersecting boundaries are straight lines while the aforementioned methods are preferred where the boundaries are curved.It is worth mentioning that combination of the proposed approach with the method introduced by Kuo et al. (1994) or Huang and Yeh (2015) could lead to more accurate results for the calculation of the drawdown.This will be an issue of future study.

Structure of the Drawdown Function
The proposed methodology can be applied to any type of aquifer under steady state or transient flow conditions.
The drawdown functions at any point, due to steady-state pumping from a single well, in infinite aquifers, confined, leaky and unconfined, are given by the following relationships respectively (Hantush, 1967;Jacob, 1946;Thiem, 1906): ) In the above equations, Q w is the flow rate of the well, r is the distance between the pumping well and the observation point, K is the hydraulic conductivity, T is the aquifer's transmissivity, R is the radius of influence, is the leaky coefficient and c the hydraulic resistance of the aquitard, K o is the modified Bessel function of second kind and zero order and H and H 1 are the distances between a reference level and the final and initial water table level, respectively.
To demonstrate the procedure, we consider steady flow in confined aquifers.

Intersecting Angle θ Between 90°and 60°I
n the following paragraphs, we present solutions for wedge-shaped aquifers, bounded either by two constant head or by two impermeable boundaries.We present the proposed methodology for constant head boundaries.The process is the same in case of impermeable boundaries, except the kind of the fictitious wells.
The exact solutions for the drawdown due to a pumping well in flow fields bounded by two boundaries intersecting at angles 90°and 60°make use of 3 and 5 fictitious wells, respectively.Based on that, we postulate that the proposed solutions should make use of 3-5 fictitious wells.As shown in Figure 2a, W i 1 and W i 2 are the images of the real well with respect to the closest and the furthest boundary respectively; we call them first-order images.W i 3 and W i 4 are the images of W i 1 and W i 2 with respect to the other boundary; we call them second-order images.W i 5 and W i 6 , resulting in a similar way, are the third-order images, and so on.The odd order images of a pumped well are injection wells, while the even order ones represent pumping wells.
From the physical point of view, the influence of the image wells should decrease, as their order increases.For instance, the first order images should affect more the conditions in the real flow field, namely they should be Water Resources Research 10.1029/2022WR034347 located closer to it than the second order ones, the second order images should be closer to the real flow field than the third order ones and so on.The largest order image wells should operate alternately in parts of the flow field, so that both boundary conditions are satisfied.According to the geometrical proof of Appendix A, if the angle θ is larger than 72°, the third order images are closer to the real flow field than the second order ones.Consequently, we conclude that for angles larger than 72°, 3 fictitious wells should be used, while for angles smaller than 72°, 5 fictitious wells should be used, provided that the conditions described by inequalities Equations A3 and A4 of Appendix A are also satisfied.Otherwise the largest order images should not be used.
From the physical point of view, we expect that the drawdown at any point of the real flow field decreases as the angle between the boundaries decreases.According to Figure 1, the pair of W r and W i 1 results into conditions equivalent to those caused by one constant head boundary.The second constant head boundary should further reduce the drawdown.Therefore the system of W i 2 , W i 3 , W i 4 , and W i 5 or W i 6 (for the pair of wells W r and W i 2 , respectively) should represent the influence of the second boundary.This reduction depends on the distance of the above mentioned system from the real flow field.As the angle θ increases, the system of the image wells diverges from the real flow field, so its influence decreases.
Therefore, the following two-branch functions could describe the drawdown distribution to the real flow field in each case. (3) In the above equations, r i (i = 1-6) are the distances between the observation point, namely "m" and the fictitious wells.The real and the fictitious wells for the aforementioned ranges of θ values are shown in Figures 2a and 2b, respectively.Water Resources Research 10.1029/2022WR034347

Intersecting Angle θ Between 60°and
f the intersecting angle of the boundaries is 60°, the exact analytical solution makes use of 5 fictitious wells.If the angle is 45°, 7 fictitious wells are placed to the equivalent infinite flow field.We postulate, then, that the proposed approximate solutions should make use of 5-7 fictitious wells.The main difference compared to the case 90°-60°i s that the largest order images are pumping wells.Taking into account the physical interpretation of the aquifer's response, the 4th order images should be located farther than the 3rd order ones from the real flow field.
According to the geometrical proof of Appendix A, if the angle is larger than 51.42°, the 4th order images are closer to the real flow field than the 3rd order ones.The largest order image wells should be used alternately in parts of the flow field, in order to ensure the observance of boundary conditions.Consequently, we conclude that for angles larger than 51.42°we should use 5 fictitious wells, while for angles smaller than 51.42°, 7 fictitious wells should be used, provided that the conditions described by inequalities Equations A3 and A4 of Appendix A are also satisfied.Otherwise the largest order images should not be used.
Verification of the drawdown reduction as the angle θ decreases is needed.We examine separately W r and W i 1 (or W r and W i 2 ) and the rest of the wells.According to Section 3.4.1 the system of wells W i 2 , W i 3 , W i 4 , W i 5 , or W i 6 (for the pair of wells W r and W i 2 ) results in additional water inflow to the real flow field.Consequently, it is proved that W i 5 and W i 6 are closer to the real flow field than W i 7 and W i 8 , respectively.According to the geometrical proof of Appendix A, the aforementioned conditions holds.
Hence, the following two-branch functions describe the drawdown distribution to the real flow field in each case.
For 60°> θ > 51.42°s It is worth mentioning that in this case the number of fictitious wells is the same as in the case of 72°> θ > 60°.
For 51.42°> θ > 45°s In the above equations, r i (i = 1-8) are the distances between the observation point, namely "m" and the fictitious wells.The real and the fictitious wells for the aforementioned ranges of θ values are shown in Figures 2b and 2c, respectively.

Intersecting Angle θ < 45°F
or the rest of the cases, we follow exactly the same methodology as described in Sections 3.4.1 and 3.4.2,for consecutive analytical solutions.The determination of the critical angle of each case as well as the operating wells of each section are based on Equation A6 and Figure A1, found in Appendix A.

Results and Discussion
MODFLOW (Harbaugh, 2005;Harbaugh et al., 2017), an established, finite-difference, groundwater flow simulation model is used for the evaluation of the proposed solutions.It is widely used by hydrogeologists to simulate real and theoretical flow fields for estimation of the response of the aquifers (Malekzadeh et al., 2019).The program has been applied extensively to aquifers bounded by two or more boundaries of irregular shape (Aghlmand & Abbasi, 2019;Karimi et al., 2019).

10.1029/2022WR034347
Here we use MODFLOW to evaluate the quality of proposed approximate solutions.Convergence of the results of approximate solutions and numerical methods is a strong indication that the solutions are valid.We consider six cases of wedge-shaped aquifers, bounded by two constant head boundaries, intersecting at the point O (0, 2,000).One boundary is described by the equation x = 0.A well with coordinates (320, 1,500) pumps at a flow rate Q W = 0.02 m 3 /s.The hydraulic parameters of the well and the aquifer are listed below.We consider steady-state flow conditions to facilitate the comparison with the approximate analytical solutions.
A grid of 5 × 5 m has been used to run MODFLOW and the results have been compared to those obtained by the approximate analytical solutions.Visualization of the results has been made through ModelMuse (Winston, 2009(Winston, , 2020)), a well-known graphical user interface for groundwater simulation models.The results for 6 θ values are discussed in the following paragraphs.
For θ = 90°, 60°, and 45°the solutions, obtained by the image well theory, are exact.These cases serve rather to check MODFLOW results.The analytical solutions give slightly larger values than MODFLOW.Discrepancies are smaller than 5% at all points of the real flow field except for the location of the well.
For θ = 80°, 65°, and 55°the solutions are approximate.The approximate analytical solutions render larger drawdown values than MODFLOW.Discrepancies are smaller than 6% at all points of the real flow field except for the location of the well.To further ensure the validity of the proposed approach, comparison between the drawdown values obtained by the two methods at an observation point with coordinates (400, 1,500), for different intersection angles, are summarized in Figure 4.

Conclusions
New functions for the calculation of the drawdown distribution of wedge-shaped aquifers via approximate analytical procedures are presented, based on the method of images.The number of fictitious wells that should be used is small and depends on the intersection angle.Observance of the boundary conditions is achieved through the use of two-branch functions.The largest order fictitious wells are activated alternately in parts of the flow field.The division of the flow field in two sections does not disrupt continuity of the drawdown.The only drawback is that the flow velocity field is discontinuous, along the straight line that separates the field sections.
To check the validity of the approximate solutions, results for six application examples have been compared to numerical ones, obtained by MODFLOW.For θ equal to 90°, 60°, and 45°, namely when the method of images is exact, the discrepancies are generally smaller than 5%, while in the other cases, discrepancies are generally smaller than 6%.So, probably, the error introduced by the approximate analytical solutions, described in the previous sections, is smaller than the error introduced by numerical tools, using reasonable grid dimension.The main attractive features of the proposed solutions are the following (a) They make use of small number of terms compared to conventional analytical and numerical solutions that involve infinite series, (b) the computational load is low, so they can be easily used in conjunction with meta-heuristic algorithms to solve groundwater resources optimization problems such as minimization of the impact of stream depletion induced by pumping on aquatic ecosystems in rivers, (c) stream depletion rate can be calculated rather accurately, and (d) the methodology is applicable to related flow or transport problems.Therefore, based on the proposed methodology, introduction of new analytical models for the estimation of stream depletion rate will be an issue of future study.
Let α and b be the angles of OW r with the field boundaries Ox 1 and Ox 2 respectively.Let angle c be given by the following equation: It follows, from Equations A1 and 1, that If b < c, then (Ox 2 , OW i 2 ) < c, namely W i 2 lies on same side of Ox n , with W r .Consequently, W i 4 , the image well of W i 2 with respect to Ox 1 lies under the side of Ox n+1 .Since Ox 2 is the bisector of Ox 1 ′ and Ox n+1 ′, the image well of W i 4 with respect to Ox 2 , namely W i 6 , lies on the opposite site of Ox 1 .The same condition holds for the consecutive mirror wells with respect to the corresponding boundary until the penultimate well, namely W i (k 2) .For the last order image well, Ox 1 is the perpendicular bisector of W i (k 2) and W i k .Therefore W i (k 2) and W i k lie on the opposite site of Ox 1 .As shown in Figure A1

Figure 1 .
Figure 1.Separation of the wells into two groups.Note: BC stands for boundary condition.

Figure 3 .
Figure 3.Comparison of equipotential lines for different intersection angles obtained by approximate solutions and MODFLOW.

Figure A1 .
Figure A1.Geometrical relationships between real and fictitious wells for arbitrary wedge angle.