Analytical Solution for 2D Inter-Well Porous Flow in a Rectangular Reservoir

: Inter-well ﬂuid ﬂows through porous media are commonly encountered in the production of groundwater, oil, and geothermal energy. In this paper, inter-well porous ﬂow inside a rectangular reservoir is solved based on the complex variable function theory combined with the method of mirror images. In order to derive the solution analytically, the inter-well ﬂow is modeled as a 2D ﬂow in a homogenous and isotropic porous medium. The resulted exact analytical solution takes the form of an inﬁnite series, but it can be truncated to give high accuracy approximation. In terms of nine cases of inter-well porous ﬂow associated with enhanced geothermal systems, the applications of the obtained analytical solution are demonstrated, and the convergence properties of the truncated series are investigated. It is shown that the convergent rate of the truncated series increases with the symmetric level of well distribution inside the reservoir, and the adoption of Euler transform signiﬁcantly accelerates the convergence of alternating series cases associated with asymmetric well distribution. In principle, the analytical solution proposed in this paper can be applied to other scientiﬁc and engineering ﬁelds, as long as the involved problem is governed by 2D Laplace equation in a rectangular domain and subject to similar source/sink and boundary conditions, i.e., isolated point sources/sinks and uniform Dirichlet or homogeneous Neumann boundary conditions.


Introduction
The Laplace equation is a very important equation commonly occurring in many branches of physics as well as in other fields of science and engineering, notably in studies of electromagnetic phenomena, hydrodynamics, heat flow, and gravitation [1,2]. In this paper, we focus on analytically solving inter-well flows through porous media in a rectangular reservoir. Such problems are governed by the Laplace equation and frequently encountered in the production of groundwater, oil, and geothermal energy [3,4].
Although geological formations are inevitably involved in 3D flows through heterogeneous [5] and anisotropic [6] reservoirs, we consider the inter-well flow in this paper as a 2D homogenous and isotropic porous flow, with the permeability of the reservoir being treated as a constant scalar, instead of a more general form of tensor with spatial variation, in order to derive an analytical solution. The significance of such an analytical solution lies in that it can clearly reveal the fundamental functional dependence of inter-well flow on its controlling factors (e.g., driving pressure, well distribution, as well as reservoir geometry and properties) and provide first order approximation to many practical engineering problems [4].
Furthermore, an analytical solution has its unique advantages over a numerical solution. First, an analytical solution does not require discretizing the computational domain, and the solution is thus free of discretization error. Second, the computational cost for evaluating an analytical solution is much lower than obtaining a numerical solution. For inter-well flows, as the size of wells are normally much smaller than typical inter-well distances, very fine meshes are required for numerically capturing the sharp variations of physical quantities near the well regions, which results in huge computations in numerical modeling. Third, an analytical solution provides an ideal tool for quantifying the accuracy of numerical modeling. It is a common practice in numerical modeling studies to verify the accuracy of the employed numerical methods and computational meshes by numerically solving some benchmark problems and comparing the obtained numerical solutions directly with available analytical solutions. With the analytical solution proposed in this paper, inter-well porous flow problems can be established as benchmark problems for verifying numerical simulations of porous flows, including 3D heterogeneous and anisotropic porous flows in industrial applications.
Numerical modeling provides a powerful approach for simulating complicated engineering problems, but the accuracy of numerical modeling results is difficult to quantify, especially for 3D modeling results (partly due to the fact that 3D simulation normally involves huge computations that almost exhausted the computer resources and thus cannot afford full examination for the modeling accuracy by carrying on standard mesh and time step independence tests). The Richardson extrapolation method [7] is frequently invoked to assist accuracy assessment, but the required condition of systematically refining meshes often cannot be satisfied and may lead to underestimating the modeling errors. In such circumstances, comparing with a relevant 2D benchmark analytical solution will be a useful and practical test for estimating the 3D modeling error. For example, the error associated with a 3D modeling of a depth-dependent, heterogeneous, and anisotropic problem can be indirectly estimated from the error associated with a 3D modeling of a depth-independent, homogenous and isotropic but otherwise similar problem solved with comparable mesh and time step, while the error associated with the 3D modeling of the depth-independent, homogenous and isotropic problem can be well quantified through a comparison with an available 2D analytical solution.
In this paper, the proposed analytical solution for inter-well porous flow is derived based on the complex variable function theory [1,2], combined with the method of mirror images [8][9][10][11][12][13]. It should be noted that using mirror wells to solve symmetrical flow problem in porous media is not new and has been successfully applied to analytically solve inter-well flows inside unbounded domains with one or two straight boundaries, as well as those inside or outside a circular domain [4]. However, to the best of our knowledge, an analytical solution for inter-well porous flow inside a rectangular domain has never been published.
In the following, the derivation of the analytical solution is detailed in Section 2. Solutions for the pore pressure, Darcy velocity, potential function, and stream function are all proposed. The application of the analytical solution to inter-well porous flow problems associated with enhanced geothermal systems [14][15][16] are presented in Section 3, together with a discussion on the convergence properties of the analytical solution.

