A new and consistent well model for one-phase flow in anisotropic porous media using a distributed source model

TL;DR

This paper introduces a novel well model for one-phase flow in anisotropic porous media, utilizing distributed source terms to improve accuracy and robustness across various orientations and anisotropy ratios.

Contribution

The paper presents a new, consistent well model based on distributed sources, independent of discretization, with analytical derivation and numerical validation for anisotropic media.

Findings

Model is robust to well and permeability orientation

Increased accuracy over Peaceman-type models

Independent of discretization method

Abstract

A new well model for one-phase flow in anisotropic porous media is introduced, where the mass exchange between well and a porous medium is modeled by spatially distributed source terms over a small neighborhood region. To this end, we first present a compact derivation of the exact analytical solution for an arbitrarily oriented, infinite well cylinder in an infinite porous medium with anisotropic permeability tensor in R3 , for constant well pressure and a given injection rate, using a conformal map. The analytical solution motivates the choice of a kernel function to distribute the sources. The presented model is independent from the discretization method and the choice of computational grids. In numerical experiments, the new well model is shown to be consistent and robust with respect to rotation of the well axis, rotation of the permeability tensor, and different anisotropy ratios. Finally, a comparison with a Peaceman-type well model suggests that the new scheme leads to an increased accuracy for injection (and production) rates for arbitrarily-oriented pressure-controlled wells.

Tables

Table 1. Table 1 : Convergence rates for E q subscript 𝐸 𝑞 E_{q} for different anisotropy ratios α 𝛼 \alpha .

	$h_{\text{max}}$
$\alpha$	$17.32m$	$8.66m$	$4.33m$	$2.17m$
1	2.0545	2.0724	1.9454	-
10	1.7715	2.0184	2.0763	-
50	1.5904	1.9747	2.0925	-
100	1.5970	1.9666	2.1218	-

Equations

$$-\nabla \cdot \left(\frac{\rho }{\mu }k\nabla p\right)=q{\Phi }_{\Lambda }\text{in }\Omega ,$$

$$\int {\text{displaylimits}}_0^{L_i}\int {\text{displaylimits}}_0^{2\pi }\int {\text{displaylimits}}_0^{\varrho (s)}{\Phi }_{{\Lambda }_i}r\text{d}r\text{d}\theta \text{d}s=L_i,$$

$$p(r)=-\frac{\mu }{k\rho }\frac{q}{2\pi }\ln r+C.$$

$$p_{\omega }=p(r_{\omega })=-\frac{\mu }{k\rho }\frac{q}{2\pi }\ln r_{\omega }+C\Rightarrow C=p_{\omega }+\frac{\mu }{k\rho }\frac{q}{2\pi }\ln r_{\omega },$$

$$p(r)=-\frac{\mu }{k\rho }\frac{q}{2\pi }\ln \left(\frac{r}{r_{\omega }}\right)+p_{\omega }.$$

$$q=2\pi r_{\omega }\frac{\rho k}{\mu }\frac{(p_{\omega }-p_{\delta })}{r_{\omega }\ln \left(\frac{\delta }{r_{\omega }}\right)}$$

