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

TL;DR

This paper introduces a method to remove singularities in coupled 1D-3D reservoir flow models, enabling more accurate numerical simulations by reformulating the equations with smooth variables.

Contribution

The authors develop a singularity removal technique for coupled flow models, allowing the use of standard numerical methods on smooth reformulated equations.

Findings

Singularity removal leads to smooth reformulated equations.

The method improves numerical stability and accuracy.

Reformulation resembles a continuous Peaceman well correction.

Abstract

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 sources. 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 singularity removal method for the model equations, resulting in a reformulated coupled 1D-3D flow model in which all variables are smooth. The singularity removal is based on a solution splitting of the reservoir pressure, where it is decomposed into two terms: an explicitly given, lower regularity term capturing the solution singularity and some smooth background pressure. The singularities can then be removed from the system by subtracting them from the governing equations. Finally, the coupled 1D-3D flow equations can be reformulated so they are given in terms of the well pressure and the background reservoir pressure. As these variables are both smooth (i.e. non-singular), the reformulated model has the advantage that it can be approximated using any standard numerical method. The reformulation itself resembles a Peaceman well correction performed at the continuous level.

Figures

Tables

Table 1. Table 1 : The reservoir pressure approximation error p e = ℐ h ​ ( p a ) − p h subscript 𝑝 𝑒 subscript ℐ ℎ subscript 𝑝 𝑎 subscript 𝑝 ℎ p_{e}=\mathcal{I}_{h}(p_{a})-p_{h} and well pressure approximation error p ^ e = p ^ a − p ^ h subscript ^ 𝑝 𝑒 subscript ^ 𝑝 𝑎 subscript ^ 𝑝 ℎ \hat{p}_{e}=\hat{p}_{a}-\hat{p}_{h} when p h subscript 𝑝 ℎ p_{h} was reconstructed using ( 82 ) with k = 1 𝑘 1 k=1 . Both errors were found to converge with optimal order, i.e., with l = 2 𝑙 2 l=2 in the L 2 superscript 𝐿 2 L^{2} -norm and l = 1 𝑙 1 l=1 in the H 1 superscript 𝐻 1 H^{1} -norm.

$h$	${\Vert p_e\Vert }_{L^2\left(\Omega \right)}$	${\Vert p_e\Vert }_{H^1\left(\Omega \right)}$	${\Vert {\hat{p}}_e\Vert }_{L^2\left(\Lambda \right)}$	${\Vert {\hat{p}}_e\Vert }_{H^1\left(\Lambda \right)}$
$1/4$	1.94e-02	1.88e-01	2.30e-3	2.51e-2
$1/8$	5.44-03	4.99e-02	6.27e-4	1.26e-2
$1/16$	1.25e-03	9.26e-02	1.55e-4	6.27e-2
$1/32$	2.77e-04	4.50e-02	7.80e-5	3.32e-2
$l$	2.0	1.0	2.0	1.0

Equations

$$\mathbf{q}+\frac{\kappa }{\mu }\nabla p=$$

$$\nabla \cdot \mathbf{q}=$$

$$\hat{q}+\frac{\hat{\kappa }}{\mu }\frac{\mathrm{d}\hat{p}}{\mathrm{d}s}=$$

$$\frac{\mathrm{d}\hat{q}}{\mathrm{d}s}=$$

$$\bar{p}(z,R)=\frac{1}{2\pi R}{\int }_0^{2\pi }p(R,z,\theta )\mathrm{d}\theta ,$$

$$q=\beta (\hat{p}-\bar{p}),$$

$${\int }_{\Omega }f{\delta }_{\Lambda }\varphi \mathrm{d}\Omega$$

$$q=\beta (\hat{p}-p_K),$$

$$p=\sum _{w=1}^{\text{wells}}E\left(\beta (\hat{p}-\bar{p})\right){\Psi }_w{G}_w+v,$$

$$H^k(\Omega )=\{u\in L^2(\Omega ):D^{\beta }u\in L^2(\Omega )\text{ for }|\beta |\le k\},$$

$$(u,v{)}_{H^k(\Omega )}=\sum _{|\beta |\le k}{\int }_{\Omega }{D}^{\beta }u{D}^{\beta }v\mathrm{d}\Omega .$$

$$H_0^k(\Omega )=\{u\in H^k(\Omega ):u{|}_{\partial \Omega }=0\}.$$

$${\int }_{\Omega }{u}^2{r}^{2\alpha }\mathrm{d}\Omega <\infty ,$$

$$(u,v{)}_{L_{\alpha }^2(\Omega )}={\int }_{\Omega }{r}^{2\alpha }uv\mathrm{d}\Omega .$$

$$H_{\alpha }^1(\Omega )=\{u\in L_{\alpha }^2(\Omega ):D^{\beta }u\in L_{\alpha }^2(\Omega )\text{ for }|\beta |\le k\},$$

$$\mathbf{q}=-\frac{\kappa }{\mu }\nabla p,$$

$${\hat{\mathbf{q}}}_w$$

$$\frac{\mathrm{d}{\hat{\mathbf{q}}}_w}{\mathrm{d}{s}_w}$$

$$q=2\pi \lambda Rf(\hat{p},\bar{p})\text{where }f(\hat{p},\bar{p})=\hat{p}-\bar{p}.$$

$$\mathbf{q}+\frac{\kappa }{\mu }\nabla p=$$

$$\nabla \cdot \mathbf{q}=$$

$$p=$$

$$\hat{q}+\frac{\hat{\kappa }}{\mu }\frac{\mathrm{d}\hat{p}}{\mathrm{d}s}=$$

$$\frac{\mathrm{d}}{\mathrm{d}s}\hat{q}=$$

$$\hat{p}=$$

$$\nabla \cdot \left(-\frac{\kappa }{\mu }\nabla p\right)=$$

