** **
28 Aug 2018

28 Aug 2018

# Modelling of multi-lateral well geometries for geothermal applications

Elisabeth Peters Guido Blöcher Saeed Salimzadeh Paul J. P. Egberts and Mauro Cacace

^{1},

^{2},

^{3},

^{1},

^{2}

**Elisabeth Peters et al.**Elisabeth Peters Guido Blöcher Saeed Salimzadeh Paul J. P. Egberts and Mauro Cacace

^{1},

^{2},

^{3},

^{1},

^{2}

^{1}TNO Applied Geosciences, Utrecht, 3584 CB, the Netherlands^{2}Helmholtz-Zentrum Potsdam – Deutsches GeoForschungsZentrum GFZ, Telegrafenberg, 14473 Potsdam, Germany^{3}Danish Hydrocarbon Research and Technology Centre (DHRTC), Lyngby, Denmark

^{1}TNO Applied Geosciences, Utrecht, 3584 CB, the Netherlands^{2}Helmholtz-Zentrum Potsdam – Deutsches GeoForschungsZentrum GFZ, Telegrafenberg, 14473 Potsdam, Germany^{3}Danish Hydrocarbon Research and Technology Centre (DHRTC), Lyngby, Denmark

**Correspondence**: Elisabeth Peters (lies.peters@tno.nl)

**Correspondence**: Elisabeth Peters (lies.peters@tno.nl)

Received: 29 May 2018 – Revised: 07 Aug 2018 – Accepted: 12 Aug 2018 – Published: 28 Aug 2018

Well inflow modelling in different numerical simulation
approaches are compared for a multi-lateral well. Specifically radial wells
will be investigated, which can be created using Radial Jet Drilling (RJD).
In this technique, powerful hydraulic jets are used to create small diameter
laterals (25–50 mm) of limited length (up to 100 m) from a well. The
laterals, also called radials, leave the backbone at a 90^{∘} angle. In
this study we compare three numerical simulators and a semi-analytical tool
for calculating inflow of a radial well. The numerical simulators are FE
approaches (CSMP and GOLEM) and an FV approach with explicit well model
(Eclipse^{®}). A series of increasingly complex
well configurations is simulated, including one with inflow from a fault.
Although all simulators generally are reasonably close in terms of the total
well flow (deviations < 4 % for the homogeneous cases), the
distribution of the flow over the different parts of the well can vary
significantly. Also, the FE approaches are more sensitive to grid size when
the flow is dominated by radial flow to the well since they do not include a
dedicated well model. In the FE approaches, lower dimensional elements (1-D
for the well and 2-D for the faults) were superimposed into a 3-D space. In
case the flow is dominated by fracture flow, the results from the FV approach
in Eclipse deviates from the FE methods.