$${\Phi }^{\text{const}}(r)=\{\begin{array}{ll}\frac{1}{\pi {\varrho }^2}& r\le \varrho ,\\ 0& r>\varrho ,\end{array}$$

$$p(r)=\{\begin{array}{ll}-\frac{\mu }{k\rho }\frac{q}{2\pi }\left[\frac{r^2}{2{\varrho }^2}+\ln \left(\frac{\varrho }{r_{\omega }}\right)-\frac{1}{2}\right]+p_{\omega }& r\le \varrho ,\\ -\frac{\mu }{k\rho }\frac{q}{2\pi }\ln \left(\frac{r}{r_{\omega }}\right)+p_{\omega }& r>\varrho .\end{array}$$

$$q=2\pi \frac{\rho k}{\mu }(p_{\omega }-p_0)\Xi ,\text{with}\Xi =\frac{(p_{\omega }-p_{\delta })}{(p_{\omega }-p_0)}\frac{1}{\ln \left(\frac{\delta }{r_{\omega }}\right)}$$

$$p_0=-(p_{\omega }-p_0)\Xi \left[\ln \left(\frac{\varrho }{r_{\omega }}\right)-\frac{1}{2}\right]+p_{\omega },$$

$$\Xi ={\left[\ln \left(\frac{\varrho }{r_{\omega }}\right)-\frac{1}{2}\right]}^{-1}.$$

$$-\nabla \cdot \left(\frac{\rho }{\mu }K\nabla p\right)=0\text{in }\hat{\Omega },$$

$$K=QD{Q}^T,$$

$$U:{\mathbb{R}}^3\to {\mathbb{R}}^3,\mathbf{x}\mapsto \mathbf{u}=\tilde{S}\mathbf{x},$$

$$-{\nabla }_u\cdot \left(\frac{\rho }{\mu }{k}_I\nabla {}_{u}p\right)=0\text{in }{\hat{\Omega }}_u=U(\hat{\Omega }),$$

$${\mathbf{x}}^T\Psi \mathbf{x}=r_{\omega }^2,\Psi =I-\mathbf{\psi }{\mathbf{\psi }}^T.$$

$${\hat{\mathbf{v}}}^{T}E\hat{\mathbf{v}}={\hat{\mathbf{v}}}^T{P}^{T}R{\tilde{S}}^{-1}\Psi {\tilde{S}}^{-1}{R}^{T}P\hat{\mathbf{v}}=r_{\omega }^2,{\tilde{S}}^{-1}{R}^{T}P\hat{\mathbf{v}}=\mathbf{x}$$

$$\quad R=2\frac{(\boldsymbol{e}_{3}+\boldsymbol{\psi}^{\prime})(\boldsymbol{e}_{3}+\boldsymbol{\psi}^{\prime})^{T}}{(\boldsymbol{e}_{3}+\boldsymbol{\psi}^{\prime})^{T}(\boldsymbol{e}_{3}+\boldsymbol{\psi}^{\prime})}-I,\quad\boldsymbol{e}_{3}=\left[\begin{array}[]{c}0\\ 0\\ 1\end{array}\right],\quad P=\left[\begin{array}[]{cc}1&0\\ 0&1\\ 0&0\end{array}\right].$$

$$V:{\mathbb{R}}^3\to {\mathbb{R}}^3,\mathbf{u}\mapsto \mathbf{v}=\tilde{R}\mathbf{u}={\hat{R}}^T{R}^T\mathbf{u},$$

$$\hat{R}=\left[\boldsymbol{\nu}_{1}\bigg{|}\boldsymbol{\nu}_{2}\bigg{|}\boldsymbol{\psi^{\prime}}\right]$$

$$Z:{\mathbb{R}}^3\to \mathbb{C},\mathbf{v}\mapsto z=\tilde{Z}\mathbf{v}=[1,i,0]\mathbf{v}=v_1+i{v}_2,$$

$$T:\mathbb{C}\to \mathbb{C},z\mapsto w=z+\sqrt{z-f}\sqrt{z+f},f=\sqrt{a^2-b^2}$$

$$p(w)=p_{\omega }-\frac{\mu }{\rho k_I}\frac{\hat{q}}{2\pi }\ln \left(\frac{|w|}{r_{\circ }}\right)\zeta ,\hat{q}=q\zeta =q\frac{ab}{r_{\omega }^2},$$

$$w=T(Z(V(U(\mathbf{x}))))=T(\tilde{Z}\tilde{R}\tilde{S}\mathbf{x}).$$

$$K_{A}=\left[\begin{array}[]{ccc}1&0&0\\ 0&5&4\\ 0&4&5\end{array}\right]$1\text{\cdot}{10}^{-10}\text{\,}{\mathrm{m}}^{2}$$

$

$$T^{-1}:\mathbb{C}\to \mathbb{C},w\mapsto z=\frac{1}{2}\left(w+\frac{f^2}{w}\right),|w|>f,$$

$$J_{T^{-1}}=\left[\begin{array}[]{cc}\frac{\partial x}{\partial u}&\frac{\partial x}{\partial u}\\ \frac{\partial y}{\partial v}&\frac{\partial y}{\partial v}\end{array}\right]=\left[\begin{array}[]{cc}\eta&\epsilon\\ -\epsilon&\eta\end{array}\right],$$

$${\Delta }_{w}p=\frac{{\partial }^{2}p}{\partial u^2}+\frac{{\partial }^{2}p}{\partial v^2}={\left|\frac{\partial T^{-1}}{\partial w}\right|}^2{\Delta }_{z}p=\left|\det (J_{T^{-1}})\right|{\Delta }_{z}p$$

$${\left|\det (J_{T^{-1}})\right|}^{-1}{\Delta }_{w}p={\Delta }_{z}p,$$

$$\left|\det (J_{T^{-1}})(w)\right|={\left|\frac{\partial T^{-1}}{\partial w}\right|}^2={\left|\frac{\partial z}{\partial u}\right|}^2=\frac{1}{4}\left(1+\frac{f^4-2f^2\Re (w^2)}{|w{|}^4}\right):={\Phi }_J^{-1},$$

Full text

1 Introduction

Well modeling is essential for various engineering applications, as for example reservoir simulation, geothermal energy production or energy storage, where injection or extraction processes strongly influence the flow behavior. Usually the well geometry is not explicitly resolved in the mesh but instead modeled as a line source with given extraction or injection rate. However, this simplified approach introduces singularities, meaning that the logarithmic solution profiles are undefined at the center-line of the well. This leads to a significant deviation between numerical and analytical solution in the near-well region. For a better approximation, locally refined meshes around the wells are needed, which however deteriorate efficiency and is therefore often not suitable for field-scale simulations, especially when multiple wells are present. Similar issues are encountered for the modeling of vascularized biological tissue perfusion [1, 2] or the modeling of root water update [3, 4].

A common approach is the use of well-index-based well models. Such well models aim to find a relation between well rate, bottom hole pressure and numerically calculated pressure (well-block pressure) for each cell (grid-block) that contains the well. In reservoir engineering such a relation is denoted as well index. The first theoretical derivation of a well index for two-dimensional structured uniform grids with isotropic permeability has been presented by Peaceman [5]. He has shown that the well-block pressure differs from the areal averaged analytical pressure and introduced a new relation by using an equivalent well radius. The equivalent well radius is defined as the distance (relative to the well location) at which the analytical and numerical pressures are equal. A generalization for structured two-dimensional non-square grids ($\Delta x\not=\Delta y$) and anisotropic but diagonal permeability tensors has been presented by the same author in [6]. However, these well models are restricted to uniform grids where the well is oriented with one of the grid axes. Furthermore, since the effective well radius relates numerical and analytical pressure values, the well index has to be calculated depending on the used discretization scheme. Thus, Peaceman’s well model is only valid for a cell-centered finite difference scheme with 5-point stencil. Discussion of other discretization schemes can be found in [7], again with the restriction of two-dimensional grids. Enhancements include, among others, three-dimensional slanted wells [8, 9], Green’s functions for the computation of well indices [10, 11, 12], or the singularity subtraction method to obtain smooth solutions in the near-well region [13, 14].

In this work, a new approach for obtaining more accurate source term for a given well bottom hole pressure is presented. The new model is, in contrast to most of the existence methods, independent of the discretization scheme and can be used for general unstructured grids. Additionally, the presented method is not restricted to diagonal tensors and thus works for general anisotropic permeabilities. In Section 2, we derive a well model, initially for isotropic porous media, for which the fluid mass injected by a well is distributed over a small neighborhood around the well, using kernel functions. The derivation follows the idea recently presented in [15], however we herein discuss the case without membrane or casing. The model yields a pressure solution without singularity, from which the source term can be reconstructed using a relation found with the analytical solution for the case of an infinite well in an infinite medium. The model generalizes to more complex problems due to the superposition principal valid for the Laplace operator, a linear operator [15]. In Section 3, the model is generalized to porous media with general anisotropic permeabilities, based on an analytical solution constructed in Section 3.1 using a series of coordinate transformations. We shown that the general model reduces to the model derived in Section 2 for isotropic permeabilities. Finally, the new well model is analyzed with several numerical experiments in Section 5. The results indicate that the model is consistent for different anisotropy ratios, robust with respect to rotations of the well relative to the computational grid, and to rotations of an anisotropic permeability tensor. A comparison with a Peaceman-type well model in a setup with a K-orthogonal grid and an embedded slanted well suggests that the new model more accurately approximates the fluid exchange between well and rock matrix.

2 A well model with distributed source for isotropic media

First, we derive a well model with distributed source for porous media with isotropic permeability tensor, not including a well casing. The derivation follows [15], where a casing is included in form of a membrane in a different but related application, that is modeling fluid exchange between the vascular system and the embedding biological tissue. Stationary single-phase flow around a well with radius $r_{\omega}$, in an isotropic porous medium with permeability $k$, can be described by the following flow equation

[TABLE]

where $p$ is the fluid pressure, $\rho$ the fluid density and $\mu$ its dynamic viscosity. Denoted by $\Phi_{\Lambda}$ is a set of kernel functions $\Phi_{\Lambda_{i}}$ that distribute $q$ ($\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-1}$) over a small tubular support region, $\mathcal{S}(\Phi_{\Lambda_{i}})$, with radius $\varrho(s)$, around a well segment $i$, such that $\Phi_{\Lambda_{i}}=0$ outside the support region. We choose kernel functions $\Phi_{\Lambda_{i}}(s)$ with the property

[TABLE]

where $r$, $\theta$, $s$ are the radial, angular, and axial coordinate in a segment-local cylinder coordinate system, and $L_{i}$ is the length of segment $i$.

Assume a radially symmetric zone (distance $\delta>r_{\omega}$ from center-line) around the well, with $p_{\delta}$ denoting the pressure at distance $\delta$ from the well center-line, and constant fluid density and viscosity. Then, the pressure for $r_{\omega}<r<\delta$ is described by the analytical solution

[TABLE]

The constant $C$ is determined by fixing a well pressure, $p_{\omega}$,

[TABLE]

so that

[TABLE]

Consequently, the source term can be expressed in terms of $p_{\omega}$ and $p_{\delta}$ as

[TABLE]

We choose a simple kernel function which regularizes the pressure solution for $r\leq\varrho$,

[TABLE]

where $\varrho\leq\delta$. The pressure for $r<\rho$ can be obtained by integration from Eq. 2.1, yielding

[TABLE]

Figure 1 graphically explains the most important symbols introduced in this section. As the regularized pressure can be evaluated at the well center-line we can reformulate Eq. 2.5,

[TABLE]

where $\Xi$ is the so called flux scaling factor. The flux scaling factor can be also expressed independent of the pressure. To this end, Eq. 2.7 is evaluated at $r=0$, so that $p_{\omega}$ is expressed in terms of $p_{0}$,

[TABLE]

where $q$ was replaced by inserting Eq. 2.8. It directly follows from Eq. 2.9 that

[TABLE]

3 A well model with distributed source for anisotropic media

In the following section, the developed well model is extended for porous media with anisotropic permeability. In Section 3.1, we derive an analytical solution for one-phase flow around an infinitely long cylindrical well embedded in an infinite porous domain in $\mathbb{R}^{3}$. This derivation motivates the choice of a suitable kernel function for anisotropic problems, presented in Section 3.3.

3.1 Analytical solution for anisotropic permeability and slanted well

In the following section, we derive an analytical solution for one-phase flow around an infinite cylindrical well $\Gamma$ with radius $r_{\omega}$ in an infinite porous domain $\hat{\Omega}=\mathbb{R}^{3}\setminus\Gamma$ with anisotropic, homogeneous permeability. We assume, without loss of generality, that the well axis passes through the origin of the Cartesian coordinate system, and denote by $\boldsymbol{\psi}$ a unit vector parallel to the well signifying the well orientation. We seek an analytical expression for the hydraulic pressure $p$ such that

[TABLE]

for a constant well pressure $p_{\omega}$ in $\mathrm{Pa}$ and some specific pumping rate $q$ in $\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-1}$ given on $\partial\Gamma$. The total mass flow over the boundary of a well segment of length $L$ is thus given by $Q=qL$.