$$p=$$

$$\frac{\mathrm{d}}{\mathrm{d}s}\left(-\frac{\hat{\kappa }}{\mu }\frac{\mathrm{d}}{\mathrm{d}s}\hat{p}\right)=$$

$$\hat{p}=$$

$$-\Delta p=f{\delta }_{\Lambda }$$

$$p=f(z)\Psi (r)G(r)+v(r,z).$$

$$G(r)=-\frac{1}{2\pi }\ln (r),$$

$$\Psi (r)=$$

$$\Psi (r)\in$$

$$\Psi (r)=$$

$$-\Delta v$$

$$v$$

$$F=f^{\prime \prime }(z)G(r).$$

Full text

The line source induces a logarithmic type singularity in the reservoir pressure that is difficult to resolve numerically. We present here a singularity removal method for the model equations, resulting in a reformulated coupled 1D-3D flow model in which all variables are smooth. The singularity removal is based on a solution splitting of the reservoir pressure, where it is decomposed into two terms: an explicitly given, lower regularity term capturing the solution singularity and some smooth background pressure. The singularities can then be removed from the system by subtracting them from the governing equations. Finally, the coupled 1D-3D flow equations can be reformulated so they are given in terms of the well pressure and the background reservoir pressure. As these variables are both smooth (i.e. non-singular), the reformulated model has the advantage that it can be approximated using any standard numerical method. The reformulation itself resembles a Peaceman well correction performed at the continuous level.

1 Introduction

Accurate well models are of critical importance for reservoir simulations. The well constitutes the driving force for reservoir flow, in addition to being the main access point of information about its state. The major challenge of well modelling is that of scale disparity; a well has a radius of $\sim 10$ cm, while the reservoir might extend several kilometres in the horizontal plane. From a computational viewpoint, this makes it exceedingly expensive to resolve the well as a 3D object in the grid representing the reservoir. For this reason, wells are typically modelled using either zero-dimensional (0D) point sources or (1D) line sources.

In this work, we take as a starting point the coupled 1D-3D flow model

[TABLE]

where $\Omega\subset\mathbb{R}^{3}$ denotes the reservoir domain and $\Lambda=\cup_{w=1}^{\text{wells}}\Lambda_{w}\subset\mathbb{R}^{1}$ a collection of line segments each representing a well. The 1D domain is parametrized by its arc-length $s$. The parameters $\kappa$, $\hat{\kappa}$ and $\mu$ denote reservoir permeability, well permeability and fluid viscosity, respectively, and are assumed to be positive and constant. The variables $p$ and $\mathbf{q}$ denote fluid pressure and flux in the reservoir, $\hat{p}$ and $\hat{q}$ fluid pressure and flux in the well, and $\bar{p}$ the reservoir pressure averaged over the surface of the borehole $r=R$:

[TABLE]

as is illustrated in Figure 1.

Physically, equations (1a)-(1b) describe a Darcy-type flow in the reservoir domain $\Omega$, and equations (1c)-(1d) a Poiseuille-type flow in the well. The latter is a 1D flow equation, where the radial and angular components have been neglected. For a description of this model reduction method for the well flow, we refer to the work of Cerroni et al. in Cerroni2019 . The mass flux $q$ between reservoir and well is modelled using a linear filtration law,

[TABLE]

which states that the connection flow between them is proportional to their pressure difference. The proportionality coefficients $\beta,\hat{\beta}\in C^{1}(\Lambda_{w})$ are assumed piecewise continuous and allowed to vary along the well. The wells are considered as concentrated line sources $\delta_{\Lambda}$ in the reservoir equation (1a)-(1b), with the line sources defined in the following manner:

[TABLE]

