A Singularity Removal Method for Coupled 1D-3D Flow Models

In reservoir simulations, the radius of a well is inevitably going to be small compared to the horizontal length scale of the reservoir. For this reason, wells are typically modelled as lower-dimensional source terms. In this work, we consider a coupled 1D-3D flow model, in which the well is modelled as a line source in the reservoir domain and endowed with its own 1D flow equation. The flow between well and reservoir can then be modelled in a fully coupled manner by applying a linear filtration law. The line source induces a logarithmic type singularity in the reservoir pressure that is difficult to resolve numerically. We present here a splitting method for the solution, where the reservoir pressure is split into a low regularity term that explicitly captures the singularity, and a well-behaved, regular component $v$. The singularities can then be removed from the system by subtracting them from the governing equations. This results in a reformulated coupled 1D-3D model in which all variables are smooth. It can therefore be approximated using standard numerical methods.


Introduction
Accurate well models are of critical importance for reservoir simulations, as the well often constitutes the driving force for reservoir flow, as well as being the main access point of information about its state. The main challenge of well modelling is rooted in scale disparity; the well has a radius of ∼ 10 cm, while the reservoir might extend several kilometres in the horizontal plane. For this reason, wells are typically modelled by the use of lower-dimensional source terms.
In this work, we take as a starting point the coupled 1D-3D flow model q + κ µ ∇p = 0 in Ω, q +κ µ dp ds = 0 in Λ, where Ω ⊂ R 3 denotes the reservoir domain and Λ = ∪ wells w=1 Λ w ⊂ R 1 a collection of line segments Λ w each representing a well. In this work, we assume the wells to be spatially separated. The wells are considered as line sources in the reservoir equation (1a)-(1b), and endowed with their own 1D flow equations (1c)-(1d). Here, s denotes the arc length s of Λ, p and q the fluid pressure and flux in the reservoir, κ, µ positive, constant values giving the reservoir permeability and fluid viscosity therein, andp,q,κ andμ the corresponding variables in the well.
The mass flux q between well and reservoir is modelled using a linear filtration law, which states that the flow between them is proportional to their pressure dif- ference: Here, the coefficient β ∈ C 2 (Λ) is allowed to vary along the well. The pressure difference is not taken directly, but using the reservoir pressure averaged around the borehole,p. Letting R denote the radius of the well, and considering the cylindrical domain and line segment Λ illustrated in Figure 1, this averaged pressure is given bȳ p(z, R) = 1 2πR 2π 0 p(R, z, θ)dθ. ( Finally, the flow from well to reservoir is modelled by means of a line source δ Λ on Λ, which we understand in the following manner: for all φ ∈ C 0 (Ω), with s w denoting the arc-length of line segment Λ w . Elliptic equations with line sources of the type (4) have been used in a variety of applications, e.g., the modelling of 1D steel components in concrete structures [26] or the interference of metallic pipelines and borecasings in electromagnetic modelling of reservoirs [35]. A coupled 1D-3D heat transfer problem was considered in the context of geothermal energy in [5], where it was used to model heat exchange between (3D) soil and a (1D) pipe. Coupled 1D-3D flow models have also been studied in the context of biological applications, such as the efficiency of cancer treatment by hyperthermia [28], the efficiency of drug delivery through microcirculation [9,32], and the study of blood flow in the vascularized tissue of the brain [18,33]. In this work, we restrict ourselves to considering its application in the context of reservoir modelling.
The main challenge with the coupled 1D-3D flow problem is that the line source induces the reservoir pressure to be singular, thereby making its analysis and approximation non-standard. Typically, reservoir simulations are performed using finite volume methods. The discretized form of the coupling in (1a)-(1d) is then given by where p K denotes the average pressure in the grid block containing the well. Due to the singularity, p K may not be representative of the reservoir pressure at the borehole; this is typically accounted for by altering β with a well index J, i.e., β → β * with β * = J β. A correction of this type was first developed by Peaceman in [31], where he considered the two-point flux approximation method on uniform, square grids when the well is aligned with one of its axes. Via an analytic solution valid for simplified cases, he gave a well index depending on the equivalent radius of the well, i.e., the radius at which the reservoir pressure equals the well block pressure. The equivalent radius depends, among other factors, on the discretization scheme, placement of the well relative to the mesh, and reservoir permeability. The problem of finding appropriate well indexes has been treated in a multitude of works; Peaceman himself treated an extension of his method to non-square gridblocks and anisotropic permeability [30]. The extension to more generalized grids was treated by e.g. Aaavatsmark in [1][2][3], to more generalized flow models by e.g. Ewing in [16], and to more generalized discretization schemes by e.g. Chen et al. in [10]. Many authors have contributed to the extension to generalized well placements, we mention here the work of King et. al in [21], Aavatsmark in [4], and of special relevance to our work, that of Wolfsteiner et al. in [36] and Babu et al. in [7]. We take here a different approach, in which the singularities are explicitly removed from the governing equations. To be more precise, we show here that the reservoir pressure p admits a split into a collection of low-regularity terms where the singularity is explicitly captured, and a well-behaved, regular component v ∈ H 2 (Ω): Here, r w is the distance from a point to the linear extension of Λ w , and ln(r w ) describes the singular nature of the reservoir pressure around the well. E w denotes an extension operator E w : Λ w → Ω, the precise meaning of which will be made clear at a later point. With the splitting in hand, we can then remove the singular terms from the system; in practice, this leads to the reformulated coupled 1D-3D flow model (33a)-(33d). The reformulated model is posed in terms of the smooth variablesp and v, and can thus be approximated using standard numerical methods.
The technique of removing singularities is commonly known for point sources; we refer here to [15, p. 14] for a more in depth explanation. It has previously been studied in the context of reservoir models by e.g. Hales, who used it to improve well modelling for 2D reservoir models [19]. A similar splitting result as the one we utilize was introduced by Ding in [14] for the point source problem, where it was to formulate grid refinement strategies to improve the approximation of the flux. We are, however, to the best of our knowledge the first to formulate a singularity removal method for the coupled 1D-3D flow problem.
The singularity removal, and subsequent reformulation of the model in terms of the smooth variables v andp, is similar to the Peaceman well correction in that it leads to an alteration of the inflow parameter β. It differs, however, in that it works on the continuous level. It is therefore easily adapted to different discretization methods, as well as generalized well placements. While we restrict ourselves here to considering wells that are spatially separated, our method can be extended to handle also wells that branch. To do so, one needs only to exchange the ln(r w ) term with a function G w that satisfies −∆G w = δ Λw . We refer here to our earlier work [17, p. 8] for a candidate G w . It was constructed by integrating the 3D Green's function of the problem over the line source, the same method as was employed by Wolfsteiner et al. and Babu et al. in [7,36] to construct analytical solutions with which to calculate the well index. The use of Green's functions to construct analytical well models has a rich history, e.g., Nordbotten et al. used Green's functions to construct analytical models to estimate leakage of CO 2 stored in geological formations in [29]. Our method differs, however, from these earlier approaches in that we use only G w to handle the line source. We do not use the Green's function to account for boundary conditions; they are handled instead by v. Our method is therefore applicable to any sufficiently smooth domain. Lastly, let us point out that since our method gives an explicit representation of the logarithmic nature of the solution, it does more than alter the well index. With the pressure splitting (6), the reservoir pressure can be accurately represented in the whole domain, i.e., including in the near-vicinity of the well.
In our presentation of the method, we limit ourselves to considering Poiseuille flow in the well and to having a linear reservoir equation. Only the former of these is critical to the methodology; the well equation could, for example, contain non-linearities. For the sake of brevity, we consider only wells that are spatially separated. Lastly, we assume constant, scalar-valued permeabilities. For an extension to tensor-valued permeabilities and non-linear reservoir equations, we suggest using the solution splitting in (6) to formulate a generalized finite element method [34], where the analytic functions capturing the solution singularity are used to enrich the finite element space itself.
For the discretization and numerical experiments, we consider herein the use of Galerkin Finite Element (FE). The FE approximation of the line source problem was studied by D'Angelo in [11] by means of weighted Sobolev spaces, using similar techniques as those known for e.g. corner-point problems [8]. He proved in this work that the approximation converges sub-optimally unless the mesh is refined around the well. These suboptimal convergence rates were found to be local to the line source by Köppl et al. in [25], meaning that they do not pollute the approximation outside the well block. The approximation of the coupled 1D-3D problem will, however, still suffer if R is smaller than the mesh parameter h, as the averaged reservoir pressurep is then defined inside the well block. In practice, one therefore needs a fine mesh around the well for the standard FE approximation of (1a)-(1d) to converge. This makes the problem computationally expensive to solve. With this in mind, Kuchta et al. studied suitable preconditioners in [24], that handle the ill-conditioning that the coupling introduces in the discrete matrix A. Holter et al. then applied this preconditioner to simulate flow through the microcirculature found in a mouse brain [20]. An alternative coupling scheme was introduced by Köppl et al. in [22], where the source term was taken to live on the boundary of the inclusions. The result is a 1D-(2D)-3D method where the approximation properties have been improved, at the expense of having to resolve the 2D boundary of the well.
The article is structured as follows. We start in Section 2 by defining the relevant function spaces. In Section 3, we then introduce in more detail the coupled 1D-3D flow model we take as a starting point. In Section 4, we introduce a splitting theorem for the reservoir pressure, showing that it admits a splitting into a low-regular term that explicitly captures the singularity, and a smooth function v. The splitting theorem can then be used to subtract the singularities from the governing equations. The result is the reformulated coupled 1D-3D flow model (33a)-(33d), posed in terms of the smooth variablesp and v. As the singularities have been removed from the system, this formulation enjoys significantly enhanced regularity, and can be approximated using standard numerical methods. The variational formulation and FE discretization of the reformulated problem are given in Sections 5 and 6, respectively, and require only standard function spaces. In Section 7, we discuss how this discretization of the reformulated model resembles a Peaceman well correction. We then conclude the article with a numerical experiment, where we test the Galerkin FE method of both the standard and reformulated coupled 1D-3D flow model. We show that the singularity removal recovers optimal convergence rates on uniform meshes, i.e., without needing to refine the mesh around the well. Moreover, in a manner similar to altering the well index, it makes the approximation robust with respect to the ratio R/h.

Background and notation
The purpose of this section is to introduce the appropriate function spaces for the coupled 1D-3D flow model. Let H k (Ω) be the Sobolev space, with β denoting a multi-index and D β the corresponding weak distributional derivative of u. H k (Ω) is a Hilbert space endowed with inner product We use a subscript to denote the subspace of H k with zero trace on the boundary, H k 0 , i.e., As we will see, the reservoir solution p in (1a)-(1d) fails to belong to H 1 (Ω) due to singular behaviour on Λ. For this reason, we consider also a weighted Sobolev space. To define it, let −1 < α < 1, and take L 2 α (Ω) to denote the weighted Hilbert space consisting of measurable functions u such that where r denotes the distance of a point to Λ, i.e., r(x) = dist(x, Λ). This space is equipped with the inner product For α > 0, the weight r α has the power to dampen out singular behaviour in the function being normed; for α < 0, the weight function can induce or worsen already singular behaviour. We therefore have the re- A practical use of this space is found, for example, considering the logarithmic grading (refinement) that is often performed on a mesh around the well. The well introduces a logarithmic type singularity in the reservoir pressure that cannot be resolved using linear elements. Consequentially, the convergence rate of standard numerical methods degrade using uniform meshes. Optimal convergence can be retrieved by a specific refinement of the mesh around the well [6,11,14]. The exact convergence rates and mesh grading requirements are closely related to the weighted Sobolev space wherein the solution exists; in fact, the graded mesh will be uniform with respect to the weight function r α .

Mathematical model
Here, we introduce in more detail the coupled 1D-3D equation we take as a starting point. Let Ω ⊂ R 3 denote a bounded domain describing a reservoir, with smooth boundary ∂Ω. We consider here steady-state, incompressible Darcy flow where q and p denote reservoir flow and pressure, µ the fluid viscosity, and κ a given, positive, scalar permeability. We consider also a collection of line segments representing wells, Λ = ∪ wells w=1 Λ w , where each well is considered to be a tube of fixed radius R with centreline Λ w . The centreline is parametrized by the arc length s w , and its normalized tangent vector denoted by τ τ τ sw . As the radius of the tube is small, we assume the radial and angular components of the well solution can be neglected, meaningp| Λw =p(s w ). Poiseuille's law for conservation of momentum and mass then readŝ (8b) withq w andp w denoting flow and pressure in the well and q the linear mass flux into or out of the well. d dsw denotes the derivative with respect to the tangent line, or equivalently, the projection of ∇ along τ τ τ , i.e., d dsw = ∇·τ τ τ sw . As the fluid flux in the well has a fixed direction, it can be given as a scalar functionq w , characterized by the propertyq w =q w τ τ τ w .
Letting now Λ = ∪ wells w=1 Λ w denote the collection of line segments Λ w , the well pressure and flux can be written as 1D variablesp,q : Λ → R. The well and reservoir flow can then be coupled together using a linear filtration law, which states that the mas flux q between the two, given as the rate of transfer per unit length, is proportional to their pressure difference: The variable λ ∈ C 2 0 (Λ) is used to account for the fact that the well may not be in perfect contact with the reservoir, giving rise to a pressure drop ∆p skin across the skin of the well. The mass flux q between well and reservoir is then given by q = 2πRλ∆p skin . To ease the readability of the next section, we assume λ = 0 at the endpoints of the line segment.
The pressure difference f (p,p) between well and reservoir is not taken directly, but uses an averaged valuē p(z; R) for the reservoir pressure given in (3). This can be interpreted physically as the reservoir pressure averaged around the borehole. The flow in well and reservoir can be then modelled, in a fully coupled manner, by the set of equations whereκ = R 2 /8, β = 2πRλ and p D ∈ C 2 (Ω) and p D (Λ) some given boundary data. The flow from well to reservoir is here modelled by means of a generalized Dirac delta function δ Λ , which we understand in the sense of (4). Finally, this system can be reduced to its conformal form by eliminating the 1D and 3D fluxes:

Splitting Properties of the Solution
In this section, we show that the line source in the righthand side of (11a) introduces a particular structure to the solution of the coupled 1D-3D flow problem. We do this by means of a splitting technique, in which the reservoir pressure is split into a low regularity term that explicitly captures the singularity, and a regular component v being the solution of a suitable elliptic equation.
To start with, we treat the case where Λ is a single line source aligned with the z-axis, κ µ = 1 and the well outflow q is a given function f ∈ C 1 0 (Λ). For this simplified scenario, the splitting is straightforward; it therefore serves well in illustrating the method itself. We then generalize it in two steps; handling first an arbitrary line segment and κ µ = 1, and finally the coupling between reservoir and well. Finally, we use the splitting to reformulate the coupled 1D-3D flow problem into the system (33a)-(33d), wherein the singularity has been removed and all variables are smooth.

Elliptic equations with a single line source
In this section, we consider the elliptic equation when Λ and Ω are as illustrated in Figure 1, and f = f (z) ∈ C 1 0 (Λ) is a given, smooth line source intensity, assumed to be zero at the endpoints of Λ. The solution p then admits a splitting into an explicit, low-regularity term f (z)G(r), and an implicit, high-regularity term v: The singular term G(r) is further given by G(r) = Ψ (r) ln(r), with Ψ (r) denoting some smooth cut-off function satisfying where 0 < R < R c . Assuming the cut-off radius R c is chosen small enough to satisfy Ψ (r) = 0 on ∂Ω, the regular component v is then defined as the solution of where F = f (z)G(r) + f (z) (Ψ (r) ln(r) + 2 ln(r) Ψ (r)) .
To see that this p indeed solves (12), let us first note that − 1 /2π ln(r) is in fact the fundamental solution of the Laplace equation in 2D. It thus has the property Considering then the Laplacian of p, all terms vanish by construction except for f (z)Ψ (r)∆ ln(r)/(2π). Integrating this term over the domain, we find that It follows that the p constructed in (13) indeed solves (12) in a suitably weak sense. Formally speaking, the splitting works by adding a logarithmic term for which the Laplacian returns the line source with the required intensity f , and using the regular component v to correct the solution so it solves the original problem. The existence of such a v follows from standard elliptic theory. Letting S z0 be the intersection of Ω with the plane perpendicular to Λ at a given point z 0 , i.e., S z0 = {(r, θ, z) ∈ Ω : z = z 0 }, it can be shown by calculation that G(r) ∈ L 2 (S z0 ). The same holds for the terms Ψ (r) ln(r) and Ψ (r) ln(r) . Note here that ln(r) does not belong to L 2 (S z0 ) due to the singular behaviour at r = 0; however, we chose Ψ (r) so that Ψ (r) = 0 in a small circle around the line source. Next, f (z) and f (z) both belong to L 2 (Λ). It thus follows that the entire right-hand side in (15a) must belong to L 2 (Ω), and that a v ∈ H 2 (Ω) exists solving (15a)-(15b).
Moreover, a straightforward calculation shows that ln(r) fails to belong to H 1 (Ω) due to the singular behaviour on Λ. Instead, we have p ∈ H 1 α (Ω) for any α > 0. The splitting (13) therefore forms a split into high and low-regularity terms, meaning that v is smoother and better behaved than the full solution p. This observation will be central to the numerical method considered in Section 6.

Elliptic equations with an arbitrary line source
In this section, we consider the elliptic problem when the right-hand side is a line source δ Λ , located on a single line segment Λ with endpoints a, b ∈ Ω. Letting again f = f (s) ∈ C 1 0 (Λ) be a given, smooth line source intensity, the solution p then admits a splitting into an explicit, low-regularity term E(f )G(r), and a high-regularity component v: Here, G(r) = Ψ (r) ln(r), as before, and E denotes an extension operator E : H 2 (Λ) → H 2 (Ω) extending f so that it can be evaluated in the entire domain Ω. In this work, we choose it to be given by where P is the orthogonal projection operator P (x) = τ τ τ ·(x−a) that maps a point x ∈ Ω to the closest point on the extension of Λ. Here, τ τ τ denotes the tangent vector of Λ, τ τ τ = (b − a)/ b − a , see Figure 2 for an illustration. The extension thus acts by mapping a point in the domain onto the arc length parameter s corresponding to the closest point on the line, and evaluating f for this s. The extension E(f ) depends then only on s. Notice that for the line segment aligned with the z-axis, as was considered in Section 4.1, the extension would Assuming again R c is chosen small enough to satisfy Ψ (r) = 0 on ∂Ω, the regular component v is then defined as the solution of where F = d 2 ds 2 E(f )G(r) + E(f ) (Ψ (r) ln(r) + 2 ln(r) Ψ (r)) .
To see that the constructed p indeed solves (18), let us consider its Laplacian. By construction, all terms disappear except E(f )Ψ ∆ ln(r)/(2π). Integrating this term over the domain, we find that for all ∀φ ∈ C 0 (Ω), where we used the property that E(f ) = f on Λ. It follows that the p constructed in (13) indeed solves (12) in a suitably weak sense. Now, ∆E(f ) and E(f ) both belong to L 2 (Λ). By the discussion made in the preceding section, it then follows that the entire right-hand side of (21a) must belong to L 2 (Ω). A v ∈ H 2 (Ω) then exists solving (21a)-(21b). As was noted before, ln(r) ∈ H 1 α (Ω) for any α > 0, but fails to belong to H 1 (Ω) due to the singular behaviour on Λ. The decomposition (19) thus forms a split into high and low-regularity terms, meaning that v is smoother and better behaved than the full solution p. For a more in-depth explanation of the splitting method, as well as a treatment of discontinuities in f , we refer to our previous work in [17].

The coupled 1D-3D flow problem
Let us now consider the coupled 1D-3D flow problem (11a)-(11d). From the discussion in the preceding section, it is natural to assume p solving this problem admits a solution splitting of the type: with f being the previously introduced pressure difference f =p − 2πp, and the regular component v defined as the solution of where F (p,v; β) = d 2 ds 2 E(βf )G(r) + E(βf ) (Ψ (r) ln(r) + 2 ln(r) Ψ (r)) .
Unlike in Sections 4.1 and 4.2, f = f (p,p) is now implicitly given fromp andp solving the coupled 1D-3D flow problem. To reformulate (11a)-(11d) in terms ofp and v, the right-hand side needs to be reformulated; to this end, let us first treat the pressure differencep−2πp.

By the splitting (24) and the definition of the averaging in (3), calculations reveal that
meaning we can write From this we state the reformulated coupled 1D-3D flow model: and β * is given by The extension to multiple wells follows naturally by applying the superposition principle. This yields the following splitting for the reservoir solution p solving (11a)-(11d): where E w : H 2 (Λ w ) → H 2 (Ω) is the same extension operator as before, G is a logarithmic type term r w denotes the distance from a point to the linear extension of Λ w , i.e., r w = dist(x,Λ w ), Ψ the same cut-off function that ensures G(r w ) = 0 on ∂Ω for all w, and v solves The system (33a)-(33c) constitutes a reformulation of the coupled 1D-3D flow model in terms of the smooth variables v andp. For an example of what the splitting might look like, the reader is invited to examine Figure  3. As the singularities have here been removed from the system, it enjoys significantly improved regularity compared to the standard formulation (11a)-(11b). Consequently, it can be approximated using standard numerical methods.

Weak formulation
In this section, we state a weak formulation of the reformulated coupled 1D-3D flow problem (33a)-(33d). As the variables in this formulation are all smooth functions, this can be done using standard Sobolev spaces. For the sake of completeness, we give also a weak formulation of the standard coupled 1D-3D flow problem (11a)-(11d). The reservoir pressure p therein contains a singularity; for this reason, its weak formulation requires the use of weighted Sobolev spaces. Consider first the reformulated coupled 1D-3D flow problem. Let V denote the product space V = V ×V , where Multiplying (33a) and (33c) with test functions φ ∈ H 1 0 (Ω) andφ ∈ H 1 0 (Λ), respectively, integrating over their respective domains, and performing an integration by parts, we arrive at the following variational formulation: where The full reservoir pressure can then be constructed from v by the relation Next, let us consider the standard coupled 1D-3D flow model, and give its variational formulation as it was proposed in [12]. Let V α denote the weighted prod- Multiplying (11a) and (11c) with test functions v ∈ H 1 −α,0 (Ω) andv ∈ H 1 0 (Λ), respectively, integrating over their domain of support, and performing an integration by parts, we arrive at the variational formulation: and the test space V −α,0 is the space of functions (φ,φ) ∈ V −α,0 with zero trace on the boundary. Notice here that the test and trial spaces are chosen with opposite weight functions; this is what ensures the continuity and coercivity of the bilinear form (46). For a proof of the well-posedness of this formulation, the reader is referred to [11,12].

Numerical Discretization
In this section, we show the block matrix resulting from a finite element discretization of weak formulation of the reformulated coupled 1D-3D problem. As the pressure difference f (p,v) =p −v now uses the regular part of the pressure, v ∈ H 2 (Ω), we introduce here also the simplificationv h = v h | Λ ; i.e., we take the trace of v h on Λ rather than the average over the cylinder. This is motivated by the fact that R is assumed negligible compared to the mesh size h, and v is regular, meaningv ≈ v| Λ . The result is a "true" coupled 1D-3D flow model, in that it considers only 1D and 3D variables, with no averaging performed over a 2D cylinder. The same approximation is not possible for the standard coupled 1D-3D flow model as the reservoir pressure is there undefined on Λ.
We will now give the discretized form of the variational formulation (38). For simplicity, let us assume Ω is a polyhedron that readily admits a partitioning T T,h into simplicial elements T : The simplicial partitioning T T,h forms a mesh, assumed conforming, which can then be characterized by the mesh size h = max T ∈T T ,h h T . Next, we associate this mesh with the usual (3D) Lagrange space of order 1, V h u , given by Here, P 1 denotes the space of polynomials of degree 1, and C 0 u (Ω) the space of continuous elements that equal the interpolation of u on the boundary, i.e., Next, we assume Λ admits a partitioning T I,h into line segments I: assumed again to satisfy all the requirements of a conforming mesh, and associated with the mesh sizeĥ = max I∈T I,h h I . For the discretization ofV , we use the (1D) Lagrange space of order 1, with C 0 u (Λ) interpreted as in (47). Considering first the reformulated system (38 where {φ 1 , φ 2 , ..., φ N } and {φ 1 ,φ 2 , ...,φN } are linear hat functions spanning V h andV h , respectively. Testing (38) with v = φ i for i = 1, ..., N andv =φ j for j = 1, ...,N , we arrive at the following discrete system: where A i,k andÂ j,l are the standard stiffness matrices Recall now that v h is a linear function used to approximate the high regularity term v ∈ H 2 (Ω). Assuming R h, the averagev can be well approximated by simply taking the trace v| Λ . The pressure differencê p −v is then given bŷ Here, T : V h →V h is the discrete trace matrix, characterized by the property φ k | Λ = N l=1 T k,lφl . Lettinĝ M denote the 1D mass matrix,M j,l = (φ j ,φ l ) Λ , and C denote the coupling block, we find the following block matrix form of the system: We will refer to this system as the Singularity Removal Based FE method. After solving (56), a discretization of the full reservoir pressure p h can be reconstructed using where I k h denotes the interpolation onto the Lagrange space of order k. As the interpolation of ln(r) is fairly cheap, the approximation property of p h can here be improved by choosing the interpolation degree k high.
A more straightforward method can be found by discretizing (45) directly; this is the finite element formulation analysed in e.g. [11]. As we will compare the performance of this method against the Singularity Removal Based FE method, we give here its discretization for the sake of completeness. Setting testing (45) with v = φ i for i = 1, ..., N andv =φ j for j = 1, ...,N , we arrive at the following discrete system: The discretization of the pressure difference between well and reservoir involves an averaging of the 3D test function. In order to formulate the discrete averaging operator (3), we consider the discontinuous Galerkin space of order 0, A straightforward calculation then reveals where Π is the discrete averaging matrix Π : V h → X h and {ψ 1 ,ψ 2 , ...,ψM } are the basis functions spanninĝ X h .
Letting now N be the mass matrix given by N m,l = (ψ m ,φ l ), we arrive at the following block system for the discretization of (45): We will refer to this system as the standard FE method.

Relation to the Peaceman well model
In this section, we show that the reformulated coupled 1D-3D flow model (33a)-(33d) under certain conditions reduces to the Peaceman well correction. We start by giving a brief summary of the methodology Peaceman introduced in his seminal work [31]. We then return to our reformulated model, and show that with G(r) chosen so that its support is the equivalent radius of the Peaceman well correction, the reformulation results in a well index that equals the one derived by Peaceman.
In reservoir simulations, the mass flux between well and aquifer, q, is usually modelled in a manner analogous to that in (11a): where p w is the flowing pressure in the well, J its well index, and p K the reservoir pressure averaged over the grid cell K. In Section 4, we showed how the line source that models the well introduces a logarithmic type singularity in the reservoir pressure. For wells with radius much smaller than the grid size h, i.e., R h, p K is therefore likely to constitute a poor representation of the reservoir pressure in the near vicinity of the well. The Peaceman well model accounts for this by altering the well index J in (62) so that q better corresponds to the numerical approximation of the pressure difference between well and aquifer. Assuming radial flow, Darcy's law in a heterogeneous reservoir is given, per unit well length, by the relation Integrating this equation to a radius r e , 2πκ qµ pe pw dp = − re R dr, we find that when p e = p(r e ). We also need to take into account the pressure drop ∆p skin across the skin of the well. To do so, let S be the skin-factor, defined by the relation Letting now r e be the radius at which the reservoir pressure equals the averaged grid cell pressure p K , Peaceman used the following relation between q and the pressure difference p w − p K [31]: (a) FE approximationsp h and v h (b)p h and the reconstructed logarithmic (c) Well pressurep and reconstructed reservoir pressure p h , found using (57). To utilize this correction, one must first identify the equivalent radius r e entering in (67). This radius generally depends on the discretization method, the location of the well within the grid, and the permeability of the rock around the well. Assuming for example square grid blocks and a well at the center of an interior grid block, Peaceman derived an equivalent radius r e = 0.2h for the two-point flux approximation [31].
The reformulation of the pressure difference f in terms ofp and v bears a strong resemblance to the Peaceman well correction in (67). In a practical sense, the reformulation into (33a)-(33d) can be interpreted as a non-local well correction, which has a support in a region around the well which may significantly exceed the grid resolution. To see more clearly the similarity with the Peaceman well correction, let us now consider a single well; we have then Next, we let nowp be the flowing well pressure p w . The term G(r) contains the logarithmic component of the solution; in a manner analogous to the Peaceman well correction, we make it local to the cylinder of radius r e by setting G re (r) = − 1 2π ln(r/r e ) for r ≤ r e , 0 otherwise.
Note that this G is not smooth enough to work for the solution split (31), we use it here for the sake of comparison. Inserting it in (68) yields the relation Here, 2πκ/µβ can be substituted by the skin factor of the well by recalling q = β∆p skin = 2πRλ∆p skin . This results in an expression that equals the Peaceman well correction given in (67), i.e., The regular component v can be interpreted as a sort of background pressure, or more precisely, the component of the reservoir pressure that can be approximated using linear functions. We see then that the singularity removal constitutes an alteration of β (which can be interpreted as a well index) so that the mass flux function q better corresponds to the numerically computed pressure difference between well and reservoir, i.e.,p−v. For this reason, we expect that the singularity removal, in a manner similar to the Peaceman well correction, will improve the stability of the FE approximation with respect to the ratio R/h.

Numerical Results
In this section, we perform numerical experiments to test the approximation properties of the Singularity Subtraction Based FE method given by (56). For the implementation, we utilized the finite element framework FEniCS [27]. We compare the results against those obtained using the standard FE given by (61). Our implementation of this method uses an earlier implementation from Kuchta [23], the same as was utilized for the results of Holter et al. in [20]. The Singularity Removal Based FE method was implemented by an extension of this code, using also the mixed-dimensional  functionality of FEniCS developed and implemented by Daversin-Catty [13]. We test the capability of each method in approximating the test problem Of special interest is the stability of the approximation when the well radius is small compared to mesh size h; for this reason, we test using four different values for the well radius; R ∈ {1.0e-1, 1.0e-2, 1.0e-3, 1.0e-4}. The rest of the problem parameters are then as follows: whereβ is the inflow parameter assigned to the 1D flow equation. Figure 4 shows the approximation errors, measured in the L 2 -norm, when the problem was solved using a sequence of increasingly fine meshes. The blue lines in Figure 4a show the approximation error of v h , measured in the L 2 -norm, i.e., v h − v a L 2 (Ω) with v a being the analytic solution in (74b). For R < h, the errors are seen to be invariant with respect to R, and the approximation of v h exhibits moderate superconvergence. To expand upon this, we expect for this approximation optimal convergence rates of order h l with l = 2.0; we see here a slight super-convergence as l = 2.2. For h > 1/8 and R = 0.1, our assumption of R < h is no longer valid, and we see a degradation of the convergence rates.
To be more precise, we made in the construction of the block matrix (56) the simplificationv = v| Λ , and this is not justified for R ∼ h. Optimal convergence rates could be restored by taking the average of v h rather than its trace.
The red lines in Figure 4a give the approximation errors for the full reservoir pressure using the standard FE method described by (61). We give here the approximation error of p h in the L 2 -norm, i.e., p h −p a L 2 (Ω) . For the standard FE method, the convergence properties strongly depend on the well radius R, with decreasing R leading to a reduction in the convergence rate. The best convergence rates are seen when R ∼ h, but even here, the convergence is sub-optimal compared to the Singularity Removal Based FE method. This can be explained by noting that the standard FE method explicitly resolves the line source in the problem; it was shown in [11] that this leads to a reduction in the convergence rate of p h . We refer here to our comments in [17, p. 14-15] for a more in-depth explanation of this, and re-   mark only that the line source problem is expected to converge with order h 1− for > 0 arbitrarily small. Thus, the convergence order h l with l = 1.4 surpasses the theoretical expectation when R ∼ h.
The blue and red lines in Figure 4b give the approximation error ofp h using the Singularity Removal Based and standard FE method, respectively. The approximation error is also here measured in the L 2 -norm, i.e., using p h −p a L 2 (Λ) . We see here that the singularity removal significantly improves the convergence properties of the problem for R < h. The convergence rates degrade when R > h. This is again due to the simplificationv = v| Λ used in the construction of the block matrix (56), and could be resolved by removing this simplification.
From 4b, is clear that the standard FE method has trouble approximating the solution when R < h. Moreover, the approximation error ofp h is seemingly more sensitive than p h with respect to the ratio R/h. This can be understood by returning to the splitting and noting that the error in p h is due to three separate issues, namely, the error in the approximation of the pressure difference, i.e., p h − 2πp h − (p a − 2πp a ) L 2 (Λ) , the error in approximating the logarithm, i.e., ln(r) h − ln(r) L 2 (Ω) , and the error in approximating v (which is comparatively small). To understand the first of these further, Figure 5 shows the pressure along the line x, x = {(x, y, z) ∈ Ω : y = 0.5, z = 0.75}, obtained using the FE method with and without the singularity removed. We see that the standard FE method has trouble resolving the logarithmic nature of the reservoir pressure around the well, leading to a large approximation error inp. This further pollutes the approximations of bothp and p. The effect is not as noticeable for p as its approximation error is dominated by the approximation error for ln(r). The well pressurep, however, is in principle a smooth function, for which the FE approximation should be comparatively small. Its approximation error is therefore dominated by the term p a −p h L 2 (Λ) . For this reason, we see that p h −p a L 2 (Λ) depends significantly on the ratio R/h, and that R < h leads to a poor approximation ofp.
In summary, we see here that the standard FE method has difficulty resolving the pressure differencep − 2πp when R < h, due to the fact thatp is then poorly approximated. This can further pollute the approximations of both the well and reservoir pressure. Explicitly subtracting the singularity in p, which results in the Singularity Removal Based FE described by (56), restores optimal convergence rates for the reservoir pressure p, and improves the robustness of the method with respect to a small well radius R.