From thermodynamic constraints, $K$ is a positive definite and symmetric, second-order tensor field. Hence, $K$ can be decomposed such that

[TABLE]

where $D=\operatorname{diag}(\lambda_{1},\lambda_{2},\lambda_{3})$ is a diagonal matrix composed of the eigenvalues $\lambda_{i}$ of $K$, $Q=[\boldsymbol{\nu}_{K,1}|\boldsymbol{\nu}_{K,2}|\boldsymbol{\nu}_{K,3}]$ is a rotation matrix with the corresponding eigenvectors as columns, and $A^{T}$ denotes the transpose of a matrix $A$. Further useful properties derived from the decomposition are $\operatorname{det}(K)=\lambda_{1}\lambda_{2}\lambda_{3}$, where $\operatorname{det}(A)$ denotes the determinant of $A$, and $K^{n}=Q\Lambda^{n}Q^{T}$, where $D^{r}=\operatorname{diag}({\lambda_{1}^{r},\lambda_{2}^{r},\lambda_{3}^{r}})$, $r\in\mathbb{R}$.

It is well known that the anisotropic one-phase flow problem can be transformed to an isotropic problem using a coordinate transformation [9, 16, 17, 6, 18]

[TABLE]

with the stretching matrix $\tilde{S}=k_{I}^{1/2}K^{-1/2}$, where $k_{I}$ is an arbitrary scalar constant, that we choose as $k_{I}=\operatorname{det}(K)^{-1/3}$ (cf. [9]), rendering the transformation isochoric. The transformation $\boldsymbol{u}=\tilde{S}\boldsymbol{x}$ deforms the well cylinder such that a cross-section orthogonal to the transformed well direction is elliptical. The solution to the isotropic problem

[TABLE]

is identical on two parallel planes perpendicular to the transformed (normalized) well direction, $\boldsymbol{\psi}^{\prime}=\tilde{S}\boldsymbol{\psi}||\tilde{S}\boldsymbol{\psi}||^{-1}$. This motivates the rotation of the coordinate system such that the first and second axis are aligned with the major and minor axis of the well-bore ellipse and third axis is aligned with $\boldsymbol{\psi}^{\prime}$. To determine the corresponding rotation matrix $\tilde{R}$, we need to characterize this well-bore ellipse. The well cylinder in $x$-coordinates is given by

[TABLE]

After stretching, the coordinate system can be rotated with the rotation matrix $R$ so that the third axis is aligned with the well direction. Then, projecting into the plane perpendicular to the well direction yields the well-bore ellipse equation

[TABLE]

in $\hat{v}$-coordinates, where

[TABLE]

The rotation matrix $R$ can be derived using Rodrigues’ rotation formula as shown in Appendix A. The length of the major and minor ellipse axis are found as $a=r_{w}\gamma_{1}^{-1/2}$ and $b=r_{w}\gamma_{2}^{-1/2}$, where $\gamma_{i}$ are the eigenvalues of $E$, and the axis orientations are given by $\boldsymbol{\nu}_{1}=P\hat{\boldsymbol{\nu}}_{E,1}$, $\boldsymbol{\nu}_{2}=P\hat{\boldsymbol{\nu}}_{E,2}$, where $\hat{\boldsymbol{\nu}}_{E,i}$ denote the corresponding eigenvectors of $E$. We assume that the eigenvalues and eigenvectors are sorted such that $a\geq b$, and oriented such that $\boldsymbol{\psi^{\prime}}=\boldsymbol{\nu}_{1}\times\boldsymbol{\nu}_{2}$. Finally the desired rotation is given by

[TABLE]

where

[TABLE]

is rotating about the well direction axis such that the coordinate system is aligned with the principal ellipse axes.