Analytical Solution
The inter-well fluid flow in a horizontal layer of porous media saturated with incompressible fluid is governed by the conservation laws of mass and momentum. The Darcy velocity satisfies the continuity equation and the Darcy's law [17] where p, k, and µ are pore pressure, permeability and fluid viscosity, respectively.
In a rectangular domain with impermeable boundaries at x = ±a and y = ±b, the normal gradient of pore pressure on the boundaries vanishes following Darcy's law, resulting in the following no-penetration boundary conditions: As the injection and production wells have a size much smaller than the inter-well distance, an injection well may be well represented by a point source, while a production well may be represented by a point sink. Therefore, the steady fluxes associated with the wells may be simulated by a number of isolated point sources or sinks with constant volumetric flow rates constrained by an overall mass conservation for steady porous flow The integration in Equation (5) is evaluated along a closed path surrounding a well, with ds being the arc differential along the integration path. A schematic of the above defined problem is shown in Figure 1a, in terms of an example case with one source and one sink.
Defining a scalar potential function ϕ that relates to Darcy velocity by i.e., Comparing Equations (2) and (7) yields Substituting Equation (8) into Equations (1)-(4) reveals that the potential function satisfies 2D Laplace equation with a homogeneous Neumann boundary condition, i.e., Furthermore, a scalar stream function ψ can be defined to relate to Darcy velocity by Such defined stream function provides a convenient tool for visualizing the velocity field, since the contours of the stream function form the streamlines. Substituting Equation (11) into Equations (1)-(4) reveals that the stream function satisfies the 2D Laplace equation and takes a uniform value along the boundary, i.e., subject to an assigned uniform Dirichlet boundary condition, As the combination of Equations (8) and (11) forms the Cauchy-Riemann Equations [2], the potential function and the stream function can be defined as the real and imaginary parts of an analytical complex function [2], i.e., Since the complex logarithm function is an analytical function [2], the superposition of multiple complex logarithm functions results in an analytical function. Accordingly, we obtain w(z) based on the method of mirror images [8][9][10][11][12][13], i.e., where is specified by the corresponding reference values of the potential and stream functions at a given point. In the series of Equation (15), the terms with m = n = 0 correspond to the real sources/sinks in the originally defined domain (Figure 1a), while the other terms correspond to the mirror images of the sources/sinks in the expanded domain, as shown in Figure 1b for the example case. In order to exactly satisfy the boundary conditions given by Equations (3) and (4), mirror images in both the x and y directions are required to be distributed symmetrically about all the four boundaries, which results in the infinite series solution given by Equation (15). However, for a practical application, the series solution does not need to be exact, and an appropriate truncation with finite M and N values is adequate to provide sufficient accuracy, i.e., Accordingly, the analytical solutions for the potential function, stream function, pore pressure, and Darcy velocity components can be expressed, respectively, by

