Well testing of radial jet drilling wells in geothermal reservoirs

Radial Jet Drilling (RJD) is a stimulation technique to increase a well’s performance through creation of multiple laterals of up to 100 m long. The technique has been used in the petroleum industry for several years and recently also for geothermal reservoirs. Interpretation of well tests of a well with multiple laterals may become problematic if the effect of the laterals is not correctly modelled. In this work it is investigated what the impact of RJD laterals is on the pressure transients for both single and dual porosity reservoirs. A semi-analytic model for the calculation of the transient well pressure that accounts explicitly for the RJD well geometry is developed and validated with a detailed numerical simulation. The results show that changes in the lateral configuration affect the pressure transients significantly. In particular for single porosity media, the laterals affect the pressure transients in ways that cannot be captured by representing the stimulation by a skin factor. Since the RJD process is unsteered and the exact configuration of the laterals uncertain, the model may potentially assist in estimating the effectiveness of the stimulation such as the achieved reach of the laterals. Although the work was carried out in the context of RJD application in geothermal reservoirs, the presented model approach can be extended to multi-lateral wells in oil or gas reservoirs.


Introduction
Wells in a geothermal system are not only the conduit through which fluids and heat are brought to surface, but also a major source of information on the subsurface. Static information on the properties of the reservoir rocks can be derived from numerous different types of well logs and cores. Dynamic information on the flow behaviour is primarily derived from well or flow testing (e.g. [1]) and, ultimately, from production behaviour. For most well tests, a controllable surface injection or production rate is changed and the pressure response in the well is measured. Through matching the pressure response with a mathematical (semi-analytical or numerical) model, certain reservoir and well parameters may be estimated [2]. Examples of such parameters are reservoir permeability, skin factor, reservoir dimensions and fracture characteristics in case of a naturally fractured reservoir. Matching well test results, requires a model that captures transient behaviour. For simple systems such as vertical wells in horizontally layered reservoirs, analytical solutions have been available for a long time [3]. For more complicated well geometries and reservoirs such as fractured reservoirs, deriving analytical solutions can be difficult.
Most existing geothermal wells are vertical or deviated for which analytical solutions for transient pressure analysis can be derived. Wells with low productivity, however, need to be stimulated, which is commonly done using acids or hydraulical stimulation [4]. Since chemical stimulation only affects the near well-bore area and hydraulic stimulation carries certain risks (in particular for geothermal applications, see [5]), other ways of stimulating wells or increasing production are required. Alternative well geometries are used increasingly for geothermal applications in which productivity needs to be enhanced. This is evidenced by a range of recent projects and papers: multi-lateral wells were drilled in Indonesia [6,7] and evaluated for an Enhanced Geothermal System (EGS) [8]; sub-horizontal geothermal wells were completed in Switzerland in 2013 [9] and a longer one in France in 2017 [10]. Another stimulation technique is Radial Jet Drilling (RJD) (e.g. [11][12][13]). The technique has been applied successfully in oil and gas wells for several years and is being extended to geothermal applications [12,14]. Stimulation with RJD was tested in a geothermal well in Klaipėda in Lithuania in 2015 ( [4,15,16], in The Netherlands [17] and in Iceland [18]. With RJD, several small-diameter laterals of length up to 100 m can be created in a vertical or deviated main well bore. The laterals are created using hydraulic jetting: a jetting nozzle is deployed on coiled tubing and via a deflecting shoe exits the well at a right angle. The achieved diameter depends on the rock properties and is generally in the order of a few centimetres. Since the jetting nozzle is not steered, the position and reach of the jetted laterals is uncertain. Typical well configurations resulting from RJD are 8 to 16 laterals which are grouped in kick-offs with 4 laterals each (Fig. 1). The resulting well configuration is challenging to simulate because of the complex geometry for both reservoir simulation [19] and well test interpretation. To our knowledge, semi-analytic models used in well test software are not capable of incorporating explicitly the complex RJD well geometry. The standard approach is to use a vertical or deviated well with a (negative) skin to incorporate well stimulation effects. However, given the extent of laterals (in the order of 100 m), such an approach might be too coarse for an accurate analysis of the early transients. In particular for fractured reservoirs which are common in geothermal reservoirs, early transients in well tests are very important because they can give information on the relative contributions of the fractures and the matrix [1,20].
In this paper, the impact of the lateral configuration on pressure transients is investigated for both fractured and nonfractured reservoirs. It is analysed what the benefits are of including the well geometry explicitly rather than as a skin factor in well test analysis. Since properties such as the length or position of the jetted laterals are uncertain, we will investigate whether it is possible to distinguish between lateral properties in the pressure transients. If this is possible, a well test might be used to improve the estimate of the lateral dimensions.
It is possible to account for the well geometry in a well test analysis numerically (e.g. [21,22,23], but this a comparatively large effort and lacks flexibility in varying well geometry. To simulate the pressure transients for a large range of reservoir conditions and RJD well configurations, a fast and flexible simulation approach is required. Therefore, a semianalytic approach is used to derive the full pressure transient solution for an RJD well in a dual-porosity reservoir including wellbore storage. The solution is for single-phase flow and isothermal conditions, thus most appropriate for production testing of low to medium temperature wells. The semi- Fig. 1 Illustration of an RJD well configuration analytic approach that is used is the Analytical Element Method (AEM). The AEM is used in [24] for steady state performance calculations for an RJD well and is extended in the present work to a transient solution. An advantage of the followed approach is that it allows for including wellbore storage and naturally extends to dual porosity reservoirs. The semi-analytic approach is validated with a numerical reservoir model with a very fine grid allowing to represent the well explicitly.

Single porosity
To calculate the full pressure transient solution the diffusivity equation needs to be solved with appropriate reservoir and well boundary conditions. In Eq. (1), η ¼ k φμc t is the diffusivity coefficient. Similarly as the steady state pressure solution for RJD wells was derived in Egberts and Peters [24], we will construct a solution for the pressure field from 3D point source solutions exploiting the principle of superposition [25]. One may derive a transient point source solution solving Eq. (1), however we choose to solve the problem in the Laplace domain as the formulation and solution simplifies. This comes at the additional cost of the numerical Laplace inversion to arrive at a pressure solution in the time domain. Before applying the Laplace transform, Eq. (1) is rephrased using dimensionless time, space and pressure variables. This is also convenient to generalize the solution approach to natural fractured reservoirs with a dual porosity description (Section 2.2). Defining the dimensionless quantities The Laplace transform e f of a function f is defined by with s a real or complex number. Applying the Laplace transform (3) to Eq. (2) gives a simpler equation as the time derivative becomes a multiplication in the Laplace domain: The full pressure solution in the Laplace domain will be constructed from point source solutions of Eq. (4). A point source solution of Eq. (4) is given by Here q p (m3/s) is the continuous rate at a point source location x ′ , e p D the pressure solution at a point x and r D = ‖x − x ′ ‖/r w with ‖•‖ the Euclidian norm. In a similar fashion as the steady state solution is constructed for an RJD well, in Egberts and Peters [24] the well is split up in a number of linear segments and for each segment a line source solution is derived through integration of point source solutions (Eq. (5)) along the segment, thereby allowing for a polynomial varying flux q along the segment. Thus, parametrizing a segment by a line l(τ) with -1≤τ ≤ 1 such that l(−1) and l(1) are the begin and endpoint of the segment, and assigning the point sources along the segment line l(τ) with q p τ ð Þ ¼ ∑ where we have defined r D (τ) = ‖x − l(τ)‖/r w . This integration is done numerically as no analytic solution was found. In this way a line source solution is calculated for each well segment having an independent set of polynomial coefficients. Using the principle of superposition, a pressure solution in the Laplace domain can be formed for the full RJD well by summing up the pressure solutions of the segments. This solution has multiple degrees of freedom, namely the coefficients of the polynomials describing the flux along the well segments, that can be exploited to satisfy a uniform pressure condition at distributed points along the well face. No-flow conditions at the top and bottom of the reservoir are obtained using the method of images. The reservoir is assumed to be described with homogeneous (average) permeability. Anisotropy is handled by transforming the reservoir to an equivalent isotropic reservoir by a volume conserving scaling of the space variables and as homogeneous permeability the geometric mean of the three directional permeabilities. The construction of a pressure solution in the Laplace domain requires inversion to the time domain in which well test data interpretation is usually done although well test interpretation may be performed in Laplace space [26]. A numerical Laplace inversion algorithm often used for well test interpretation is the Gaver-Stehfest method. This method expresses a function value at time t as a finite sum of weighted transformed function values, evaluated for real values s see e.g. Abate and Whitt [27]. Typically, a function value can be calculated by 8-14 transformed function evaluations using double precision arithmetic. This method works fine if the function has no discontinuities or singularities otherwise it may be inaccurate ( [28,29], Ajmi et al. 2008). It was indeed found that for injection tests to be an efficient algorithm, however for a well test consisting of a period of constant injection followed by a shut-in the Gaver-Stehfest algorithm fails, instead an inversion approach by Den Iseger [28,30] may be used (see Appendix B Fig. 10).

Dual porosity
Fractured reservoirs are systems in which the flow is characterized by main flow in the fractures and the main fluid storage in the rock matrix [3]. This conceptual model of a fractured system with primary porosity in the matrix and secondary porosity in the fracture system as presented by Warren and Root [20] is still mainly used today [1], although more complicated models with triple porosity have been proposed for shale [31]. For dual porosity reservoirs, flow in the matrix is exclusively from the matrix to fractures. For dual porositydual permeability reservoirs, flow can also occur within the matrix. Reservoirs with dual porosity behaviour are common in geothermal systems, in particular in volcanic rocks [32] or Enhanced Geothermal Systems [33]. Some (brittle) carbonate reservoirs and reservoirs in which the flow is dominated by high-permeability streaks can also exhibit dual porosity behaviour. For interpretation of well tests, accounting for the distribution of the storage over these different systems is important, because the transient behaviour of flow in a well depends strongly on this [20]. Also the exchange of fluids between the fractures and matrix, which depends mainly on the geometry and density of the fractures and the permeability of the matrix to the fractures, has a strong impact on the transient behaviour important for well test interpretation.
Not all aspects of flow in fractured media are accounted for in a dual porosity continuum approach. For example, flow near the well exhibits a linear rather than a radial flow regime as assumed in the solution of the diffusivity equation. This means that the pressure drop will be smaller than expected from the semianalytical solution, which will result in a negative skin [34,35]. In the Warren and Root model, the flow from the matrix blocks to the fractures is assumed to be in pseudo-steady state [1]. Other approaches assume a transient flow between matrix and fractures [36], which is not included in this paper.
It will be shown that the chosen solution approach described in Section 2.1 can rather straightforwardly be adopted for solving the flow equation for a dual porosity reservoir by constructing a solution from point source solutions in the Laplace domain. We follow the common dual porosity description of Warren and Root [20]. They derive a flow equation similar to the equation for single porosity with an additional matrix-fracture transmissibility term assuming pseudo steady state flow condition between matrix and fractures. This flow condition describes a flow between fractures and matrix restricted by skin giving a negligible pressure gradient in the matrix [1]. The dual porosity model with unrestricted inter-porosity flow assuming transient matrixfracture flow is not considered in this paper. Warren and Root [20] introduced two dimensionless parameters that characterizes the dual porosity system, that is, the storativity ratio and the inter-porosity flow parameter Values for the storativity ratio ω are generally below 0.1 or even 0.01. The higher values occur in cases were matrix porosity is very small as in crystalline rocks, or when fracture or secondary porosity is large as for example in case of highperm streaks or dissolution effects [1,3].
The matrix to fracture transmissibility factor σ can be expressed as σ = αk mb . Where α is a purely geometrical factor describing the fracture network. Density or fracture spacing is the most important factor influencing α. Values for the interporosity flow parameter or coefficient λ thus depend on the ratio of matrix over fracture permeability. This results in a wide range of values for λ since both matrix permeability and fracture density are highly variable.
Equations for dual porosity (see Appendix A) can be conveniently written down using the following dimensionless variables and Analogue to Eqs. (4) and (5) for single porosity the Laplace transformed equation for a dual porosity and its point source solution can be derived, see Eqs. (20) and (22) in Appendix A. The same solution approach for calculating the transient flow for single porosity reservoir can be generalized to dual porosity.

Wellbore storage
For early times, wellbore storage effects should be accounted for to perform a proper well test interpretation [2]. A well test for which the rate at surface conditions is taken constant will not instantly result in a constant rate at reservoir conditions due to the presence of a compressible liquid inside the well bore. Effectively this means that the reservoir pressure equation for the reservoir needs to be solved with a time varying pressure boundary condition at the well-reservoir interface. This requires the coupling of the wellbore and reservoir pressure equations.
We will assume that the controllable surface volume well rate q s w is taken constant. The reservoir volume well rate q w (at the reservoir depth) is time dependent due to wellbore storage effects and is given by where B w (reservoir Vol./surface Vol.) is the formation volume factor, V w the volume of the liquid in the well and c w the compressibility of the liquid inside the well bore. In dimensionless quantities, Eq. (11) converts into Where p wD is defined according Eq. 10 and The constraint that the well operates at a reservoir volume rate is in the AEM translated into a constraint through Eq. (13) for the unknowns describing the polynomial-shaped well influx along the well segments.

Base case RJD well
The base case for the validation of the transient AEM implementation is an RJD well consisting of a fully penetrating vertical well with 4 orthogonal laterals of length 100 m each with the same kickoff location at one-third of the well. The configuration is shown in Fig. 2. The radius r w of the vertical well (backbone) is 0.1 m and the radii of the laterals are 0.075 m. The well is producing during the test period at a constant surface volume rate of 3600 Sm 3 /day. It is further assumed that the pressure at the well face at reservoir depth is uniform. In case well bore storage (WBS) is included in the calculation a well volume V w of 500 m 3 is taken.
At the lateral boundary of the reservoir, a constant pressure is taken equal to the initial reservoir pressure. At the top and bottom of the reservoir a no-flow boundary condition is imposed. In Table 1 the reservoir and water properties are listed.

Comparison with a numerical reservoir simulator
The transient AEM is compared/validated with a reservoir simulation model with a fine scale grid permitting an explicit representation of the well in which the grid size is adjusted to the well diameter. The reservoir simulator used is the Eclipse® software [37]. The reservoir model and the well modelling is described in Appendix C Fig. 11 and Table 2.
The transient solution results of AEM for a single porosity reservoir are compared to a numerical reservoir simulation in Fig. 3a. In this figure the dimensionless pressure and its Bourdet derivative are shown in a log-log plot. These are typical curves used for analysing well tests, see e.g. [2,38]. For the Laplace inversion of the pressure build up, the Gaver-Stehfest algorithm was used, using 8 transformed values to approximate a pressure value at a single time. It was found that taking 8 transformed values is computationally efficient and gives sufficiently accurate results.
In Fig. 3b wellbore storage (WBS) is included in the simulations. Three time-regions can be identified from the pressure derivative plot; initially a region where WBS dominates followed by a transition period between the WBS region and the region where radial flow has developed. The radial flow time region starts when the pressure derivative flattens. For the three time-regions, two linear slopes can be identified. The linear slope for early times is the WBS effect giving a linear relation between well pressure and time, see e.g. [36]. The second linear slope is when radial flow is fully developed, but not hindered yet by boundary conditions (Infinite Acting Radial Flow).
We observe from Fig. 3 that the validity of the transient solution extends to earlier times for the AEM than for the numerical reservoir simulations. In Fig. 4 the pressure (a) and its Bourdet derivative (b) are compared for fixed interporosity flow parameter λ and varying storativity ratio ω. The S-shaped well pressure in the semi-log plot and the valley in the pressure derivative log-log plot as shown in the figure are typical for dual porosity reservoirs. This behaviour of the pressure will be discussed in Section 4.2.

Well-test response of RJD wells
In this section we present well test response curves obtained from the AEM simulations to illustrate the impact of the lateral configuration on the pressure transients and the benefit of using an explicit treatment of the RJD well geometry for interpretation of well tests results. Because the jetting process is not steered, the position of the laterals in the reservoir is highly uncertain [12]. By taking into account the laterals explicitly rather than as a (negative) skin of a vertical well, the well test response can be predicted more accurately and potentially information can be derived from the well test about the lateral geometry. For the simulations discussed below, we consider the base case RJD well as described in Section 3.1.

Single porosity
In Fig. 5a response curves for increasing lateral length are shown. The figure clearly shows that the development of radial flow, the time when the pressure derivative flattens, is later for longer laterals as expected. These response curves can assist in analysing well tests for RJD wells and can give independent information about how far the laterals have been jetted away from the well. Differences between a well test prior and after jetting can be interpreted to identify how far the laterals reached away from the well. It should be noted that the observed dependency requires an explicit treatment of the geometry of the well Fig. 2 The base case RJD well configuration  In addition to the reach of the laterals, also the number of laterals strongly impacts the well test response. To illustrate this, response curves have been generated for an increasing number of laterals, see Fig. 5b. In Fig. 6 the used RJD well configurations are shown.
The development of radial flow occurs nearly at the same time for all well configurations except for the well without laterals that shows a much earlier developed plateau of the pressure derivative. Obviously, the length of the laterals determines the timing of radial flow development rather than the number of laterals. Note also that the response curves of the two well configurations with two laterals are on top of each other.

Dual porosity
For better understanding of the pressure response of an RJD well in a dual porosity reservoir we make a comparison with the pressure behaviour of a vertical well. In Fig. 7 the transient well responses of a vertical and RJD well are shown for ω = 0.01 and λ varying from 10 −4 to 10 −7 .
The specific form of pressure and pressure derivative curves for the vertical well in a dual porosity reservoir are well-known [1,39]. Theoretically, for a vertical well, the pressure semi-log plot shows an S-shaped curve with two tangent lines for early and late times that are parallel with a slope equals 0.5 (Fig. 7, top left). This manifests in the pressure derivative plot (Fig. 7, bottom left) as constant value (=0.5) for early and late time. In practice the early time tangent line may not be observable as it may be obscured by wellbore storage effects especially for small λ. However, if it is there it may be used to estimate ω, as illustrated in Fig. 7 (top left).
The S-shaped pressure reveals three regimes, typical for a dual porosity reservoir with restricted inter-porosity flow [1,36]. In the first regime the flow in the fracture system is dominant and the reservoir acts as a homogeneous reservoir with fracture compressibility. In the second regime a transition  The fracture fluid pressure increases (assuming injection), the matrix will act as fluid storage and the matrix fluid pressure starts to equilibrate with the fracture fluid pressure. Fracture and matrix pressure at later times equalizes and the system will behave as a homogenous reservoir with total compressibility.
The RJD well does not show a parallel line at early times in the pressure semi-log plot shown in Fig. 7 (top right). The pressure plot does reveal a tangent line of slope 0.1 at very early times (t D < 0.01) when the RJD well act as five individual wells as interference is not yet developed. As, in this example, lateral length and well backbone length are equal, the pressure slope is 0.1, that is, 1/5 of the slope of a single vertical well. The deviation of the linear slope starts as soon as the well segments starts to interfere.
The typical occurrence of a valley in the pressure derivative plot for dual porosity reservoirs for a vertical well is also clearly seen for the RJD well (Fig. 7, bottom left and right). The pressure derivative log-log plot shows for the RJD well besides a valley also a peak. The depth of the valley in the pressure derivative for the RJD well appears to be depending on λ in contrast to a vertical well; the times associated to the depth of the valley are comparable.
For smaller λ, that is, less matrix-fracture flow coupling, a later onset of the late time behaviour is seen. For a vertical well, late time behaviour occurs typically for t D ≫ 1 λ . For larger λ, the start of late time behaviour is later for an RJD well as can be seen by comparing the pressure derivative curves for λ = 10 −4 . The late time pressure behaviour (the pressure derivative value) of an RJD well is the same as for a vertical well because the RJD can then be considered as a vertical with a (negative) skin.
We observe from Fig. 8 that the onset of late time behaviour takes place for t D ≫ 1 λ ¼ 10 6 , independent of ω and is the same for the vertical and RJD well. For a vertical well, the parallel early and late time straight lines may be used to estimate ω [36]  For the considered values for λ and ω, the early homogeneous behaviour is less noticeable for an RJD well; only for a short time a constant slope can be seen in the response that is not parallel to the late time tangent. The response curves for an RJD well for the varying parameters are less congruent than seen for the vertical well, only for the late time the curves are similar.
The effect of the RJD well design on the well response is investigated for a dual porosity reservoir with λ = 10 −6 and ω = 0.01 in Fig. 9 in which pressure and its derivative is shown for (a) varying lateral length for the base case RJD configuration and (b) varying number of laterals. The derivative curves give hardly any distinction for later times; the chosen λ and ω determine the shape. In early times, there are differences due to the varying configuration of the well, but these are much more difficult to detect in real well testing conditions due to well bore storage effects and noise.

Conclusions
In this paper, the impact on pressure transients of laterals created using Radial Jet Drilling (RJD) is investigated for both fractured and non-fractured reservoirs. It is analysed what the benefits are of including the well geometry explicitly for well test analysis. To achieve this, an Analytic Element Method (AEM) was developed to derive the full pressure transient solution for an RJD well including wellbore storage effects for estimating well and reservoir characteristics through interpretation of well test data. The formulation in Laplace domain was found convenient to describe the solution for single and dual porosity reservoirs in a uniform manner. The semianalytic approach was successfully validated with a reservoir simulation model with a very fine grid with detailed representation of the well.
Using the AEM approach, it was shown that the details of the RJD lateral configuration have impact on the pressure transients and that incorporating the explicit geometry of an RJD well when analysing the early pressure transients has benefits compared to a simplified approach based on skin. In particular for single porosity flow, a clear impact of the length of the laterals on the pressure transients could be identified: Verical well RJD well Fig. 7 Pressure and pressure derivative for a vertical well (left) and an RJD well (right) in a dual porosity reservoir. ω = 0.01 is fixed and λ varying from 10 −4 to 10 −7 . No wellbore storage included longer laterals impact the pressure transients for much longer periods than short ones. Since the reach of the laterals is uncertain due to the jetting process, a well test can therefore provide information on the achieved reach if the pressure transients are analysed using the AEM approach with explicit lateral geometry. This cannot be captured in a pressure transient analysis where skin is used to represent the laterals. A simplified approach with skin thus limits the possibility to Vertical well RJD well estimate RJD well characteristics and thus to test the effectiveness of the jetting of the laterals. For dual porosity reservoirs, the impact of differences in lateral configurations on the pressure transients is more difficult to identify. Due to the strong imprint of the characteristics of the dual porosity system, the difference caused by for example the length of the laterals is smaller. Therefore it is possible that for real well test data with noise the impact of the lateral configuration cannot be identified.
Although the AEM approach was developed and tested for RJD wells in this paper, the method can also be used to analyse other multi-lateral well configurations and as such provides a valuable addition to currently available well test interpretation tools.
The standard Gaver-Stehfest Laplace inversion, often applied in (semi-) analytic well models, is erroneous for tests with discontinuities in the well rate. The pressure response is not correctly inverted. This is demonstrated for a pressure falloff test after a well shut-in of an injector. A much less known and more recent inversion method by Den Iseger can perform the inversion accurately.

APPENDIX A: DUAL POROSITY PRESSURE EQUATION
The flow equation described by Warren and Root [20] is given by Here η fb ¼ k fb φ fb c tf μ and the extra term accounting for matrix-fracture flow is given by σ mf ¼ σ μ ðp m −p f ). The equation is complemented with the mass balance equation in the matrix Introducing the storativity ratio (7), the inter-porosity flow parameter (8) and the dimensionless variables (16) the equations can be conveniently rephrased as and Applying the Laplace transform (3) to Eq. (16) and (17) gives the following simpler equations as the time derivative becomes a multiplication in the Laplace domain: Eliminating e p m in (18) using (19) gives with p fD the dimensionless fracture pressure and Note that the single porosity description is a special case of the dual porosity description as for ω = 1 and f(s, ω, λ) = 1 Eq. (20) reduces to Eq. (4).
The full pressure solution in Laplace domain will be constructed from point source solutions of Eq. (20): in a similar way as described in Section 2.1 for single porosity.
Here q p (m3/sec) is the continuous rate at a point source location x ' , e p fD the solution at a point x and r D = ‖x − x ′ ‖/r w .

APPENDIX B: NUMERICAL LAPLACE INVERSION
The usage of the Gaver-Stehfest Laplace inversion algorithm was found accurate for the injection tests discussed in the previous sections in which the pressure increase is modelled resulting from injection with a constant rate. Well tests with pressure response of a well that is shut-in after a period of injection of production may sometimes be more desirable as the flow rate is zero and thus precisely known while a constant rate unequal zero may be difficult to realize. Other well test that may be relevant is where the well rate is stepwise changed thereby avoiding shut-in and hence less production loss. These tests have in common discontinuities in the well rate that cannot be handled by the Gaver-Stehfest algorithm. An alternative Laplace inversion algorithm by Den Iseger [28,30] can deal with discontinuities. The method calculates the function values at a predefined equidistant time grid (t k + 1 = t k + Δ) at once with increasing accuracy for smaller grid size Δ. This differs from the Gaver-Stehfest algorithm that can perform inversion at arbitrary time points but one at a time.
Also, in contrast to the Gaver-Stehfest algorithm, the inversion by Den Iseger requires the evaluation of the Laplace transformed function for complex Laplace variable values. For our implementation of the solution approach this poses no problem as the validity of the solution in Laplace space extends to the complex plane.
The Laplace inversion method of Den Iseger uses the Poisson summation formula that relates a Fourier series of function values to an infinite sum of transformed function values and as such falls under the Fourier series method class of inversion methods. See [27,40] for a discussion on Laplace inversion methods. Den Iseger applied a Gaussian quadrature rule to approximate the (slow converging) infinite sum of transformed values by a finite sum. Application of the FFT to the finite sum of transformed values gives the function values [28,30].
Both Laplace inversion methods are compared for a falloff test of the base case RJD well that is injecting for 2 days with constant rate and then shut-in. To simulate a pressure falloff test, i.e. pressure build up (with constant injection rate q) followed by a pressure falloff due to shut-in is a well is injecting with the same q also beyond shut-in time and an identical well (same configuration and location) is introduced coming on stream at moment of shut-in which produces with rate −q. By superposition the sum of the solution gives the correct pressure field as the sum of the rates is zero after shut-in time.
In Fig. 10 the two inversion methods are compared, for reference the fine scaled Eclipse simulation result is shown as well (green line). Gaver-Stehfest inversion result (red line) is shown with tuning parameter value N = 8 (equals the number of summed transformed values). Gaver-Stehfest fails to perform an accurate Laplace inversion for the pressure falloff. Increasing values of N does not improve the solution, in fact too large N values gives non-monotonous results due to numerical errors. Using larger N would require higher (than double) precision arithmetic which will be computationally infeasible for our method. Applying Gaver-Stehfest independently for the pressure build up and falloff period yields the best results with a large N = 20 (red dots); the pressure buildup period is accurate in agreement with the findings in Section 3 but it fails for the falloff period. The observed erratic behaviour Gaver-Stehfest, of the Laplace inversion algorithm, has also been reported in [29]. The inversion method by Den Iseger [28,30] gives indeed accurate results both for the buildup and falloff phase. (dashed blue line), see also [29].

APPENDIX C: NUMERICAL RESERVOIR MODEL
The configuration of the RJD well consists of a fully penetrating vertical well (backbone) with 4 orthogonal and horizontal laterals at 33 m from the top. The vertical backbone and the laterals are modelled explicitly using fine grid blocks with high permeability (for details see the third paragraph of this appendix). The size of the grid blocks representing the well determines the well diameter and thus must be taken small resulting in a large number of grid blocks for the entire grid. To reduce CPU time, only ¼ of the near well area is simulated, as indicated in blue in Fig. 11. Note that the flow in this part of the model domain is symmetrical to the flow in the other three parts. From the backbone ¼ of the diameter is simulated and for the laterals half the diameter.
The first layer of the model is reserved for representing the well above the reservoir. Only the grid block connected to the well is active. All other grid blocks are inactive. In this active grid block a sink is located to allow the flow to leave the model. The grid block containing the sink is only connected to the reservoir via the grid blocks representing the well. The well representation will be explained further below. Both the backbone and laterals are simulated explicitly with grid blocks of 0.10 × 0.10 m. Thus, the diameter of the backbone is 0.20 m (four grid blocks of which 1 is simulated). The diameter of the laterals is 0.15 m (two grid blocks of 0.10 × 0.10 m of which 1 is simulated). The transmissibility between the sink and the block is taken very large to avoid pressure drop. The grid block containing the sink can be used to simulate wellbore storage by increasing the block volume. For simulations with wellbore storage, the pore volume of the well block is multiplied by 1000.
Permeability in the gridblocks representing the well need to be set large to ensure low pressure drop in the well. For comparison with the results of AEM, which does not take into account pressure drop in the well, the pressure drop in the well is minimized. For that, the permeability in the well grid blocks is set to 10 million D vertically and 1 million D horizontally. In the top part of the backbone, where the rate is 3600 Sm 3 /d, the pressure drop is 56 Pa/m. In the backbone, the assumption of no pressure drop is normally justified, unless the tubing has scaling or corrosion. For the laterals, the pressure drop is highly uncertain, because the rate in the laterals and roughness of the laterals are uncertain. Earlier analysis showed that pressure drop is likely to be important in the laterals [14]. However, the impact of potential pressure drop in the laterals has not been investigated in this paper.
The boundary conditions of the model are set as follows: A constant pressure at 758 m and a no-flow top and bottom boundary condition.
Different grid size distributions were tested in both horizontal and vertical direction, but this had little impact on the pressure drop and well test results. In Table 2 the applied numerical grid is specified.
Reservoir and fluid properties are listed in Table 1. The fluid properties are consistent with a brine of 150.000 ppm at 75°C and 25 MPa. Both viscosity and compressibility are constant with pressure. The well flow rate used in the model is 900 Sm 3 /d to take into account that only ¼ of the well and reservoir is simulated.
The results are sensitive to time stepping: due to the large contrasts in permeability, the time steps have to be relatively small during the entire period to avoid numerical instability. At the start, time steps have to be very small to capture the detailed pressure response near the well. Time step size starts at 0.864 s and grows to a maximum time step size of 110.6 s after 0.686 days.
Acknowledgements 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.
Data availability The data that support the findings of this study are available from the corresponding author, Egberts P., upon reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Fig. 11 Overview of the model setup with the RJD well consisting of a vertical well with 4 perpendicular laterals