Following the derivations from above, we now have to solve a two-dimensional isotropic Laplace problem with boundary conditions prescribed on an ellipse. To this end, we note that the transformation of a harmonic function $f$ (a function satisfying Laplace’s equation $\Delta f=\nabla\cdot\!\nabla\!f=0$) with a conformal (angle-preserving) mapping yields another harmonic function [19] (see Appendix B). Using a Joukowsky transformation, a conformal mapping well-known from aerodynamics [20], the isotropic problem with a well with elliptic cross-sections, can be transformed to an isotropic problem with circular cross-sections [16]. Transforming into the complex plane (parametrizing the well-bore ellipse plane)

[TABLE]

the (inverse) Joukowsky transformation

[TABLE]

transforms elliptic isobars into circular isobars, where $a$ and $b$, $a\geq b$, are the major and minor axis of the well-bore ellipse, as derived above. In particular, the well-bore ellipse (where $p=p_{\omega}$) is mapped onto a circle with radius $r_{\circ}=a+b$. Finally, in the new coordinate system we find the (now) radially symmetric analytical solution to problem Eq. 3.1

[TABLE]

where the source scaling factor $\zeta$ is necessary to recover the original source $q$ on $\partial\Gamma$. This can be derived from simple geometric considerations as shown in Appendix C. The $w$ corresponding to some $x\in\hat{\Omega}$ in original coordinates is obtained by using all above-mentioned transformations after each other as follows

[TABLE]

Such a solution for a slanted well ($30\text{\,}\mathrm{\SIUnitSymbolDegree}$ with respect to vertical axis) and anisotropic permeability tensor

[TABLE]

is visualized in Fig. 2.

3.2 Properties of the conformal mapping

To construct a suitable kernel function for anisotropic problems, we first have a closer look at the properties of the employed Joukowsky transformation. The effect of the mapping $T$, Eq. 3.11, is shown in Fig. 3. Points on the exterior of a line on the real axis between $f$ and $-f$ are mapped onto the exterior of a circle with radius $f$. The ellipse with major axis $a$ and minor axis $b$ is mapped onto a circle with radius $r_{\circ}=a+b$. Going further away from the well, the deformation due to the mapping is less and less pronounced. This matches the expectations for the physical flow problem, since isobars at large distance from an elliptical well-bore become increasingly circular in isotropic media.

The inverse transformation is given by

[TABLE]

where the restriction on $|w|$ is necessary to obtain a one-to-one mapping. The transformation $T^{-1}$ can equally be interpreted as an $\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ mapping. The Jacobian of the transformation $T^{-1}$ for $z=x+iy$ and $w=u+iv$ has the form

[TABLE]

which follows from the the Cauchy–Riemann equations [21]. Since the transformation can be viewed as the composition of a scaling and a rotation, it is angle-preserving. The transformation is associated with a spatially dependent volume deformation characterized by the determinant of $J_{T^{-1}}$. Furthermore, it can be shown that the Laplace operator behaves as follows under the transformation $z=T^{-1}(w)$,

[TABLE]

by computing the derivative of the real and the imaginary part of $p$ separately, applying the chain rule and the Cauchy–Riemann equations, as shown for completeness in Appendix B. From Eq. 3.17 follows that

[TABLE]

for the transformation $w=T(z)$. The determinant can be explicitly computed, using Eq. 3.17 and complex differentiation (shown in Appendix D) as

[TABLE]

where $\Re(w^{2})$ is real part of $w^{2}$, and $|w|$ the absolute value of $w$. We note that $\Phi_{J}$ quickly converges to the value $4$ with increasing $|w|$, that is with increasing distance from the well. The function $\Phi_{J}^{-1}$ is plotted in Fig. 4 exemplarily for $f=1.2$.

3.3 A kernel function for anisotropic media

Instead of excluding the well domain $\Gamma$ from $\Omega=\mathbb{R}^{3}$ and modeling infiltration or extraction by a flux boundary condition, we will now model the action of the well on the flow field by a spatially distributed source term, as presented for the isotropic problem,

[TABLE]

From the above derivations, we know that solving Eq. 3.20 in $w$-coordinates is straight-forward. Hence, we choose kernel functions in $w$-coordinates and then transform to $x$-coordinates so that the pressure solution satisfies Eq. 3.20. Motivated by the properties of the Joukowsky transform (see Fig. 3), we choose a local kernel that is constant on the annulus with inner radius $f<\varrho_{i}\leq r_{\circ}$ and outer radius $\varrho_{o}>r_{\circ}$,

[TABLE]

In $w$-coordinates, we can find a solution to the problem

[TABLE]

for a given constant well pressure $p_{\omega}$, $\hat{q}=q\zeta$ and constant density and viscosity.

By means of integration (cf. [15]), we get

[TABLE]

where $\xi^{2}=\varrho_{o}^{2}-\varrho_{i}^{2}$. This shows that outside the kernel support region ($|w|>\varrho_{o}$), we obtain the exact analytical solution derived in Section 3.1. Moreover, the source term can be reformulated, cf. [15],

[TABLE]

where $p_{0}:=p(|w|=0)=p(|w|=\varrho_{i})$ is the fluid pressure evaluated on the well center-line. Note that for $f=0$ and $\varrho_{i}=f=0$, the isotropic solution with for a circular constant kernel (Eqs. 2.7 and 2.10) is obtained. From the transformation of the Laplace operator, Eq. 3.18, we see that the problem

[TABLE]

with altered kernel function $\Phi_{\Lambda}=\Phi_{A}\Phi_{J}$ is equivalent to Eq. 3.22.

The transformation $T^{-1}$ changes shape of the kernel support $\mathcal{S}(\Phi_{A})$ from an annulus to an ellipse $E_{\Phi,v}$. Inverting $Z$ extrudes the solution along the well center-line, and inverting the rotation and stretch described by $V$ and $U$ results in a kernel support region in the shape of an elliptic cylinder. Moreover, each ellipse $E_{\Phi,v}$ with normal vector $\boldsymbol{e}_{3}$ is transformed to an ellipse $E_{\Phi,x}(s)$, that is the intersection the elliptic cylinder with a plane with the normal vector $\boldsymbol{n}_{E_{x}}=\tilde{S}\tilde{R}^{T}\boldsymbol{e}_{3}$, centered at $s$ on the well center-line. The transformation and the normal vector $\boldsymbol{n}_{E_{x}}$ are visualized in Fig. 5. We note that if none of the principal axes of the permeability tensor are aligned with the well direction, $\boldsymbol{n}_{E_{x}}$ is not parallel to the well direction $\boldsymbol{\psi}$ in $x$-coordinates. The integral of the right-hand side of Eq. 3.20 for a well segment $\Lambda_{i}$ with length $L_{i}$ is equal to the integral over the kernel support $\mathcal{S}(\Phi_{\Lambda,i})$ which has the shape of the elliptic cylinder given by