Geothermal resources can vary widely in their characteristics, depending on
the downhole temperature, water composition, depth and the reservoir rock
(see e.g. http://geothermalcommunities.eu). For all these applications,
production occurs via geothermal wells, which form the link between the
subsurface and surface. Planning and operation of the geothermal resource
relies on predictions of the performance of the geothermal wells, which is
generally done using numerical methods. Different numerical approaches can be
used such as finite element (FE), finite difference (FD) or finite
volume (FV) or combinations thereof, which rely on different discretisations
of the subsurface domain. FE approaches are particularly suitable for
mechanical simulations and are generally combined with control volumes for
flow simulation (Zienkiewicz and Taylor, 2000). FD is the oldest method and
mostly used for flow simulation. In the FV approach, the domain is
discretised in volumes rather than distances as in FD, which makes it
strongly mass-conservative. This method is particularly suitable for flow
modelling on unstructured grids (Eymard et al., 2003). Irrespective of the
numerical approach, accurate well inflow modelling is essential. With more
complex well geometries becoming possible, accurate simulation of well inflow
becomes more difficult (Wolfsteiner et al., 2003). In this study we focus on
the well inflow modelling of different simulation approaches for a specific
type of well, namely a radial well. This type of well can be created using
Radial Jet Drilling (RJD) (e.g. Kamel, 2017). With this technique, powerful
hydraulic jets are used to create small diameter laterals (2.5 to 5 cm) of
limited length (up to 100 m) from a well. The laterals leave the backbone at
a 90^{∘} angle and are called “radials”. Generally 10 to 15 of such
laterals are jetted in a single well. Currently the technique is best
suitable for sedimentary rock with porosity > 5 %, and most
implementations are in hydrocarbon reservoirs. For such a radial well
configuration, accurate well inflow modelling is more difficult, because of
multiple laterals being close together and close to the backbone and with
varying diameters. In this study we compare three numerical simulators and a
semi-analytical tool for calculating inflow into a radial well. The numerical
simulators are FE approaches (CSMP and GOLEM) and an FV approach with
explicit well model (Eclipse^{®}). The
semi-analytical tool is based on the Analytical Element Method (AEM). A
series of increasingly complex well configurations is simulated, including a
case with inflow from a fault. The first two cases which are simulated have
homogeneous reservoir properties, because these can be simulated by the
semi-analytical tool. In the third case, laterals are used to connect to a
fault. The laterals are simplified as straight wells although it is unlikely
that they are straight in the subsurface (Reinsch et al., 2018).

## 2.1 Description of the simulation codes

### 2.1.1 Semi-analytical approach (AEM)

The semi-analytic method to estimate the productivity (or injectivity) of a radial well as applied in this paper does not require a spatial discretization and is as such entirely distinctive from the other discussed numerical methods. To allow for a semi-analytic approach, certain simplifications are needed for the geometry and geology of the reservoir. We assume that the radial well is positioned in a cylindrically shaped homogeneous and anisotropic reservoir with laterally a constant pressure boundary condition and a no flow boundary condition at the top and bottom. Due to the complex geometry of a radial well it is likely not possible, even with the above simplifications, to obtain a closed formula for the well's performance (e.g. the Productivity Index that relates drawdown to production rate). The semi-analytic method is based on the AEM (Analytic Element Method) (see Egberts et al., 2013, and Egberts and Peters, 2015, for a more in-depth explanation of the method).

The semi-analytic method provides a solution of the pressure field
constructed from analytic solutions for the well segments (backbone and
laterals) of the radial well. An important component to arrive at the
semi-analytic solution is to model each well segment by a line sink (or line
source in case of injection) of the (steady state) pressure equation with a
variable influx along the line. The variable influx profile along a well
segment is described by a polynomial of an order *n*, typically less than 10.
The coefficients of the polynomials are then exploited to satisfy a
prescribed boundary condition at the well such as a constant bottom hole
pressure. This is done by distributing so-called control points along the
well bore face of the segments and requiring that the pressure solution
satisfies the bottom hole pressure at these points.

### 2.1.2 ECLIPSE

Eclipse^{®} is an industry-standard simulator used extensively in
the petroleum industry and to a lesser extent for geothermal applications
(Schlumberger, 2016). The simulator uses a finite volume approach and employs
a well model to estimate the pressure drop between the well to the grid block
in which it is located based on the approach of Peaceman (1983). This allows
for much larger (> 10 m) grid blocks than well radii
(∼ 0.1 m) and simulation of multiple wells in a single
reservoir model. A well is discretised in a number of well connections (well
nodes) with each connection associated to a grid block intersected by the
well. We will call a grid block intersected by the well a well block and the
intersection a well segment. The well index WI_{i} of well segment *i* is
defined as (e.g. Wolfsteiner et al., 2003):

with *μ* the viscosity. WI_{i} relates the difference of well pressure
*p*_{w,i} and gridblock pressure *p*_{b,i} to the flow rate *q*_{i} for a
segment *i* and depends on the length and diameter of the well segment and
the well block size and permeability. Also the orientation of the well with
respect to the axis of the grid influences the value of the WI. The WI can be
calculated by Eclipse from given well geometry and grid properties or can be
calculated in a pre-processing step and be provided to Eclipse. The latter
has been done in this paper, where the WI is calculated by
Petrel^{®}, a package for geological modelling.
In a well block *i*, at the point of the backbone where the laterals emerge,
the well indices WI_{j,i} of the different well segments, labeled by
index *j*, need to be amalgamated to a single total well index
(WI_{t,i}). In this case WI_{t,i} is calculated as:

### 2.1.3 GOLEM

The numerical simulator GOLEM is an open source software for solving parallel tightly coupled non-linear THM processes in fractured reservoir (Cacace and Jacquey, 2017). It is based on the MOOSE framework (Gaston et al., 2009) and its internal architecture relies on state of the art libraries for finite element analysis (libMesh, Kirk et al., 2006) and nonlinear iterative algebraic solvers (PETSc, Balay et al., 2016). GOLEM is based on the finite element method using irregular tetrahedral meshes. The major complication of using irregular meshes is to maintain the internal geometric consistency between well, fracture and matrix elements. This requires each well element to be mapped onto an edge (1-D representation) and each fracture/fault element mapped onto a face (2-D representation) of at least a single 3-D matrix element. In the current study, the well is simulated as a 1-D element superimposed in the 3-D domain, which consists of tetrahedrons. The refinement of the 1-D well-path strongly determines the pressure drop near the well and thus the productivity/injectivity of the well.

### 2.1.4 CSMP

Complex Systems Modelling Platform (CSMP) is an object-oriented application programme interface (API), for the simulation of complex geological processes and their interactions (formerly CSP, cf. Matthäi et al., 2001). CSMP is based on the finite element method (FEM) and finite element – finite volume (FE-FV) method. CSMP has been utilised to develop coupled thermo-hydro-mechanical-chemical (THMC) model for simulation of subsurface processes (Salimzadeh et al., 2017, 2018). Fractures/faults and wells are modelled using lower dimension elements, i.e. surface elements for fractures/faults and line elements for wells, while the rock layers are modelled using 3-D volume elements. In CSMP, the spatial discretization can be achieved using both structured and unstructured meshes, providing flexibility for different applications. Since the wells are modelled using 1-D elements, like in GOLEM the pressure response is highly mesh dependent. To remove this well-known issue, the size of the tetrahedral connected to the well has been reduced such that the position of the nearest integration point, to the well, where the integrand is being evaluated numerically, resembles the well radius. With this strategy, CSMP results provide a very good fit to the analytical results for the simple geometry wells as discussed later in the results section. A basic well model based on the Peaceman approach (Peaceman, 1983) is available within CSMP, but was not used in this study.

## 2.2 Description of the cases

### 2.2.1 Case 1

The first case is a vertical well with a single kickoff with four orthogonal
and horizontal laterals in a homogeneous, anisotropic reservoir of 100 m
thick (from 2500 to 2600 m depth) (see Fig. 1). The initial pressure is
25 MPa (at 2500 m depth) with a hydrostatic pressure distribution in the
vertical. Horizontal permeability is 200 mD, vertical permeability is 20 mD
and porosity is 0.2. The lateral boundary condition is a constant pressure
boundary on the edge of the model, which is at 1000 m from the well. Top and
bottom of the reservoir are no flow boundaries. The reservoir fluid is water
with a salinity of 150.000 ppm, which
has, at reference pressure, a viscosity of 0.54 cP and density of
1110.2 kg m^{−3}. The simulations are isothermal. The well is operated as
an injector with a fixed injection pressure of 26 MPa @ 2500 m depth. Thus,
the steady state flow rate for a pressure difference of 1 MPa is calculated.
Because gravity and water and rock compressibility are not accounted for in
AEM, the settings of the numerical simulators has to be adjusted in order to
be comparable. Therefore compressibility is set extremely low: rock
compressibility is 10^{−9} 1∕MPa and water compressibility
10^{−6} 1∕MPa and the gravitational acceleration is set to 0.

For GOLEM an element size of 0.32 m for both the backbone and the laterals was selected, because this gave results that most closely matched the semi-analytical solution. For CSMP an element size of 1.5 m for the backbone and 0.2 for the radials was selected for the same reasons. Grid size for the Eclipse model was 10 × 10 × 2.5 m and was not tuned to the results of the semi-analytical solution. The well and laterals in the Eclipse model are defined in the middle of the grid blocks were possible.

### 2.2.2 Case 2

For case 2, the reservoir is identical to the reservoir in case 1. Only the
well configuration is changed to a deviated well (see Fig. 1). The backbone
has an inclination of 35^{∘} with respect to the vertical. The laterals
leave the backbone at a 90^{∘} angle. This case is numerically more
challenging because the well is not aligned with the main direction of the
anisotropic permeability and for Eclipse is not aligned to the grid.

### 2.2.3 Case 3

For case 3 a low permeability reservoir is selected with a single large
fault. The case is based on the Groß Schönebeck case and details of
the reservoir properties and the reservoir fluids can be found in
Blöcher et al. (2010) and Blöcher et al. (2015). For simplicity
reasons only 1 fault (F21n) and 1 well were considered for the simulation.
The well itself has 8 radials with 2 kickoff points at −4000 and
−4050 m, respectively (Fig. 2). Each kickoff point has 4 radials with a
length of at least 100 m each. The two radials in the direction of the fault
are increased in length to ensure they intersect the fault. Although
extensions of more than 100 m are currently not technically feasible, this
scenario has been assumed to test the simulation capabilities. For the fault,
two scenarios are defined with different transmissivity to mimic: (Scen. 1) a
highly conductive fluid pathway and (Scen. 2) a hydraulically invisible fault
zone. The transmissivity *T* (Table 1) is obtained by the product of the
permeability *k* and the cross-sectional flow area *A* = *w*⋅*a*,
where *w* denotes the fault width (assumed to be 1 m) and *a* denotes the
fault aperture. The permeability of fault is approximated using the
Hagen-Poiseuille solution of the Navier Stokes equation for laminar flow
between two parallel plates separated by a constant aperture *a*:

which expresses the so-called cubic law (Witherspoon et al., 1980).

## 3.1 Case 1 and 2

In Table 2, the calculated flow is presented for the vertical well with four
horizontal laterals. The relative difference is reported with respect to the
semi-analytical tool AEM. The flow rates obtained by AEM, Eclipse and CSMP
are 3241.8, 3214.4, and 3161.0 sm^{3} d^{−1}, respectively. AEM,
Eclipse and CSMP are within 2.5 % from each other in terms of the total
flow. The results obtained by GOLEM (3374.0 sm^{3} d^{−1}) shows a
4.1 % larger flow rate in comparison to AEM. Table 3 summarizes the
results for case 2. For the deviated well and laterals the difference between
the simulators becomes larger (up to 21 %), although the total flow of
the well is within 3 % of each other for all simulators. Figures 3 and 4
show the inflow into the downward pointing (track 13) and horizontal lateral
(tracks 12 and 14) respectively. The inflow by the numerical simulators is
not always smooth as a result of the spatial discretization (e.g. GOLEM and
Eclipse for the downward pointing lateral in Fig. 3), which can give rise to
inaccuracies. Table 4 shows the increase in flow rate (for a fixed pressure
difference) achieved as a result of stimulation of the well with four
laterals. This is calculated by comparing the results of a run with and
without laterals. For case 1 the increase of flow rate predicted by AEM is
57.4 %. Only the estimated increase by CSMP is somewhat lower at
54.1 %. For case 2, the increase of flow rate predicted by AEM is
60.4 %. The numerical simulators Eclipse and GOLEM are quite close with
an increase in flowrate of 63.3 %, 57.5 %, while CSMP deviates a bit
more at 65.0 %.

## 3.2 Case 3

For case 3, only Eclipse and GOLEM are used to simulate the flow, because the semi-analytical tool cannot handle fracture flow and the results of CSMP and GOLEM are very similar if the same gridding is used. The differences between these simulators mostly arise from differences in the elements and refinement of the elements near the well and laterals. Overall the injectivity is lower for Eclipse than for GOLEM (Table 5). The difference between GOLEM and Eclipse is largest for the low permeability fault with laterals: 48 %. In general, connecting to a previously unconnected fault is highly beneficial for the injectivity in the well. It should be noted however, that for commercial rates a considerable number of laterals should be achieved since the diameter of the laterals are small and thus the flow through the laterals is limited.

Although all simulators generally are reasonably close in terms of the total well flow (deviations < 4 % for the homogeneous cases), the distribution of the flow over the different parts of the well can vary up to 20 % for some laterals. For the homogeneous cases (1 and 2), the predictions of increase of flow as a result of stimulation by RJD show a range of variation up to 5 % just from differences between numerical solutions even for a simple setup. In realistic implementations with heterogeneous reservoir properties, larger uncertainty from the numerical solution can be expected for all simulators: for Eclipse because of inaccuracies in the calculation of the well index and for the FE approaches because of difficulties in determining the correct mesh size and large number of elements. In case the flow is dominated by fracture flow (case 3), the results deviate more with up to 50 % difference in the predicted flow rate in the case of radials. Even though these uncertainties are considerably smaller than those arising from uncertainty in the properties and uncertainty in the radial path, it is a source of errors that is often ignored.

No data sets were used in this article.

EP wrote the largest part of the paper, defined cases 1 and 2 and did the Eclipse simulations. GB and MC conducted the FEM simulations performed with GOLEM and OGS. Furthermore case 3 was provided by GB and the related simulations were performed and illustrated. SS conducted the simulations performed with CSMP. PE coded and ran the AEM method.

The authors declare that they have no conflict of interest.

This project has received funding from the European Union's Horizon 2020
research and innovation programme under grant agreement no. 654662. The
content of this paper reflects only the authors' view. The Innovation and
Networks Executive Agency (INEA) is not responsible for any use that may be
made of the information it contains.
Edited by: Luke Griffiths

Reviewed by: Estanislao Pujades and one anonymous referee

Balay, S., Abhyankar. S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp, W. D., Kaushik, D., Knepley, M. G., May, D. A., Curfman McInnes, L., Tran Mills, R., Munson, T., Rupp, K., Sanan, P., Smith, B. F., Zampini, S., and Zhang, H.: PETSc Users manual, Tech. rep., ANL-95/11 – Revision 3, Argonne National Laboratory, available at: http://www.mcs.anl.gov/petsc (last access: 27 August 2018), 2008.

Blöcher, M. G., Zimmermann, G., Moeck, I., Brandt, W., Hassanzadegan, A., and Magri F.: 3D numerical modeling of hydrothermal processes during the lifetime of a deep geothermal reservoir, Geofluids, 10, 406–421, https://doi.org/10.1111/j.1468-8123.2010.00284.x, 2010.

Blöcher, G., Cacace, M., Reinsch, T., and Watanabe, N.: Evaluation of three exploitation concepts for a deep geothermal system in the North German Basin, Comp. Geosci., 82, 120–129, https://doi.org/10.1016/j.cageo.2015.06.005, 2015.

Cacace, M. and Jacquey, A. B.: Flexible parallel implicit modelling of coupled thermal-hydraulic-mechanical processes in fractured rocks, Solid Earth, 8, 921–941, doi10.5194/se-8-921-2017, 2017.

Egberts, P. J. P. and Peters, E.:A Fast Simulation Tool for Evaluation of Novel Well Stimulation Techniques for Tight Gas Reservoirs, Society of Petroleum Engineers, https://doi.org/10.2118/174289-MS, 2015.

Egberts, P. J. P., Shatyrbayeva, I., and Fokker, P. A.: Well Inflow Modelling for Wells not Aligned to a Numerical Grid. Society of Petroleum Engineers, https://doi.org/10.2118/165986-MS, 2013.

Eymard, R., Gallouët, T., and Herbin, R.: Finite Volume Methods. 2003. Update from Handbook of Numerical Analysis, P. G. Ciarlet, J. L. Lions eds, vol 7, 713–1020, available at: https://old.i2m.univ-amu.fr/~herbin/PUBLI/bookevol.pdf (last access: 7 August 2018)

Gaston, D., Newman, C. Hansen, G., and Lebrun-Grandié, D.: MOOSE: A parallel computational framework for coupled systems of nonlinear equations, Nucl. Eng. Design, 239, 1768–1778, 2009.

Kamel, A. H.: Radial Jet Drilling: a technical review, Society of Petroleum Engineers, https://doi.org/10.2118/183740-MS, 2017.

Kirk, B. S., Peterson, J. W., Stogner, R. H., and Carey, G. F.: libMesh: A C$++$ Library for Parallel Adaptive Mesh Refinement/Coarsening Simulations, Eng. Comp., 22, 237–254, 2006.

Matthäi, S. K., Geiger, S., and Roberts, S. G.: The complex systems platform csp3.0: Users guide. Technical report, ETH Zürich Research Reports, 2001.

Peaceman, D. W.: Interpretation of Well-Block Pressures in Numerical Reservoir Simulation with Nonsquare Gridblocks and Anisotropic Permeability, Society of Petroleum Engineers, https://doi.org/10.2118/10528-PA, 1983.

Reinsch, T., Paap, B., Hahn, S., Wittig, V., and van den Berg, S.: Insights Into the Radial Water Jet Drilling Technology – Application in a Quarry, https://doi.org/10.1016/j.jrmge.2018.02.001, 2018.

Salimzadeh, S., Paluszny, A., and Zimmerman, R. W.: Finite Element Simulations of Interactions between Multiple Hydraulic Fractures in a Poroelastic Rock, Int. J. Rock Mech. Min. Sci., 99, 9–20, 2017.

Salimzadeh, S., Paluszny A., Nick, H. M., and Zimmerman, R. W.: A three-dimensional coupled thermo-hydro-mechanical model for deformable fractured geothermal systems, Geothermics, 71, 212–224, 2018.

Schlumberger: Eclipse (R) version: Industry standard reference manual, 2016.2, 2016.

Witherspoon, P. A., Wang, J. S., Iwai, K., and Gale, J. E.: Validity of cubic law for fluid flow in a deformable rock fracture, Water Resour. Res., 16, 1016–1024, 1980.

Wolfsteiner, C., Durlofsky, L. J., and Aziz, K.: Calculation of Well Index for Nonconventional Wells on Arbitrary Grids, Comp. Geosci., 7, 61–82, 2003.

Zienkiewicz, O. C. and Taylor, R. L.: The Finite Element Method, Fifth Edition, Butterworth-Heinemann, Oxford, 707 pp., 2000.