for all $\phi\in C^{0}(\Omega)$, with $s_{w}$ denoting the arc-length of line segment $\Lambda_{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 llau2016 or the interference of metallic pipelines and bore-casings in electromagnetic modelling of reservoirs weiss2017 . A coupled 1D-3D heat transfer problem was considered in the context of geothermal energy in alkhoury2005 , 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 nabil2016 , the efficiency of drug delivery through microcirculation Cattaneo2014 ; zunino2018 , and the study of blood flow in the vascularized tissue of the brain Reichold2009 ; Grinberg2011 . 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

[TABLE]

where $p_{K}$ denotes the average pressure in the grid block containing the well. Due to the singularity, $p_{K}$ will not be representative of the reservoir pressure at the bore-hole; this is typically accounted for by multiplying $\beta$ with a well index $J$. A correction of this type was first developed by Peaceman in peaceman1978 , 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 grid-blocks and anisotropic permeability Peaceman-rectangle . The extension to more generalized grids was treated by e.g. Aaavatsmark in Aavatsmark2016-0 ; Aavatsmark2016 ; Aavatsmark2016-2 , to more generalized flow models by e.g. Ewing in Ewing1999 , and to more generalized discretization schemes by e.g. Chen et al. in Chen2009 . Many authors have contributed to the extension to generalized well placements, we mention here the work of King et. al in king , Aavatsmark in aavatsmark2003index , and of special relevance to our work, that of Wolfsteiner et al. in Wolfsteiner2003 and Babu et al. in babu .

In this work, we take a different approach, in which the singularities are explicitly removed from the governing equations. We start by showing that the reservoir pressure $p$ admits a splitting

[TABLE]

where $G_{w}$ is a given logarithmic function that captures the near-well behaviour of the reservoir pressure, $E$ is an extension operator $E\colon H^{2}(\Lambda)\rightarrow H^{2}(\Omega)$, $\Psi_{w}$ some smooth cut-off function, and $v\in H^{2}(\Omega)$ some higher-regularity remainder term. The key point here is that the singular nature of the solution is explicitly captured by the logarithmic terms $G_{w}$. With the splitting (6) in hand, we can therefore remove the singular terms from the system by straightforward subtraction. Finally, we reformulated coupled 1D-3D flow model can then be reformulated so it is given with respect to the high-regularity variables $\hat{p}$ and $v$. The main contribution of this article is the reformulation of the coupled 1D-3D flow model into equations (34a)-(34d), for which the solution is smooth (non-singular). On a practical level, this means the solution can be approximated using any standard numerical method.

The technique of removing singularities is commonly known for point sources; we refer here to (ewing-book, , 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 hales . A splitting of the type (6) was introduced by Ding in Ding2001 for the point source problem, where it was used to formulate grid refinement strategies. We are, to the best of our knowledge, the first to formulate a singularity removal method for the coupled 1D-3D flow problem. Central to this method is the construction of a function $G_{w}$ capturing the solution singularity; we use here a function $G_{w}$ found by integrating the Green’s function for the reservoir equations (1a)-(1b) over the line $\Lambda$; we refer here to our earlier work in (Gjerde2018, , Section 3.2). This use of Green’s functions to construct analytical and semi-analytical well models has a rich history. Of special relevance to our work, we mention that of Wolfsteiner et al. and Babu et al. in babu ; Wolfsteiner2003 , in which the Green’s function was used to construct analytical solutions with which to calculate the well index $J$. More recently, Nordbotten et al. used Green’s functions to construct analytical models to estimate leakage of CO2 stored in geological formations Nordbotten2009 .

The singularity removal, and subsequent reformulation of the model in terms of the smooth variables $v$ and $\hat{p}$, is similar to the Peaceman well correction in that it leads to an alteration of the inflow parameter $\beta$. We discuss this in more detail in Section 7. It differs, however, in that it works on the continuous level. It is therefore easily adapted to different discretization methods, generalized well placements within the domain and different types of boundary conditions. Moreover, since our method gives an explicit representation of the logarithmic nature of the solution, it allows us to accurately represent the reservoir pressure in the whole domain (including in the near-vicinity of the well).

In our presentation of the method, we limit ourselves to considering a linear reservoir equation with constant, scalar-valued permeabilities and Poiseuille flow in the well. The latter restriction is not critical to the methodology; the well equation could for example be taken non-linear as long as the well pressure remains sufficiently regular. To be more precise, the method requires $\hat{p}$ to be piecewise $C^{1}$ on $\Lambda$. As for the reservoir equation, the reservoir pressure could be replaced with a potential expression $\phi$ so that the effect of gravity can be included. The singularity removal and reformulation can be extended to handle spatially varying, scalar-valued permeabilities as shown in Gjerde2018 . For an extension to tensor-valued permeabilities and non-linear reservoir equations, we suggest using the solution splitting in (6) to formulate a multiscale finite volume method such as in hamdi2006 , or a generalized finite element method gfem , where the analytic functions capturing the solution singularity are used to enrich the set of basis functions.

For the discretization and numerical experiments, we consider herein the Galerkin Finite Element (FE) method. The FE approximation of the line source problem was studied by D’Angelo in dangelo2012 by means of weighted Sobolev spaces, using similar techniques as those known for e.g. corner-point problems babuska1972 . D’Angelo proved that the approximation of the coupled 1D-3D flow problem (1a)-(1b) converges sub-optimally unless the mesh is sufficiently refined around the well. The sub-optimal convergence rates were found to be local to the line source by Köppl et al. in koppl2015 , meaning that they only pollute the pressure approximation inside the well block. However, this means the approximation of the coupled 1D-3D flow problem will suffer until the mesh parameter $h$ is smaller than the well radius $R$. In practice, one therefore needs a very fine mesh around the well for the FE approximation of (1a)-(1d) to converge. This makes the problem computationally expensive to solve. Different strategies have been proposed to remedy this, e.g., Kuchta et al. studied suitable preconditioners in miro2016-3D1D . Holter et al. then applied this preconditioner to simulate flow through the microcirculature found in a mouse brain Holter2018 . An alternative coupling scheme was introduced by Köppl et al. in koppl2016 , 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 for the problem. In Section 3, we introduce in more detail the coupled 1D-3D flow model we take as a starting point. In Section 4, we show that the reservoir pressure $p$ admits a splitting into lower-regularity terms that capture the solution singularities, and a higher-regularity remainder term $v$. With the splitting in hand, the singularities can then be subtracted from the governing equations. The result is the reformulated coupled 1D-3D flow model (34a)-(34d), posed in terms of the smooth variables $\hat{p}$ and $v$. As the solutions then enjoy significantly improved regularity, this system 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 two numerical experiments, 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$.

2 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}(\Omega)$ be the Sobolev space,

[TABLE]

with $\beta$ denoting a multi-index and $D^{\beta}$ the corresponding weak distributional derivative of $u$. $H^{k}(\Omega)$ is a Hilbert space endowed with inner product

[TABLE]

We use a subscript to denote the subspace of $H^{k}$ with zero trace on the boundary, $H^{k}_{0}$, i.e.,

[TABLE]

As we will see, the reservoir solution $p$ in (1a)-(1d) fails to belong to $H^{1}(\Omega)$ due to singular behaviour on $\Lambda$. For this reason, we consider also a weighted Sobolev space. To define it, let $-1<\alpha<1$, and take $L^{2}_{\alpha}(\Omega)$ to denote the weighted Hilbert space consisting of measurable functions $u$ such that

[TABLE]

where $r$ denotes the distance of a point to $\Lambda$, i.e., $r(\mathbf{x})=\text{dist}(\mathbf{x},\Lambda)$. This space is equipped with the inner product

[TABLE]

For $\alpha>0$, the weight $r^{\alpha}$ has the power to dampen out singular behaviour in the function being normed; for $\alpha<0$, the weight function can induce or worsen already singular behaviour. We therefore have the relation $L^{2}_{-\alpha}(\Omega)\subset L^{2}(\Omega)\subset L^{2}_{\alpha}(\Omega)$ for $\alpha>0$. Letting now $H^{1}_{\alpha}(\Omega)$ be the Sobolev space

[TABLE]

we will later find that the reservoir pressure solving (1a)-(1d) belongs to $H^{1}_{\alpha}(\Omega)$ for $\alpha>0$.

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. Consequently, 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 apel2011 ; dangelo2012 ; Ding2001 . 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^{\alpha}$.

3 Mathematical model

Here, we introduce in more detail the coupled 1D-3D equation we take as a starting point. Let $\Omega\subset\mathbb{R}^{3}$ denote a bounded domain describing a reservoir, with smooth boundary $\partial\Omega$. We consider here steady-state, incompressible Darcy flow

[TABLE]

where $\mathbf{q}$ and $p$ denote reservoir flow and pressure, $\mu$ the fluid viscosity, and $\kappa$ a given positive and scalar permeability. We consider also a collection of wells, each considered to be a thin tube with fixed radius $R$ and centreline $\Lambda_{w}$. The centreline is parametrized by the arc length $s_{w}$. We denote by $\boldsymbol{\tau}_{s_{w}}$ its normalized tangent vector. As the radius of the tube is small, we assume the radial and angular components of the well pressure can be neglected, meaning $\hat{p}|_{\Lambda_{w}}=\hat{p}(s_{w})$. The well flow domain $\Lambda$ will then consist of a collection of line segments, $\Lambda=\cup_{w=1}^{\text{wells}}\Lambda_{w}$. We consider on this domain Poiseuille-type flow,

[TABLE]

with $\hat{\mathbf{q}}_{w}$ and $\hat{p}_{w}$ denoting flow and pressure in the well and $q$ the linear mass flux into or out of the well. $\frac{\mathrm{d}}{\mathrm{d}s_{w}}$ denotes the derivative with respect to the tangent line, or equivalently, the projection of $\nabla$ along $\boldsymbol{\tau}$, i.e., $\frac{\mathrm{d}}{\mathrm{d}s_{w}}=\nabla\cdot\boldsymbol{\tau}_{s_{w}}$. As the fluid flux in the well has a fixed direction, it can be given as a scalar function $\hat{q}_{w}$, characterized by the property $\hat{\mathbf{q}}_{w}=\hat{q}_{w}\boldsymbol{\tau}_{w}$. Note that the assumption of Poiseuille flow is not critical; (8a) could for example contain certain non-linearities.

Letting now $\Lambda=\cup_{w=1}^{\text{wells}}\Lambda_{w}$ denote the collection of line segments $\Lambda_{w}$, the well pressure and flux can be written as 1D variables $\hat{p},\hat{q}\colon\Lambda\rightarrow\mathbb{R}$. The well and reservoir flow can then be coupled together using a linear filtration law, which states that the mass flux $q$ between them is proportional to their pressure difference:

[TABLE]

The mass flux is given as the rate of transfer per unit length, and the variable $\lambda\in C^{2}(\Lambda)$ denotes the permeability of the borehole lateral surface. It accounts for the fact that the well may not be in perfect contact with the reservoir, leading to a pressure drop across the borehole. Letting $\Delta p_{\text{skin}}$ denote this pressure drop, this can be expressed by the following relation: $q=2\pi R\lambda\Delta p_{\text{skin}}$.

The pressure difference $f(\hat{p},\bar{p})$ between well and reservoir uses an averaged value $\bar{p}(z;R)$ for the reservoir pressure given in (2). 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

[TABLE]

where $\hat{\kappa}=R^{2}/{8}$, $\beta=2\pi R\lambda$, $\hat{\beta}=\beta/\pi R^{2}$. The functions $p_{D}\in C^{2}(\bar{\Omega})$ and $\hat{p}_{D}(\bar{\Lambda})$ denote given boundary data. The connection flow from well to reservoir is modelled by means of a generalized Dirac delta function $\delta_{\Lambda}$, 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:

[TABLE]

with $f(\hat{p},\bar{p})=\hat{p}-\bar{p}$.

4 Splitting Properties of the Solution

In this section, we will show that the line source in the right-hand 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 discuss in detail the splitting when $\Lambda$ is assumed a single line segment aligned with the $z$-axis, $\frac{\kappa}{\mu}=1$ and the well outflow $q$ is a given function $f\in C^{1}_{0}(\Lambda)$. The splitting is then especially simple; this case therefore serves to illustrate the splitting method itself. We then generalize it in two steps, handling first an arbitrary line segment and $\frac{\kappa}{\mu}\neq 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 (34a)-(34d), wherein the singularity has been removed and all variables are smooth.

4.1 Elliptic equations with a single line source

In this section, we consider the elliptic equation

[TABLE]

when $\Lambda$ and $\Omega$ are as illustrated in Figure 1, and $f=f(z)\in C^{1}_{0}(\Lambda)$ is a given, smooth line source intensity (assumed zero at the endpoints of $\Lambda$). The solution $p$ then admits a splitting into an explicit, low-regularity term $f(z)\Psi(r)G(r)$, and an implicit, high-regularity term $v$:

[TABLE]

Here, $G(r)$ captures the singular part of the solution, and is given by

[TABLE]

and $\Psi(r)$ denotes some smooth cut-off function satisfying

[TABLE]

Assuming the cut-off radius $R_{c}$ is chosen small enough to satisfy $\Psi(r)=0$ on $\partial\Omega$, the regular component $v$ can then be defined as the solution of

[TABLE]

where

[TABLE]

To see that $p$ given by (13) indeed solves (12), let us first note that $G=-\nicefrac{{1}}{{2\pi}}\ln(r)$ was so chosen because it satisfies $-\Delta G=\delta_{\Lambda}$. To be more precise, $G$ is the fundamental solution of the Laplace equation in 2D, and thus has the property

[TABLE]

Considering then the Laplacian of $p$ given by (13), a straightforward calculation shows that all but one term vanish by construction, i.e.,

[TABLE]

By (18), we then find that

[TABLE]

and it follows that the $p$ constructed in (13) indeed solves (12) in a suitably weak sense.

Formally speaking, the splitting works by introducing first the logarithmic term $G$ for which the Laplacian returns the line source with the required intensity $f$. The higher-regularity term $v$ is then used to correct the solution so it solves the original problem. The existence of such a function $v$ follows from standard elliptic theory. As $\ln(r)\in L^{2}(\Omega)$, and $f^{\prime\prime}(z)\in L^{2}(\Lambda)$ by assumption, one can show that the entire right-hand side $F$ in (16a) belongs to $L^{2}(\Omega)$ (Gjerde2018, , Section 3.1). Consequently, there exists $v\in H^{2}(\Omega)$ solving (16a)-(16b). The full solution $p$, meanwhile, fails to belong to $H^{1}(\Omega)$. This can be shown by straightforward calculation, as one has $\ln(r)\in L^{2}(\Omega)$ but $\nabla\ln(r)\,\cancel{\in}\,L^{2}(\Omega)$. Instead, one has $p$ belonging to the weighted Sobolev space $H^{1}_{\alpha}(\Omega)$ for any $\alpha>0$. It follows that $v$ is indeed the higher-regularity term in the splitting (13). Formally, this means 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.

4.2 Elliptic equations with an arbitrary line source

In this section, we consider the elliptic problem

[TABLE]

when the right-hand side is a line source $\delta_{\Lambda}$ located on a single line segment $\Lambda$ with endpoints $\mathbf{a},\mathbf{b}\in\Omega$. The line $\Lambda$ can be described by the parametrization $\mathbf{y}=\mathbf{a}+\boldsymbol{\tau}s\quad\text{for }s\in(0,L),$ where $L=\|\mathbf{b}-\mathbf{a}\|$ denotes the Euclidean norm and $\boldsymbol{\tau}=(\mathbf{b}-\mathbf{a})/L$ is the normalized tangent vector of $\Lambda$. Letting again $f=f(s)\in C^{1}(\Lambda)$ be a given 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$:

[TABLE]

The function $G$ is now given by

[TABLE]

with $r_{b}(\mathbf{x})=\|\mathbf{x}-\mathbf{b}\|$ and $r_{a}(\mathbf{x})=\|\mathbf{x}-\mathbf{a}\|$. This function was constructed by integrating the 3D Green’s function for (11a) (when posed in $\mathbb{R}^{2}$) over the line segment $\Lambda$. It thus satisfies the property $\nabla\cdot(-\frac{\kappa}{\mu}\nabla G)=\delta_{\Lambda}$ (Gjerde2018, , Section 3.2). Next, $E$ denotes an extension operator $E\colon H^{2}(\Lambda)\rightarrow H^{2}(\Omega)$ extending $f$ so that it can be evaluated in the entire domain $\Omega$. Assuming again that the cut-off function $\Psi$ satisfies $\Psi=0$ on $\partial\Omega$, the regular component $v$ is then defined as the solution of

[TABLE]

where

[TABLE]

To see that the constructed $p$ indeed solves the right problem, let us start by inserting it into (20). construction, all terms disappear except $E(f)\Psi\Delta G$. Integrating this term over the domain, we find that

[TABLE]

for all $\forall\phi\in C^{0}(\Omega)$, where we used the property that $E(f)=f$ on $\Lambda$. It follows that the $p$ constructed in (13) indeed solves (12) in a suitably weak sense.

By a similar argument as the one given in (Gjerde2018, , Section 3.2), one finds that $F$ given by (24) belongs to $L^{2-\epsilon}(\Omega)$ for arbitrarily small $\epsilon>0$. It follows that there exists $v\in H^{2-\epsilon}(\Omega)$ solving (23a)-(23b). Moreover, a straightforward calculation shows that $G$ again fails to belong to $H^{1}(\Omega)$. In fact, one has $G\in H^{1-\epsilon}(\Omega)$. It follows that $v$ constitutes the higher-regularity component of the solution split (21), meaning that $v$ is smoother and better behaved than the full solution $p$.

4.3 The coupled 1D-3D flow problem

Let us now consider the coupled 1D-3D flow problem (11a)-(11d). To start with, let us again consider a single line segment $\Lambda$ with endpoints $\mathbf{a},\mathbf{b}\in\Omega$. From the discussion in the preceding section, it is natural to assume $p$ solving (11a)-(11d) admits a solution splitting of the type:

[TABLE]

with $G$ being as in (22), $\Psi$ being some smooth cut-off function, $f$ being the previously introduced pressure difference $f=\hat{p}-\bar{p}$, and $v$ defined as the solution of

[TABLE]

with

[TABLE]

Unlike in Sections 4.1 and 4.2, $f=f(\hat{p},\bar{p})$ is now implicitly given from $\hat{p}$ and $\bar{p}$ solving the coupled 1D-3D flow problem. To reformulate (11a)-(11d) in terms of $\hat{p}$ and $v$, the right-hand side therefore needs to be reformulated. To this end, let us first treat the pressure difference $\hat{p}-\bar{p}$. By the splitting (26) and the definition of the averaging in (2), calculations reveal that

[TABLE]

Here we used the simplifications $\overline{E(f)}=f|_{\Lambda}$ and $\overline{\Psi}\approx 1|_{\Lambda}$. This is motivated by the fact that the well radius $R$ is assumed negligible. From this, we can state the reformulated coupled 1D-3D flow model:

[TABLE]

where

[TABLE]

$\beta^{*}$ is given by

[TABLE]

and $\hat{\beta}^{*}=\beta^{*}/\pi R^{2}$.

The extension to multiple wells follows naturally by applying the superposition principle. Considering now $\Lambda=\cup_{w=1}^{\text{wells}}\Lambda_{w}$, with each line segment $\Lambda_{w}$ having endpoints $(\mathbf{a}_{w},\mathbf{b}_{w})\in\Omega$, we can formulate a solution splitting

[TABLE]

where $E\colon H^{2}(\Lambda)\rightarrow H^{2}(\Omega)$ is the same extension operator as before,$G_{w}$ is given by (22) with $\mathbf{a}=\mathbf{a}_{w}$ and $\mathbf{b}=\mathbf{b}_{w}$, $\Psi_{w}$ is some smooth cut-off function with respect to line segment $\Lambda_{w}$, and $v$ solves

[TABLE]

with right-hand side

[TABLE]

and

[TABLE]

The system (34a)-(34d) constitutes a reformulation of the coupled 1D-3D flow model in terms of the smooth variables $v$ and $\hat{p}$. For an example of what the splitting might look like, the reader is invited to examine Figure 2. As the singularities have here been removed from the system, it enjoys significantly improved regularity compared to the standard formulation (11a)-(11d).

5 Weak formulation

In this section, we state a weak formulation of the reformulated coupled 1D-3D flow problem (34a)-(34d). 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 $\mathbf{V}$ denote the product space $\mathbf{V}=V\times\hat{V}$, where

[TABLE]

normed by

[TABLE]

Multiplying (34a) and (34c) with test functions $\phi\in H^{1}_{0}(\Omega)$ and $\hat{\phi}\in H^{1}_{0}(\Lambda)$, respectively, integrating over their respective domains, and performing an integration by parts, we arrive at the following variational formulation:

Find $(v,\hat{p})\in\mathbf{V}$ such that

[TABLE]

for all $(\phi,\hat{\phi})\in\mathbf{V}_{0}$, where

[TABLE]

and

[TABLE]

The full reservoir pressure can then be constructed from $v$ and $\hat{p}$ by the relation

[TABLE]

Next, let us consider the standard coupled 1D-3D flow model, and give its variational formulation as it was proposed in dangelo2008 . Let $\mathbf{V}_{\alpha}$ denote the weighted product space $\mathbf{V}_{\alpha}=V_{\alpha}\times\hat{V}$, where

[TABLE]

normed by

[TABLE]

Multiplying (11a) and (11c) with test functions $v\in H^{1}_{-\alpha,0}(\Omega)$ and $\hat{v}\in H^{1}_{0}(\Lambda)$, respectively, integrating over their domain of support, and performing an integration by parts, we arrive at the variational formulation:

Find $(p,\hat{p})\in\mathbf{V}_{\alpha}$ such that

[TABLE]

for all $(\phi,\hat{\phi})\in\mathbf{V}_{-\alpha,0}$, where

[TABLE]

and the test space $\mathbf{V}_{-\alpha,0}$ is the space of functions $(\phi,\hat{\phi})\in\mathbf{V}_{-\alpha,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 (48). For a proof of the well-posedness of this formulation, the reader is referred to dangelo2008 ; dangelo2012 .

6 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(\hat{p},\bar{v})=\hat{p}-\bar{v}$ now uses the regular part of the pressure, $v\in H^{2}(\Omega)$, we introduce here also the simplification $\bar{v}_{h}=v_{h}|_{\Lambda}$; i.e., we take the trace of $v_{h}$ on $\Lambda$ 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, meaning $\bar{v}\approx v|_{\Lambda}$. 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 $\Lambda$.

We will now give the discretized form of the variational formulation (40). For simplicity, let us assume $\Omega$ is a polyhedron that readily admits a partitioning $\mathcal{T}_{T,h}$ into simplicial elements $T$:

[TABLE]

The simplicial partitioning $\mathcal{T}_{T,h}$ forms a mesh, assumed conforming, which can then be characterized by the mesh size $h=\max_{T\in\mathcal{T}_{T,h}}h_{T}$. Next, we associate this mesh with the usual (3D) Lagrange space of order 1, $V^{h}_{u}$, given by

[TABLE]

Here, $\mathbb{P}^{1}$ denotes the space of polynomials of degree $1$, and $C^{0}_{u}(\Omega)$ the space of continuous elements that equal the interpolation of $u$ on the boundary, i.e.,

[TABLE]

Next, we assume $\Lambda$ admits a partitioning $\mathcal{T}_{I,h}$ into line segments $I$:

[TABLE]

assumed again to satisfy all the requirements of a conforming mesh, and associated with the mesh size $\hat{h}=\max_{I\in\mathcal{T}_{I,h}}h_{I}$. For the discretization of $\hat{V}$, we use the (1D) Lagrange space of order 1,

[TABLE]

with $C^{0}_{\hat{u}}(\Lambda)$ interpreted as in (49).

Considering first the reformulated system (40), let

[TABLE]

where $\{\phi_{1},\phi_{2},...,\phi_{N}\}$ and $\{\hat{\phi}_{1},\hat{\phi}_{2},...,\hat{\phi}_{\hat{N}}\}$ are linear hat functions spanning $V^{h}$ and $\hat{V}^{h}$, respectively. Note next that $v_{h}$ is a linear function used to approximate the high regularity term $v\in H^{2}(\Omega)$. For $R\ll h$, its average $\bar{v}_{h}$ can be well approximated by simply taking the trace $v|_{\Lambda}$. The pressure difference $\hat{p}-\bar{v}$ is then given by

[TABLE]

Here, $T\colon V^{h}\rightarrow\hat{V}^{h}$ is the discrete trace matrix, characterized by the property $\phi_{k}|_{\Lambda}=\sum_{l=1}^{\hat{N}}T_{k,l}\hat{\phi}_{l}$.

Testing (40) with $v=\phi_{i}$ for $i=1,...,N$ and $\hat{v}=\hat{\phi}_{j}$ for $j=1,...,\hat{N}$, we arrive at the following discrete system:

[TABLE]

where $A$ and $\hat{A}$ are the standard stiffness matrices

[TABLE]

$\hat{M}$ denotes the standard 1D mass matrix,

[TABLE]

and $C$ denotes the coupling block,

[TABLE]

We will refer to this system as the Singularity Removal Based FE method. After solving (52), a discretization of the full reservoir pressure $p_{h}$ can be reconstructed using

[TABLE]

where $\mathcal{I}_{h}^{k}$ denotes the interpolation onto the Lagrange space of order $k$. As the interpolation of $G(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 (47) directly; this is the finite element formulation analysed in e.g. dangelo2012 . 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

[TABLE]

The pressure difference $\hat{p}-\bar{p}$ is then given by

[TABLE]

where $\Pi$ is the discrete averaging matrix $\Pi:V^{h}\rightarrow X^{h}$ and $\{\hat{\psi}_{1},\hat{\psi}_{2},...,\hat{\psi}_{\hat{M}}\}$ are the basis functions spanning $\hat{X}^{h}$.

Testing now (47) with $v=\phi_{i}$ for $i=1,...,N$ and $\hat{v}=\hat{\phi}_{j}$ for $j=1,...,\hat{N}$, we arrive at the following block system for the discretization of (47):

[TABLE]

Here, $N$ denotes the mass matrix given by

[TABLE]

for $\hat{\psi}_{m}$ belonging to the discontinuous Galerkin space of order 0:

[TABLE]

We will refer to this system as the standard FE method.

7 Relation to the Peaceman well model

In this section, we show that the reformulated coupled 1D-3D flow model (34a)-(34d) 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 peaceman1978 . 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):

[TABLE]

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\ll 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

[TABLE]

Integrating this equation to a radius $r_{e}$,

[TABLE]

we find that

[TABLE]

when $p_{e}=p(r_{e})$. We also need to take into account the pressure drop $\Delta p_{\text{skin}}$ across the skin of the well. To do so, let $S$ be the skin-factor, defined by the relation

[TABLE]

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}$ peaceman1978 :

[TABLE]

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 peaceman1978 .

The reformulation of the pressure difference $f$ in terms of $\hat{p}$ and $v$ bears a strong resemblance to the Peaceman well correction in (67). In a practical sense, the reformulation into (34a)-(34d) 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

[TABLE]

Next, we let now $\hat{p}$ 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

[TABLE]

Note that this $G$ is not smooth enough to work for the solution split (33), we use it here only for the sake of comparison. By the definition of the averaging (2), we have $\overline{G}=-\nicefrac{{\mu}}{{2\pi\kappa}}\ln(\nicefrac{{R}}{{r_{e}}})$ Inserting it in (68) yields the relation

[TABLE]

Here, $\nicefrac{{2\pi\kappa}}{{\mu\beta}}$ can be substituted by the skin factor of the well by recalling $q=\beta\Delta p_{\text{skin}}$. This results in an expression that equals the Peaceman well correction given in (67), i.e.,

[TABLE]

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 $\beta$ (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., $\hat{p}-\bar{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$.

8 Numerical Results

In this section, we perform numerical experiments to test the approximation properties of the Singularity Subtraction Based FE method given by (52). For the implementation, we utilized the finite element framework FEniCS LoggMardalEtAl2012a . For the first test case, we consider a single well with smooth lateral well permeability $\beta$, and compare the results against those obtained using the standard FE given by (60). Our implementation of this method uses an earlier implementation from Kuchta miro-git , the same as was utilized for the results of Holter et al. in Holter2018 . 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 cecile-git . For the second test case, we consider a discontinuous lateral permeability $\beta$, and an extension operator that uses radial basis function interpolation. We show here that the reconstructed reservoir pressure $p_{h}$ converges optimally when the Singularity Removal Based FE method is applied.

8.1 Convergence test for well with smooth lateral permeability

In this section, we take $\Omega=(0,1)^{3}$ and $\Lambda=\{(x,y,z)\in\Omega:x=y=\nicefrac{{1}}{{2}}\}$. We want to test the capability of each method in approximating the test problem

[TABLE]

with the following parameters:

[TABLE]

The solution, along with the splitting terms, are shown in Figure 2.

In order to test the stability of the approximation when the well radius is small compared to mesh size $h$, we test using four different values for the well radius:

[TABLE]

Furthermore, we set $\Psi=1$ and choose as the extension operator

[TABLE]

In this case, the reformulated FE method will approximate the analytic solution for $v_{a}$ given in (73b), meaning we can compute its error directly using $\|v_{a}-v_{h}\|$.

Figure 3 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 3(a) show the approximation error of $v_{h}$, measured in the $L^{2}$-norm, i.e., $\|v_{h}-v_{a}\|_{L^{2}(\Omega)}$ with $v_{a}$ being the analytic solution in (73b). 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 (52) the simplification $\bar{v}=v|_{\Lambda}$, and this is not justified for $R\sim h$. Optimal convergence rates could be restored by taking the average of $v_{h}$ rather than its trace.

The red lines in Figure 3(a) give the approximation errors for the full reservoir pressure using the standard FE method described by (60). We give here the approximation error of $p_{h}$ in the $L^{2}$-norm, i.e., $\|p_{h}-p_{a}\|_{L^{2}(\Omega)}$. 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\sim 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 dangelo2012 that this leads to a reduction in the convergence rate of $p_{h}$. We refer here to our comments in (Gjerde2018, , p. 14-15) for a more in-depth explanation of this, and remark only that the line source problem is expected to converge with order $h^{1-\epsilon}$ for $\epsilon>0$ arbitrarily small. Thus, the convergence order $h^{l}$ with $l=1.4$ surpasses the theoretical expectation when $R\sim h$.

The blue and red lines in Figure 3(b) give the approximation error of $\hat{p}_{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 $\|\hat{p}_{h}-\hat{p}_{a}\|_{L^{2}(\Lambda)}$. 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 simplification $\bar{v}=v|_{\Lambda}$ used in the construction of the block matrix (52), and could be resolved by removing this simplification.

From 3(b), is clear that the standard FE method has trouble approximating the solution when $R<h$. Moreover, the approximation error of $\hat{p}_{h}$ is seemingly more sensitive than $p_{h}$ with respect to the ratio $R/h$. This can be understood by returning to the reservoir pressure splitting $p=\beta(\hat{p}-\bar{p})G+v$, where $G=-\nicefrac{{1}}{{2\pi}}\ln(r)$, 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., $\|\hat{p}_{h}-\bar{p}_{h}-(\hat{p}_{a}-\bar{p}_{a})\|_{L^{2}(\Lambda)}$, the error in approximating the logarithm, i.e., $\|\ln(r)_{h}-\ln(r)\|_{L^{2}(\Omega)}$, and the error in approximating $v$ (which is comparatively small). The standard FE method has trouble resolving the logarithmic nature of the reservoir pressure around the well, leading to a large approximation error in $\bar{p}$. This further pollutes the approximations of both $\hat{p}$ 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 pressure $\hat{p}$, 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 $\|\bar{p}_{a}-\bar{p}_{h}\|_{L^{2}(\Lambda)}$.

In summary, we see here that the standard FE method has difficulty resolving the pressure difference $\hat{p}-\bar{p}$ when $R<h$, due to the fact that $\bar{p}$ is then poorly approximated. This further pollutes 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 (52), restores optimal convergence rates for the reservoir pressure $p$, and improves the robustness of the method with respect to a small well radius $R$.

8.2 Convergence test for well with discontinuous lateral permeability

Let $\Omega=(0,1)^{3}$ and

[TABLE]

In this section, we will test the ability of the Singularity Remov1l Based FE method in approximating the analytic test problem

[TABLE]

when $G$ is given as in Section 4.2:

[TABLE]

The problem parameters are then as follows:

[TABLE]

Physically, this can be interpreted as modelling a well that passes through the domain but is only in contact with the reservoir when $\nicefrac{{1}}{{4}}<z<\nicefrac{{3}}{{4}}$. This translates to a jump in the lateral permeability, with discontinuities at the points $(\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}},\nicefrac{{1}}{{4}})$ and $(\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}},\nicefrac{{3}}{{4}})$.