[TABLE]

Using $\hat{q}=q\zeta$, and exploiting that $k_{I}$ was chosen such that $\operatorname{det}(\tilde{R}\tilde{S})=1$, it can be shown that

[TABLE]

where $0\leq\hat{s}\leq\hat{L}_{i}$ is a local coordinate along the transformed well direction, and the last equality is proven in Appendix C. This is the desired property of the kernel function for the anisotropic case corresponding to Eq. 2.2 for the isotropic case.

4 Numerical method

We discretize Eq. 3.20 using a cell-centered finite volume method with multi-point flux approximation (MPFA) [22]. The domain $\Omega$ is decomposed into control volumes $K_{\Omega}\in\Omega_{h}$ such that the computational mesh $\Omega_{h}$ is a discrete representation of $\Omega$. Furthermore, each control volume boundary, $\partial K_{\Omega}$, can be split into a finite number of faces $\sigma\subset\partial K_{\Omega}$, such that $\sigma=K_{\Omega}\cap L_{\Omega}$, with $L_{\Omega}$ denoting a neighboring control volume. Integrating Eq. 3.20 over a control volume $K_{\Omega}$ and applying the Gauss divergence theorem on the left hand side yields

[TABLE]

where $\boldsymbol{n}_{K_{\Omega},\sigma}$ is the unit outward-pointing normal on face $\sigma\subset\partial K_{\Omega}$. The exact fluxes are approximated by numerical fluxes

[TABLE]

which are computed using the MPFA-O method described in [22]. The discrete source term is computed as

[TABLE]

where $Q_{\mathcal{I}}$ is a numerical approximation of the source term integral over the intersection $\mathcal{I}=K_{\Omega}\cap\Lambda$,

[TABLE]

and $\mathcal{S}(\Phi_{\Lambda,\mathcal{I}})$ is the kernel support associated with $\mathcal{I}$ as depicted in Fig. 6. In summary, the discrete form of Eq. 4.1 is

[TABLE]

We note that due to the dependency of $Q_{\mathcal{I}}$ on $p_{0}$ the proposed method is non-local in the sense that non-neighbor cells $M_{\Omega}\in\Omega_{h}$ (where $M_{\Omega}\cap K_{\Omega}$ is the empty set or a single point) may have an associated degree of freedom that depends on the degree of freedom of $K_{\Omega}$.

4.1 Kernel integration

The kernel integral in Eq. 4.3

[TABLE]

is not easily approximated with a quadrature rule, since the intersection $K_{\Omega}\cap\mathcal{S}(\Phi_{\Lambda,\mathcal{I}})$, that is the intersection of an elliptic cylinder with for example a hexahedron is difficult to compute. However, we use the same idea as in [15], and remark that the integral over the entire support $\mathcal{S}(\Phi_{\Lambda,\mathcal{I}})$ is known exactly; see Eq. 3.27. Hence, the integration problem can be reformulated as the distribution of the known integral over all intersected control volumes $K_{\Omega}$ weighted with the respective support volume fractions. Following [15], we create $n_{\mathcal{I}}$ integration points $\boldsymbol{x}_{i}\in\mathcal{S}(\Phi_{\Lambda,\mathcal{I}})$ with known volume elements $V_{i}$ of similar size and shape, so that

[TABLE]

Computing the weights for each cell $K_{\Omega}$ is a pre-processing step that only has to be done once for each computational mesh and well geometry.

5 Numerical experiments and discussion

We present numerical experiments using the presented method in different setups. All experiments are conducted with constant fluid density $\rho=$1000\text{,}\mathrm{kg}\text{,}{\mathrm{m}}^{-3} and viscosity $\mu=$1\text{\cdot}{10}^{-3}\text{\,}\mathrm{Pa}\text{\,}\mathrm{s}. The well pressure is constant, $p_{\omega}=$1\text{\cdot}{10}^{6}\text{,}\mathrm{Pa}, and the well radius is $r_{\omega}=$0.1\text{\,}\mathrm{m} if not specified otherwise. The permeability tensor is given as

[TABLE]

where

[TABLE]

are rotation matrices rotating vectors about $\boldsymbol{e}_{1}$, $\boldsymbol{e}_{2}$ by the rotation angle $\gamma_{1}$, $\gamma_{2}$, respectively, and $\alpha$ is a given dimensionless K-anisotropy ratio $\alpha=\frac{K_{33}}{K_{11}}=\frac{K_{33}}{K_{22}}$. The domain $\Omega_{0}=[-100,100]\times[-100,100]\times[-50,150]~{}${\mathrm{m}}^{3} is split in two regions, $\Omega=[-100,100]\times[-100,100]\times[0,100]~{}${\mathrm{m}}^{3} and $\Omega_{D}=\Omega_{0}\setminus\Omega$. The well center-line $\Lambda$ is given by the line through the origin and $\boldsymbol{\psi}=R_{1}(\beta_{1})R_{2}(\beta_{2})\boldsymbol{e}_{3}$, where $R_{1}$, $R_{2}$ are given in Eq. 5.2 and $\beta_{1}$, $\beta_{2}$, are rotation angles. The analytical solution for all cases is given in Eq. 3.23, $q=$1\text{,}\mathrm{kg}\text{,}{\mathrm{s}}^{-1}\text{,}{\mathrm{m}}^{-1}$$, and $L=|\Lambda\cap\Omega|$ (in $\mathrm{m}$). For all setups the inner kernel radius is chosen as $\varrho_{i}=f$. In all of $\Omega_{D}$ and on the boundary $\partial\Omega$ the analytical solution is enforced by Dirichlet constraints, modeling the infinite well. The computational mesh $\Omega_{h}$ is a structured grid composed of regular hexahedra $K_{\Omega}$. Furthermore, we define two error measures.

[TABLE]

is the relative discrete $L^{2}$-norm of the pressure, where $p_{e,\boldsymbol{x}_{K_{\Omega}}}$ is the exact pressure evaluated at the cell centroid $\boldsymbol{x}_{K_{\Omega}}$ and $p_{K_{\Omega}}$ the discrete numerical cell pressure, and

[TABLE]

is the relative discrete $L^{2}$-norm of the source term, where $|\mathcal{I}|=|K_{\Omega}\cap\Lambda|$ is the length of the intersection of cell $K_{\Omega}$ and the well center-line $\Lambda$, $Q_{\mathcal{I}}$ is the discrete source term given in Eq. 4.4. All setups are implemented in DuMux [23], an open-source porous media simulator based on Dune [24, 25].

5.1 Grid convergence for different anisotropy ratios