Applications
The above obtained analytical solution can be readily applied to solve inter-well porous flow problems. Here, we demonstrate it for nine cases of enhanced geothermal systems, as shown in Figure 2. For all the cases, circles with positive sign represent injection wells, while those with negative sign represent production wells. Case 9 corresponds to a well array operating enhanced geothermal system [19], with four quarter production wells being located at the corners of the computational domain.
Among the nine well layout cases shown in Figure 2, Cases 3, 4, 5 and 8 have been numerically investigated by Chen and Jiang [18], while Case 9 is corresponding to a well array operating enhanced geothermal system [19], and the other cases are used for a comparison purpose. For all the cases, the geothermal reservoir is considered as a cube of 500 m by 500 m by 500 m with a permeability value of 10 −14 m 2 , and all the wells are vertical and have a radius of 0.15 m with a total water injection rate of 50 kg/s. These reservoir and flow conditions are typical for enhanced geothermal systems.
We adopt FORTRAN real*16 type computation to evaluate the analytical solution so that sufficient accuracy can be attained. Figure 3 shows the predicted pore pressure contours and streamlines for the nine cases. It is clear for each case that the four boundaries of the reservoir form a single streamline, while the pressure contours are all perpendicular to the boundaries, indicating that the uniform Dirichlet boundary conditions for the stream function, as well as the homogeneous Neumann boundary conditions for the pore pressure and the potential function, are well satisfied. Figure 3 indicates for all the cases that the near well regions are characterized by a higher velocity and pressure gradient. The pressure contours are essentially concentric circles near both the injection and production well walls, reflecting that the well walls are essentially boundaries of uniform pressure and velocity.  Figure 4 compares for the nine cases the predicted inter-well pore pressure variations along a straight segment connecting an injection well and a production well. The pressure difference shown in the ordinate axis of Figure 4 stands for the difference of the pore pressure at a given point along the inter-well segment subtracting the pore pressure at the wall of the corresponding production well, while the abscissa axis of the figure corresponds to the value of a nondimensional variable X, defined by the distance of a point on the inter-well segment to the injection well center divided by the corresponding distance between the injection well center and the production well center. The analytical results shown in Figure 4 can be used to quantify modeling errors associated with numerical simulation results. For example, the relative difference for Case 9 between the analytical solution shown in this figure and a finite element solution based on 1.2 million degrees of freedom [19] is less than 10 −9 . The convergence properties of the inter-well pressure difference predicted by the series solution along with increasing M and N values are quantatively investigated for the nine well layout cases. It turns out that the truncated series solutions for Cases 4, 7, 8, and 9 converge rapidly and monotonically, while those for Cases 1, 2, 3, 5, and 6 are associated with alternating series and converge slowly. This result suggests that the convergence behavior of the analytical solution is closely related to the symmetric nature of the well layouts. As is clear from Figure 2, the wells for the four cases with rapid and monotonic convergence are all distributed symmetrically about both the x and y axes, while the wells for the other five cases are all distributed asymmetrically about the x axis.
By invoking Euler transform [20,21] N ∑ n=1 (−1) n−1 a n = to the alternating series solutions for Cases 1, 2, 3, 5, and 6, the convergence processes for all the five cases are accelerated dramatically, as evident in Figure 5, indicating that the adoption of Euler transform significantly accelerates the convergence of alternating series cases associated with asymmetric well distribution. However, the five cases exhibit two different trends: the convergent rates for Cases 2, 3, and 6 constantly increase with M and N values, as shown in Figure 5a  The convergence processes of inter-well pressure difference for the nine cases are compared in Figure 6. Although Cases 4, 7, 8, and 9 all converge with nearly constant rates, the convergent rates of Cases 7, 8, and 9 are much higher than Case 4. This difference can be well explained by the additional symmetric property of the well layouts associated with Cases 7, 8, and 9: the wells for these three cases are also symmetric about the x = y and x = −y diagonals. Moreover, it is interesting to observe that the convergent curves of Cases 1 and 5 for M = N > 20 approach to a nearly constant slope that is comparable to that of Case 4. These results are clear evidence for the dependence of the convergent rate of the analytical solution upon the symmetric level of the well distribution inside the reservoir. Figure 6. Convergence processes of inter-well pressure difference for the nine cases, with Euler transform being invoked for Cases 1, 2, 3, 5, and 6.

Concluding Remarks
An analytical solution for inter-well porous flow in a rectangular reservoir is proposed in this paper. The applications of the solution to nine cases of enhanced geothermal systems are demonstrated as examples, and the convergence properties of the analytical solution are discussed accordingly. It should be noted that the nine ideal cases investigated in Section 3 are chosen just as a demonstration. In fact, the analytical solution proposed in this paper has sufficient flexibility to deal with different aspect ratios of the rectangular reservoir and different distributions of isolated sources/sinks, as long as the condition of overall mass conservation for steady porous flow, Equation (6), is satisfied. In principle, the analytical solution proposed in this paper can be readily applied to other scientific and engineering fields, as long as the involved problem is governed by 2D Laplace equation in a rectangular domain, with uniform Dirichlet or homogeneous Neumann boundary conditions, and subject to isolated point sources/sinks. This work may be furthered in the following directions. First, the method of mirror images may be employed to derive analytical solutions for the 2D Laplace equation in other polygon regions. In principle, the procedure for deriving such solutions is essentially the same as that for a rectangular domain, but the mirror image distributions will follow more complicated patterns for polygons with non-orthogonal boundaries. Second, the methodology adopted in this work may be extended straightforward to derive analytical solutions for problems governed by 3D Laplace equation. Finally, the methodology of superposition of analytical functions may be applied to solve problems with regularly arranged sources/sinks in unbounded domains. For example, Case 9 of this paper represents an infinite well array with staggered injection and production wells arranged in orthogonal straight lines. If we replace the locations of the mirror images in Equation (15) by those of true sources/sinks distributed in unbound domains, the requirement of orthogonal distribution of the sources/sinks can be eliminated as a consequence of that the symmetric boundary conditions no longer need to be satisfied. Then, modified forms of Equation (15) might be applied to solve similar problems, such as the lattice sums problems [22,23] and the Helmholtz equation problems with vanishing index of refraction [24,25].