As the cut-off function, we use the Gaussian function

[TABLE]

with $c=0.04$. For the extension operator $E$, we choose spline interpolation with radial basis functions as given in splineinterpolation . Given a discretized solution pair $(v_{h},\hat{p}_{h})$ to (52), we can then reconstruct the discretized full reservoir pressure by the relation

[TABLE]

where $\mathcal{I}_{h}^{k}$ denotes the interpolation operator onto the Lagrange elements of order $k$. Finally, the numerical error associated with $v_{h}$ can be computed as

[TABLE]

where $p_{a}$ is interpolated onto the Lagrange elements with the same order as the solution $v_{h}$.

The results of applying the SRB-FE method to solve this problem are plotted in Figure 4 for $h=\nicefrac{{1}}{{8}}$. The errors and convergence rates are given for different mesh sizes in Table 1. As is evident from this table, the SRB-FE approximation of $p_{h}$ and $\hat{p}_{h}$ both converge with optimal order. I.e., we find that

[TABLE]

9 Conclusion

In this work, we have developed a singularity removal method for the coupled 1D-3D flow model. This type of model can be used to model the interaction of wells with a reservoir. The well is endowed with its own 1D flow equation, and modelled as a 1D line source in the reservoir domain. This line source introduces a logarithmic type singularity in the reservoir solution that negatively affects the approximation properties of the problem. We provide here a method for identifying and removing this singularity from the governing equations. The result is a reformulated coupled 1D-3D flow model in which all variables are smooth.