In the first numerical experiment grid convergence is investigated for different anisotropy ratios $\alpha$. To this end, $h_{\text{max}}:=\operatorname{max}\limits_{K_{\Omega}\in\Omega_{h}}h_{K_{\Omega}}$, where $h_{K_{\Omega}}$ is defined as the maximum distance between two vertices of the cell $K_{\Omega}$. Starting at a grid resolution for $\Omega_{h}$ of $20\times 20\times 10$ cells ($h_{\text{max}}=10\sqrt{3}~{}$\mathrm{m}$$), the grid is refined uniformly. Figure 7 shows the errors $E_{p}$ and $E_{q}$ for different grid resolutions and values of $\alpha$, for $\beta_{1}=\beta_{2}=20^{\circ}$ and $\gamma_{1}=\gamma_{2}=-20^{\circ}$, so that $K$ is a full tensor and none of the principal axis of $K$ is aligned with the well direction. For all $\alpha$, the method shows second order convergence for the pressure in the given norm, as expected for the MPFA-O method [26] (super convergence at cell centers). The source term $q$ is a linear function of the pressure $p$ and also exhibits second order convergence. We note that the errors for different $\alpha$ are not directly comparable since the analytical solution for $p$ changes with $\alpha$, although $q$ is constant.

However, the convergence rates are shown to be independent of $\alpha$ with increasing grid resolution. The convergence rates for $E_{q}$ (slope of the lines in Fig. 7) are presented in Table 1. It can be seen that rates for large grid cells and large $\alpha$ are slightly smaller. This is because the kernel support is still under-resolved by the computational grid. For example, for $\alpha=100$, the kernel ellipse in $x$-coordinates has major and minor axis of $a_{x}\approx$55.9\text{,}\mathrm{m}, $b_{x}\approx$5.6\text{\,}\mathrm{m}, respectively, while $h_{\text{max}}\approx$17.32\text{,}\mathrm{m}$$ for the lowest grid resolution.

5.2 Influence of the outer kernel radius $\varrho_{o}$

In [15], it is suggested that increasing the kernel support region (increasing $\varrho_{o}$), has a similar effect on $E_{q}$ as refining the grid. However, the pressure solution is then regularized in a larger region, so that there is a trade-off between the accuracy of the source term and the accuracy of the pressure field with respect to the unmodified problem ($\varrho_{o}\rightarrow\varrho_{i},\Phi_{\Lambda}\rightarrow\delta_{\Lambda}$). However, every discrete cell $K_{\Omega}$ can be also interpreted as a kernel support region, such that the choice of $\Phi_{\Lambda}$ enables us to better control the discretization error as soon as $\mathcal{S}(\Phi_{\Lambda},\mathcal{I})$ becomes larger than $K_{\Omega}$.

As shown in [6] and Fig. 4, isobars become circular, in the transformed domain $U(\Omega)$, with increasing distance to the well. Therefore, a reasonable simplification is $\Phi_{J}\approx 4$ if $\varrho_{o}\gg r_{\omega}$. This is completely analogous to the assumption of circular isobars in [6], where an estimate of the error introduced by the assumptions is given for the two-dimensional case.

In the following numerical experiment, we step-wise increase the kernel radius $\varrho_{o}$, for the same $20\times 20\times 10$ grid. This is done once for the case, where $\varrho_{o}\gg r_{\omega}$ and for the case for which $\varrho_{o}$ is only slightly larger than $r_{\omega}$. Furthermore, $\beta_{1}=\beta_{2}=20^{\circ}$ and $\gamma_{1}=\gamma_{2}=-20^{\circ}$. The results are shown in Fig. 8. First, it can be seen that doubling $\varrho_{o}$ leads to a $4$-times smaller error $E_{q}$. This can be explained by the fact that the larger the kernel, the more grid cells resolve the kernel support, and the better is the approximation of $p_{0}$. Furthermore, the result is consistent with the results in [15]. Moreover, Fig. 8 suggests that for $\varrho_{o}\gg r_{\omega}$ the simplification of the kernel function ($\Phi_{J}\approx 4$) is not visible in $E_{q}$, while for kernel radii slightly smaller than the well radius, the simplification increases $E_{q}$ by an order of magnitude in comparison to the case using the exact kernel function as derived in Section 3.3. The results show that the presented method is also applicable in cases where the grid resolution is very close to the well radius. An adaption of the presented method for other applications, such as the simulation of flow in vascularized tissue, where such ratios of vessel radius to cell size are typical, cf. [15], is therefore well-conceivable.

5.3 Robustness with respect to rotation

In the following numerical experiment, we use a single computational mesh with a given resolution for $\Omega_{h}$: $20\times 20\times 10$. First, the well direction is fixed, and the permeability tensor is rotated by varying $\gamma_{1}$ and $\gamma_{2}$. Then the permeability tensor is fixed and the well is rotated by varying $\beta_{1}$ and $\beta_{2}$. The results are shown in Fig. 9.

It can be seen that the presented well model is rather robust with respect to rotations. Possible effects influencing the approximation error $E_{q}$, include the different quality of the kernel integral for different angles with respect to the grid axes, and differences in the flux approximation quality of the MPFA-O method depending on the face co-normal $\boldsymbol{d}_{K_{\Omega},\sigma}=K\boldsymbol{n}_{K_{\Omega},\sigma}$. Additionally, for different well angles the number and size of intersections $K_{\Omega}\cap\Lambda$ can have an influence on the discrete error.

5.4 Comparison with a Peaceman-type well model

In particular for petroleum engineering applications, commercial codes typically use Peaceman-type well models [27, 28]. In [6], Peaceman extended his well known well-index-based well model for anisotropic diagonal permeability tensors and non-cubic but structured rectangular grids. The discrete source term in a computational cell $K_{\Omega}$ is approximated by

[TABLE]

where $\Delta x$ and $\Delta y$ are the horizontal dimensions of the cell containing the well, $L_{K_{\Omega}}=|\mathcal{I}|$ the length of the well segment contained in $K_{\Omega}$, and $\gamma$ the Euler–Mascheroni constant. The Peaceman model has several known limitations. Its derivation only applies to K-orthogonal structured grids, where the well is oriented along one of the grid axes, and perfectly horizontally centered within a vertical column of computational cells $K_{\Omega}$. Furthermore, the derivation is specific to cell-centered finite difference schemes with 5-point stencil. Moreover, computational cells may have to be significantly larger than the well radius (depending on the degree of anisotropy) for optimal accuracy. The Peaceman model has been generalized for slanted wells with arbitrary orientation, for example in [8]. The Alvestad model [8] has been adapted for finite volumes, for example in [9] (formula given in Appendix E, subsequently referred to as pm well model). Such extensions usually constitute a reasonable directional weighting of the original Peaceman model but are not directly derived from the mathematical analysis of the underlying problem [9].

The herein presented model has none of the above-mentioned limitations. In particular, the presented model is valid for arbitrary positive definite and symmetric permeability tensors, unstructured grids, and is independent of the discretization scheme. Moreover, the presented model is consistent and we show grid convergence in the numerical experiments in Section 5.1. However admittedly, the Peaceman-type models are cell-local, thus computationally cheaper and easier to implement.

Several limitations of the Peaceman well models make it difficult to fairly compare it with our new model. For cases for which all assumptions of Peaceman are valid, our numerical studies (not shown here) suggest that the Peaceman well model is generally superior to the presented model with distributed sources. This is because it takes the analytical solution as well as the spatial discretization method into account. For cases where some assumptions are violated, for example off-center wells or slanted wells, it is difficult to construct cases where the analytical solution is readily constructed but does not feature a singularity on the boundary. Our preliminary numerical studies for such cases (for example the slanted well case in Section 5.1 without rotation of the permeability tensor) show large deviations ($>$10\text{,}\mathrm{\char 37\relax}$$ error in total source term) from the analytical solution for the pm well model. However, these errors may be distorted by errors made in the discrete approximation of the singular boundary condition, where the well intersects the boundary. Finally, for the general case of unstructured grids, simplex grids, and full permeability tenors it is unclear how to apply the original Peaceman model. However, we know that the presented method is consistent (at least for a single straight well), and thus, the numerical solution converges to the exact solution with grid refinement. Therefore, we expect that the numerical solution on a very fine grid using the distributed source model is a reasonable reference solution.

We compare our model to the pm well model in a numerical experiment. The computational domain $\Omega=[-50,50]\times[-100,100]\times[0,100]~{}${\mathrm{m}}^{3} contains a slanted straight well $\Lambda$ with end points at $\boldsymbol{x}_{\Lambda,1}=[-20,-50,25]^{T}~{}$\mathrm{m}, $\boldsymbol{x}_{\Lambda,2}=[20,50,75]^{T}~{}$\mathrm{m}. The permeability tensor is a diagonal tensor $K(\gamma_{1},\gamma_{2})$, with $\gamma_{1}=0^{\circ},\gamma_{2}=90^{\circ},\alpha=0.1$. The structured cube grid $\Omega_{h}$ is successively, uniformly refined starting with $10\times 20\times 10$ cells ($h_{\text{max}}\approx$17.32\text{\,}\mathrm{m}). The well radius is $r_{\omega}=$0.1\text{,}\mathrm{m} ($\Delta x/r_{\omega}=200$ for the coarsest grid). The kernel support region (chosen as $\varrho_{o}/r_{\omega}=100$) only extends over few cells in the coarsest grid, so that the regularization effect is minimized. On the boundary $\partial\Omega$, we specify Neumann no-flow boundary conditions, that is $(K\nabla\!p)\cdot\boldsymbol{n}=0$, except for the planes perpendicular to the $x_{1}$-axis, where the Dirichlet boundary condition $p_{D}(x_{2}=-100)=$1\text{\cdot}{10}^{5}\text{\,}\mathrm{Pa}, $p_{D}(x_{2}=100)=$3\text{\cdot}{10}^{5}\text{,}\mathrm{Pa} are enforced. The reference solution is computed with $160\times 320\times 160$ cells ($h_{\text{max}}\approx$1.08\text{\,}\mathrm{m}). The computational domain with pressure iso-surfaces of the reference solution are shown in Fig. 10.

In Fig. 11, the relative integral source error

[TABLE]

with respect to the reference solution $Q_{\text{ref}}$ is shown for grids with different refinement. In a variant of the distributed source model (ds), the extent of the kernel support is adapted to the grid size. This is to keep the regularization effect of the kernel function minimal in order to get, in addition to a good approximation of the source term, a better approximation of the pressure solution close to the well. For $\varrho_{o}/r_{\omega}=100$, the extent of the kernel ellipse $E_{\Phi,x}$ is given by its major and minor axes, $16.12\text{\,}\mathrm{m}$ and $12.54\text{\,}\mathrm{m}$. For the reference solution this extent is kept constant with grid refinement. While this ensures a very good approximation of the source term, the pressure solution is regularized in a larger neighborhood of the well. In the variant, the kernel support is adapted proportional to $h_{\text{max}}$, so that for the finest grid shown in Fig. 11 ($80\times 160\times 80$ cells, $h_{\text{max}}\approx$2.17\text{,}\mathrm{m}$$), the $E_{\Phi,x}$ major and minor axes measure $2.01\text{\,}\mathrm{m}$ and $1.57\text{\,}\mathrm{m}$.

It is evident that the numerical solution for the distributed source model converges to the reference solution. More importantly, the relative error is small ($<$0.5\text{,}\mathrm{\char 37\relax}) even for the coarsest grid. In comparison, the difference to the Peaceman-type model is large ($>$5\text{\,}\mathrm{\char 37\relax}). In particular, the error grows with grid refinement (to $>$8\text{,}\mathrm{\char 37\relax}$$), signifying that the generalization of Peaceman’s model for arbitrarily-oriented wells is not consistent. The result is comparable with the observations in [9, Table 2], where Alvestad’s well indices are compared to a new numerically computed well index, and it is shown that the difference between those two well indices grow, the higher the $r_{\omega}/\Delta x$ ratio. In the variant of the ds model, the error in the source term with respect to the reference solution also grows with larger $r_{\omega}/\Delta x$ ratio. However, the error is consistently smaller (by a factor $>3$) than for the pm well model.

Figure 12 shows the numerical pressure solutions along the $x_{1}$ and the $x_{2}$-axis, for the reference grid resolution. It can be clearly seen that for the reference solution (ds) the pressure solution is regularized. For the variant of ds, the regularization is minimized, however in the far field the solution matches the reference solution better than the pm well model, which is due to the better approximation of the source term (see Fig. 11). We also note that the regularized solution leads to an altered solution in the near-field of the well but to a better approximation of the source term and thus the far field pressure (outside the kernel support), whereas the poor approximation of the source term in the pm method leads to a globally poor pressure solution.

6 Summary