As the reformulated model is posed in terms of smooth variables, it has the advantage that it can be approximated using any standard numerical method. In this work, we have shown that the singularity removal restores optimal convergence rates for the Galerkin FE method. Moreover, it makes the approximation stable with respect to the ratio $R/h$ between well radius and mesh size.

A natural development of this work consists of extending the singularity removal method to apply to (i) different control modes for the wells, (ii) tensor-valued permeability and (iii) a mixed formulation of the flow, where both pressure and flux are approximated. We believe these extensions would be particularly valuable in the context of subsurface flow applications, as it would allow one to capture the interaction between well and reservoir using coarse grids. The extension to different control modes for the wells, i.e., rate controlled or pressure controlled wells, is straightforward; it can be achieved by altering the boundary conditions for the well flow equations. As the singularity subtraction is performed at the continuous level, it is likewise straightforward to adapt the method to different discretization methods ecmor . The extension to tensor-valued permeability is more challenging, and will be treated in future work.

Acknowledgements.

This work was partially supported by the Research Council of Norway, project number 250223. The participation of the first author at ECMOR XVI 2018 was funded by the Academia Agreement between the University of Bergen and Equinor. The authors thank B. Wohlmuth for her help in formulating the original splitting approach. The authors also thank P. Zunino for the interesting discussions on coupled 1D-3D flow problems. Finally, the authors would like to thank M. Kuchta, K. E. Holter and C. Daversin-Catty for developing and sharing with us code for the implementation of mixed-dimensional models in FEniCS.