A new well model was presented for which the mass exchange between a well and an embedding porous medium is modeled with a source term spatially distributed by a local kernel function. In the spirit of well-index-based well models the source term for a well with given bottom hole pressure is computed based on the numerical pressure in cells intersecting the well. However, the presented derivation of the new model is independent of the discretization method and the type of computational grid. The new model was shown to be consistent in a numerical experiment and exhibited grid convergence with the expected rates. In the same experiment it is shown that the absolute error with respect to an analytical solution is relatively small, even for coarse computational grids and small kernel support. It was shown, that the error in the source term can be decreased by increasing the region over which the source term is distributed. However, coincidentally, the pressure profile close to the well (inside the kernel support) becomes increasingly regularized. A comparison with a Peaceman-type well model generalized for arbitrarily oriented wells, suggested that even if the region is chosen to be very small (only covering the neighboring cells of cells with well intersection), thus minimizing the regularization effect, the source term can be approximated with good accuracy ($<$ $2\text{\,}\mathrm{\char 37\relax}$ error with respect to a reference solution), whereas the Peaceman-type well model for the same case showed larger differences ($>$ $8\text{\,}\mathrm{\char 37\relax}$) which also had a negative global effect on the pressure solution. The example showed that the new model gives the choice between a more accurate representation of the near-well pressure and a more accurate representation of the source term. Additionally, it was shown that if the source term is accurately approximated, the far-field pressure (outside kernel support) is equally well-approximated while the regularization of the pressure profile only happens locally in the well neighborhood. On the other hand, a bad approximation of the source term leads to global errors in the pressure profile. Finally, the new model was shown to be robust with respect to well rotation, as well as robust with respect to rotations of the anisotropic permeability tensor.

In this work, the well model derivation is restricted to one-phase flow and possible modifications for multi-phase flow are yet to be explored. An extension of the well model for wells with casing is easily conceivable, combining the findings in this work with the derivations presented in [15]. Moreover, the herein presented cases considered the case of a given constant bottom hole pressure. However, the results of [15] indicate, that the presented model may be extended for the coupled 1d-3d case where the well pressure solves an additional one-dimensional partial differential equation in the well domain. Finally, using the superposition principle as shown in [15], which equally applies for the case of anisotropic permeabilities, the presented well model is also expected to provide good approximations when multiple wells are present in the domain.

Acknowledgements

This work was financially supported by the German Research Foundation (DGF), within the Cluster of Excellence in Simulation Technology (EXC 310), and the Collaborative Research Center on Interface-Driven Multi-Field Processes in Porous Media (SFB 1313, Project Number 327154368).

Appendix A Rodrigues’ rotation formula

We want to rotate a given basis $\mathcal{B}=\{\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3}\}$ such that $\boldsymbol{e}_{3}$ is aligned with a vector $\boldsymbol{\psi}$. The Rodrigues’ rotation formula [29] describes a rotated vector $\boldsymbol{x}_{\text{rot}}$ obtained by rotating a vector $\boldsymbol{x}$ by the angle $\theta$ about an axis given by the unit normal vector $\boldsymbol{k}$

[TABLE]

The desired rotation can be described by a rotation by an angle $\theta=\pi$ about the axis $\boldsymbol{k}=\frac{(\boldsymbol{e}_{3}+\boldsymbol{\psi})}{2|\boldsymbol{e}_{3}+\boldsymbol{\psi}|}$,

[TABLE]

Equation A.2 can be expressed in matrix notation as $\boldsymbol{x}_{\text{rot},\pi}=R\boldsymbol{x}$ with

[TABLE]

Appendix B Transformation of Laplace operator

The Joukowsky transformation is a complex function $z=T^{-1}(w)$ that can be decomposed in its real and imaginary parts, $x=\Re{(z)}$ and $y=\Im{(z)}$. Furthermore, let $u=\Re{(w)}$ and $v=\Im{(w)}$. As a conformal mapping, $z$ satisfies the Cauchy Riemann equations [19]

[TABLE]

The Jacobian of the transformation $J_{T^{-1}}$ is given by Eq. 3.16. We investigation the effect of the transformation on the Laplace operator

[TABLE]

where $p$ is analytic in $\Omega_{w}$. Applying the chain rule yields

[TABLE]

Analogously, we arrive at a similar expression for $\partial^{2}p/\partial v^{2}$. Using Eq. B.1 and

[TABLE]

we find that

[TABLE]

With the complex derivative of $T^{-1}$ [21],

[TABLE]

From the determinant of the Jacobian of the transformation, we find

[TABLE]

using Eq. B.1. Hence,

[TABLE]

which also proofs that any harmonic function ($\Delta f=0$) yields another harmonic function after a coordinate transformation with a conformal mapping. \qed

Appendix C Source scaling factor in $w$-coordinates

We want to construct a pressure solution in $w$-coordinates such that the total mass flux over the well boundary matches the specified boundary condition in $x$-coordinates. Hence, the total mass flux over the boundary of a well segment with length $\hat{L}$ in $w$-coordinates needs to match the total mass flux over the boundary of a well segment with length $L$ in $x$-coordinates. A $\hat{q}$ has to be chosen such that $qL=\hat{q}\hat{L}$. A relation between $L$ and $\hat{L}$ can be derived by looking at two related volume integrals. The Joukowsky transformation only affects the two-dimensional well-bore plane such that length $\hat{L}$ of a well segment is not affected. The volume of that well segment in $v$-coordinates (an elliptic cylinder) is given by

[TABLE]

The volume of the same well segment in $x$-coordinates (a circular cylinder with slanted parallel planar elliptic caps) is given by

[TABLE]

where $|E_{\omega,x}|$ is the area of the well-bore ellipse described by Eq. 3.6 transformed to $x$-coordinates (as shown in Fig. 5) and the last equality uses the fact that the integral can be transformed to an integral over a regular cylinder with radius $r_{\omega}$ and length $L$. From the transformation theorem, we know that $V_{x}=V_{v}\operatorname{det}(\tilde{S})$. As the parameter $k_{I}$ is chosen such that $\operatorname{det}(\tilde{S})=1$,

[TABLE]

and if the source term is chosen as

[TABLE]

then $qL=\hat{q}\hat{L}$. \qed

Appendix D Determinant of the Joukowsky transformation

In Appendix B, we show that $\left|\operatorname{det}(J_{T^{-1}})\right|=\left|\frac{\partial z}{\partial w}\right|^{2}$. Using complex differentiation,

[TABLE]

where we used the identities $w^{-1}=\overline{w}|w|^{-2}$, $\overline{w}$ denoting the complex conjugate of $w$, and $\overline{w^{2}}=\overline{w}^{2}$. Furthermore,

[TABLE]

such that

[TABLE]

\qed

Appendix E Extension of Peaceman well model for slanted wells

For a given well direction $\boldsymbol{\psi}=[\psi_{1},\psi_{2},\psi_{3}]^{T}$ and a cell $K_{\Omega}$ with dimensions $\Delta x\times\Delta y\times\Delta z$, to obtain the generalized well model due to [8], reformulated for cell-centered finite volume schemes in [9], replace $k=\sqrt{K_{11}K_{22}}$ in Eq. 5.5 by

[TABLE]

and the expression for $r_{0}$ by

[TABLE]

We note that the formula reduces to Eq. 5.5, if $\boldsymbol{\psi}$ is aligned with a coordinate axis.