POLITECNICO DI TORINO Repository ISTITUZIONALE An effective criterion to prevent injection test numerical simulation from spurious oscillations /

An Effective Criterion to Prevent Injection Test Numerical Simulation from Spurious Oscillations — Injection/fall-off tests are one of the most promising alternatives to the conventional production/build-up sequence because they eliminate surface emissions and can significantly reduce testing costs. This kind of test is characterized by the presence of twomobile phases, the fluid originally in place (hydrocarbon) and the injected fluid (Diesel, brine or nitrogen). The conventional analytical approach used to describe the transient pressure behavior is no longer suitable due to the variations in fluid saturations during the test. Although applicable in theory, the analytical approach often implies excessive simplifications of the real system behavior, such as piston-like displacement. Thus only numerical simulations can thoroughly describe the phenomena occurring during the injection process. However, the pressure and pressure derivative response calculated numerically often shows nonphysical oscillations during the radial flow phase, when the pressure derivative is expected to be horizontal. It was found that these spurious oscillations arise in convection-dominated problems and are associated with sharp saturation fronts. In this paper, an effective methodology, based on an adaptive time-step calculation, is presented so as to avoid pressure oscillations. The proposed time-step selection is both computationally efficient and suitable to capture the physics of the system. Oil & Gas Science and Technology – Rev. IFP Energies nouvelles, Vol. 69 (2014), No. 4, pp. 633-651 Copyright 2014, IFP Energies nouvelles DOI: 10.2516/ogst/2013137


Abstract -An Effective Criterion to Prevent Injection Test Numerical Simulation from Spurious
Oscillations -Injection/fall-off tests are one of the most promising alternatives to the conventional production/build-up sequence because they eliminate surface emissions and can significantly reduce testing costs. This kind of test is characterized by the presence of two mobile phases, the fluid originally in place (hydrocarbon) and the injected fluid (Diesel, brine or nitrogen). The conventional analytical approach used to describe the transient pressure behavior is no longer suitable due to the variations in fluid saturations during the test. Although applicable in theory, the analytical approach often implies excessive simplifications of the real system behavior, such as piston-like displacement. Thus only numerical simulations can thoroughly describe the phenomena occurring during the injection process. However, the pressure and pressure derivative response calculated numerically often shows nonphysical oscillations during the radial flow phase, when the pressure derivative is expected to be horizontal. It was found that these spurious oscillations arise in convection-dominated problems and are associated with sharp saturation fronts. In this paper, an effective methodology, based on an adaptive time-step calculation, is presented so as to avoid pressure oscillations. The proposed time-step selection is both computationally efficient and suitable to capture the physics of the system. Absolute permeability (m 2 ) k h

LIST OF SYMBOLS
Horizontal absolute permeability (m 2 ) k r Relative permeability (-) k z Vertical absolute permeability (m 2 ) DL Distance connecting two nodes (m) M Mobility ratio (-) n o Oil relative permeability exponent (-) n r Number of cells in the radial direction (-) n w Water relative permeability exponent (-) p Pressure (Pa) q Rate for unit volume (1/s) Q Water injection rate (Sm 3 /day) r Nodal coordinate in the radial direction (m) r 1 Radial position of first cell node (m) r e External radius (m) r w Wellbore radius (m) S Saturation (-) Dt Time-step length (s) Dt 0 Initial time-step length (s) V Volume (m 3 ) z Vertical coordinate (m)

INTRODUCTION
Well testing is a process widely used in the oil industry for evaluating the well productivity and the formation damage, estimating the reservoir characteristics such as initial pressure, fluid type and effective permeability, and identifying the reservoir heterogeneities, which are key information for field development and facilities design (Coelho et al., 2005). Well tests, usually performed during the exploration and appraisal phases of a reservoir, consist in introducing abrupt changes in the surface production rates and recording the associated changes in bottomhole pressure. The pressure disturbance induced by production travels into the formation and is affected by the rock features and fluid properties in various ways. Therefore, a record of the pressure response over time produces a curve whose shape is defined by the reservoir's unique characteristics (Schlumberger, 1998).
The selection of the test type to assess the fluid nature and the reservoir potential must be balanced against operational risk, environmental constraints and value derived from affecting early decisions on project appraisal or development. Typically, conventional testing methods involve surface production of fluids. However, in exploration and often in appraisal scenarios, surface facilities to store the reservoir fluid are not available and hence the fluid is discharged or flared. Burning hydrocarbons produces significant amounts of emissions, which in turn produce acid rain, smog, ozone at ground levels and greenhouse gases in the upper atmosphere. The demands to reduce emissions during well testing put enormous pressure to avoid these tests altogether. Alternative testing procedures have thus been investigated and attempted for reservoir appraisal so as to have sufficient information to evaluate the investment risk and make a decision whether to sanction a project or to develop the field. A valid contribution to the review and discussion of technologies such as wireline formation tests, closed chamber tests, production/reinjection tests and injection tests as viable alternatives to conventional well testing can be found in the technical literature (Coelho et al., 2005;Woie et al., 2000;El-Khazindar et al., 2002;Hollaender et al., 2002;Banerjee et al., 1998;Beretta et al., 2007;Levitan, 2002;Verga and Rocca, 2010).
Among these unconventional methodologies injection testing is considered very interesting. An injection test consists substantially in injecting a fluid, typically brine, in a potential pay zone and monitor the pressure response during the injection period and the subsequent fall-off period, in which the well is closed and the pressure tends to return to the equilibrium value. Although an injection/fall-off test is similar to a conventional drawdown/build-up test, a distinction between the two is necessary when the properties of the injected and reservoir fluids are different (Gunawan et al., 2002). In fact, the physics of injection tests is characterized by the presence of two phases in the reservoir, the hydrocarbons originally in place and the injected fluid; therefore, fluid saturations change dynamically during injection in both space and time and the permeability of the reservoir rock to each fluid will be dependent on saturations.
If convection dominates the immiscible displacement of hydrocarbons by the injected fluid (i.e. capillary pressures are negligible), the saturation profile in the porous medium is characterized by a steep transition between the flushed and the unflushed zone (Buckley and Leverett, 1941). When the process is simulated numerically, such sharp transition can induce non-physical oscillations in the pressure response (Forester, 1977). Such oscillations, called spurious in the technical literature (Le Veque, 1990), are due to the truncation error which arises from the numerical discretization of the partial derivatives (Lantz, 1971). Besides being responsible of the spurious oscillations in the pressure response, this truncation error acts as an artificial dispersion term, often denoted as numerical diffusion (Lantz, 1971), which tends to artificially flatten the sharp saturation profile. If the magnitude of the artificial dispersion term is comparable with the convection one, the description of the transition zone can be inaccurate, compromising the correct computation of the fluid mobilities at the interface and, in turn, the pressure response computed during the simulated injection test.
Several methods to reduce the truncation error were presented in the literature. Some methods act on the differentiation scheme, such as the transfer overshoot (Peaceman and Rachford, 1962), the two-points upstream approximations , the filtering techniques (Van Leer, 1977) the higher order variational approximations (Settari et al., 1977) or the addition of terms for truncation error cancellation (Laumbach, 1975). Others focus on a local grid refinement, in some cases with moving grids which follow the frontal advance (Han et al., 1987). Eventually, others adjust the time-step size in order to restrict the variation of each variable over the time-step (Jensen, 1980) or to cancel the first order term of the truncation error (Aziz and Settari, 2002). This paper presents an effective and computationally efficient criterion to select the proper time-step size in order to avoid spurious oscillations of the pressure derivative in injection testing simulation. A comprehensive numerical 2D axial-symmetric near wellbore model was implemented to reproduce the pressure response registered by a gauge positioned in the wellbore at the depth of the tested zone (Verga et al., 2011). The model simulates the flow of two immiscible phases as an implicit finite difference problem. The simulated pressure transient curves are subsequently post-processed in order to extract the contained reservoir information, i.e. mimicking a test interpretation.

NUMERICAL SIMULATION OF INJECTION TESTS
Simulating an injection test consists in reproducing the pressure response registered by a gauge positioned in the wellbore just above the tested zone. The simulated pressure signal is subsequently processed in order to extract reservoir information such as the formation permeability-thickness product (kh) and the well damage, expressed in terms of mechanical skin (S m ), which are key parameters for the well productivity estimate. The pressure increments (Dp) and the pressure derivative (Dp') are displayed on a log-log diagnostic plot. Details on the calculation of the pressure derivative are given in Appendix A.
Although the model has the extension of the typical drainage area of a well, the attention was focused on the near wellbore zone, thus the bi-phase flow model for simulation of injection tests was developed under the assumption of radial flow geometry. This hypothesis limits the application of the model to the case of vertical wells intercepting a hydrocarbon-bearing formation of constant thickness. However, exploration and appraisal wells, for which injection tests are of great interest, generally have a vertical or a slightly slanted trajectory, thus this assumption is not very restrictive. Fluid injection or production can occur along the entire thickness of the reservoir or through a limited interval in order to simulate partial penetration or perforation. The injection fluid and the reservoir dominant fluid are immiscible.
The model, described in detail in Appendix B, is made up of two convection-diffusion equations, one for each fluid phase. The equations arise from the conservation equations, assuming Darcy's flow in the porous medium. Several non-linearities occur, due to the dependency of the model coefficients on pressure and saturation.
The 2D Cartesian axial symmetric model was formulated in finite differences, where a multiplication for the block volume was introduced in order to conserve flow rates (Aziz and Settari, 2002). The simulator follows a "SSimp" scheme (Aziz and Settari, 2002). In fact, the algorithm is fully implicit in pressure and saturation and the flow equations are solved iteratively with the Simultaneous Solution (SS) method. The detailed numerical formulation is shown in Appendix C. The use of an implicit formulation with simultaneous solution in pressure and saturation avoids stability restrictions on the time-step size, i.e. numerical errors are not amplified from time-step to time-step, independently from the chosen time-step length (Aziz and Settari, 2002). Furthermore, an upstream weighting scheme was adopted for the evaluation of the saturation dependent variables, such as the relative permeabilities, at the cell interfaces. Conversely, a centered approximation was applied to the pressure dependent variables, such as the PVT properties. The midpoint weighting of the relative permeability, although most appropriate from the numerical analysis standpoint because inducing a second order truncation error, was discarded because it could give erroneous solutions. When the capillary pressure values are small or negligible, the SSimp equations result in an almost hyperbolic problem, with the true solution very close to the Buckley-Leverett solution; conversely, the numerical solution using midpoint weighting converges to a different solution, which is mathematically possible but physically incorrect (Aziz and Settari, 2002). Due to the non-linear nature of the scheme, a Newton-Raphson iterative solution was adopted to solve the set of coupled equations. Finally, an automatic time-stepping selection was adopted so that shorter intervals are used when rate changes occur and, then, the time-steps increase with a geometric progression during the flow period. This option provides a reliable calculation of the pressure derivative in a reasonable computational time. Criteria for the selection of the minimum time-step length (Dt 0 ) are discussed in Appendix C.
The chosen discretization scheme leads to a first order truncation error on a uniform grid. Since the solution of the studied PDE (Partial Differential Equation) is calculated in correspondence of a non-uniform grid, with the cell dimension increasing with a fixed geometric progression, the analytical evaluation of the corresponding truncation error is not trivial. The complexity of the evaluation is increased by the presence of diffusion and advection coefficients which vary with the independent variables (pressure and saturation). Thus an effective theoretical support was impossible to achieve.

SPURIOUS OSCILLATION ANALYSIS
The equations describing immiscible displacement in a porous medium (i.e. convection-diffusion equations) convey some numerical issues when convection predominates. In fact, when the equations are solved with differentiation schemes, non-physical oscillations can be observed near steep gradient regions (Forester, 1977).
A numerical simulation of an injection test can be affected by these spurious oscillations when capillary pressures are negligible. In fact, the saturation distribution during injection can be described by the presence of three zones (Sosa et al., 1981), as shown in Figure 1: -water zone near the wellbore, where water has completely displaced the oil and the saturation reaches the constant value 1 -S hr ; -transition zone, where the water saturation progressively decreases from the maximum value, 1 -S hr , to the irreducible value, S wi ; -undisturbed zone, where the water saturation is equal to the irreducible value, S wi . When convection is dominant (i.e. negligible capillary pressure), the saturation profile of an immiscible displacement in a porous medium is characterized by a steep transition between the flushed and the unflushed zone (Buckley and Leverett, 1941). This transition becomes steeper and steeper if the mobility ratio (M) between fluids is lower than 1 or if the relative permeability curves are described by polynomials of high order (n > 2). A moving front is observed during injection, while only little changes in the saturation values are registered during the pressure fall-off subsequent to the injection period (Levitan, 2002). Azarkish et al. (2006) compared Levitan's analytical solution for injection test (Levitan, 2002) with a numerical model, discretized with an areal 2D, unstructured Voronoy grid and solved with finite elements, where the time-step duration had an almost constant progression. The authors observed some oscillations in the late time period of the pressure derivative of the injection period. These oscillations had not been observed on the pressure derivative calculated with Levitan's model, thus their nature is clearly numerical. The authors verified that spurious oscillations tend to disappear when the mobility ratio between the injected and reservoir fluids increases.
Spurious oscillations are due to the local truncation error (Lantz, 1971), which arises from the numerical discretization of the partial derivatives appearing in the mass conservation laws describing the physical system. When a uniform time and space discretization is assumed in a numerical differentiation scheme, it is possible to analytically calculate the time-step length that improves the solution accuracy, i.e. reduces the discrepancy between the theoretical and the numerically approximated solutions (Higham, 1996). This is obtained by calculating the local discretization error S w = 1-S hr S wi < S w < 1-S hr S w = S wi r z Figure 1 Fluid distribution at the end of injection period along a reservoir cross section.
via the Taylor series and by imposing a time-step which makes the term of lowest order to cancel (Aziz and Settari, 2002). Unfortunately, the physics of the displacement process taking place during an injection test requires irregular time and space discretization. In fact, a centimetric gridding around the well is needed so as to obtain a detailed description of the solution in the near wellbore zone, while the domain affected by the pressure disturbance is almost kilometric. On such a domain, the use of a centimetric uniform grid is strongly inefficient from the computational point of view, while a grid with local refinement defined by a logarithmic progression is more appropriate (Ertekin et al., 2001). Analogous problems are encountered for time discretization. Very small time-steps (few seconds) are necessary at the beginning of the flow period in order to obtain a refined pressure derivative plot on the log-log scale, but the overall test can last hours or even a day. As a consequence, a constant time-step duration would be impractical. Thus, since the analytical approach for truncation error calculation can be applied only under simplified assumptions (i.e. uniform space and time discretization), the optimal time-step length able to reduce the local truncation error was obtained as a result of sensitivity analyses to both physical and numerical parameters. In order to do this the 2D axial symmetric Cartesian finite difference discretization introduced in Section 2 was exploited. Details of the discretization are given in Appendices B and C. An undersaturated oil-bearing homogeneous reservoir, the properties of which are summarized in Table 1, was selected as the base case for sensitivity analysis. Three reservoir fluids were considered (light, medium and heavy oil, respectively) with the PVT properties summarized in Table 2; the injected water properties are summarized in Table 3. Relative permeability curves were assumed Corey-shaped with the parameters listed in Table 4. The rate history consisted of 8 h of injection followed by 24 h of fall-off; rates, summarized in Table 5, were chosen so as to induce a pressure drop of about 30 bar in each case. It should be noted that these oscillations do not jeopardize the interpretability of the simulated pressure response since they merely affect the injection period (Fig. 2).
Furthermore, it is possible to reduce, and sometimes eliminate, the oscillation amplitude on the log-log plot making use of conventional smoothing algorithms in the pressure derivative calculation. However, an insight of these phenomena was necessary in order to verify whether the numerical errors affect the reliability of the interpretation results.
Firstly, sensitivities to numerical parameters such as grid refinement and time-step length were performed on a base case (M = 1.19, n w = n o = 2) discretized into 40 cells (a r = 1.259) with a time-step length starting from 3 seconds and growing by a factor a t = 1.01. Results confirmed that spurious oscillations are related to the local truncation error, which is decreased by a finer grid and increased by finer a time-step (Todd and Longstaff, 1972). In fact, once the time progression is fixed, a finer grid discretization (n r ) reduces the oscillation amplitude and increases the oscillation frequency (Fig. 3); in addition the frequency is not constant in time but decreases logarithmically (Fig. 3c), thus looking constant on a log-log scale (Fig. 3b). At the same time, once the grid is fixed, reduced amplitudes correspond to faster timestep progression (a t ) (Fig. 4). No significant influence of the minimum time-step length (Dt 0 ) was observed, except at the very beginning (t 0.01 h) of the test (Fig. 5).  Secondly, sensitivities on reservoir and fluid properties were performed in different scenarios, starting from a reference case (M = 0.25, n w = n o = 2, n r = 40, Dt 0 = 3 s, a t = 1.01). Sensitivities on reservoir properties such as permeability, thickness and porosity, along with the injection rate, showed a negligible impact: -the oscillation amplitude and frequency were not affected by permeability (Fig. 6a) nor pay (Fig. 6b), neither was a time shift observed at late time; -amplitude and frequency were not influenced by porosity (Fig. 6c) nor injection rate (Fig. 6d); only a time shift was observed.
On the contrary, from sensitivities on fluid mobilities a relevant impact on the oscillation behavior was noticed. It was observed that spurious oscillations due to numerical errors can appear at each mobility ratio (Fig. 7b). Oscillations are generally not visible on the pressure response, except for the case of high Corey exponents, i.e. n w = n o > 2 (Fig. 7a, 8a). Conversely, oscillations could be clearly detected from the derivative of the pressure response in all cases, except for linear relative permeability curves (Fig. 7b, 8b). Sensitivity analyses, along with the complete absence of oscillations in the fall-off derivative, showed that there is a link between spurious oscillations and water front velocity. In fact, oscillations appear later if the water velocity is reduced, either due to a low permeability or to mobility ratio (M) lower than 1. Furthermore, larger oscillations were observed in the case of slow oil velocity either due to reduced mobility in the transition zone (i.e. Corey exponent n o > 2) or to mobility ratio M > 1. The following results were found: -fixing the injection pressure drop, the lower the mobility ratio the higher the oscillation amplitude (Fig. 7); -oscillations appeared later in the case of higher mobility ratio (Fig. 7); however, the oscillation frequency did not change; -the higher the Corey exponents of the relative permeability curves, the larger the oscillation amplitude (Fig. 8b,c). A time delay was observed for low exponents (Fig. 8b), but the oscillation frequency remained exactly the same (Fig. 8c).
Since the grid has a geometric progression, the transition zone moves in time from very fine cells to larger ones. As a consequence, at fixed time-step length, the truncation error is expected to affect the  solution accuracy more severely when the transition zone reaches coarse grid cells. Thus, for the fixed time-step cases, oscillations are expected to appear at different times, depending on the velocity of the saturation front advance. In turn, the front velocity depends on the fluid mobility and thus on absolute permeability, mobility ratio and relative permeability. Figures 6a, 7b and 8b seem to confirm such considerations. Furthermore, the mobility ratio and the curvature of the relative permeability functions influence the shape of the displacing front, which is characterized by smaller saturation gradients for mobility ratio M > 1 and low Corey exponents (n 2) and by higher saturation gradients (almost piston like displacement) for mobility ratio M < 1 and high Corey exponents (n > 2) as shown in Figure 7c and 8d. For a specified time-step size, the magnitude of the truncation error increases for increasing front velocity and saturation gradient at the front (Todd and Longstaff, 1972), thus a higher truncation error is expected for Time (h) dp and dp' (bar) M < 1 and n > 2. This is confirmed by the increasing oscillation amplitude observed in such cases (Fig. 7, 8).

PREVENTING SIMULATION FROM SPURIOUS OSCILLATIONS
Spurious oscillations can be avoided either by reducing the oscillation amplitude to zero or by imposing a time-step length equal to the oscillation semi-period. As previously discussed in Section 2, the frequency, which decreases logarithmically in time, is univocally determined by the imposed grid discretization (Fig. 3). Furthermore, as the time-step progression factor increases, the oscillation amplitude tends to zero and the time-step length tends to overcome the oscillation semi-period (Fig. 4). The goal was achieved by setting a time progression factor (a t ) equal to the nodal grid progression (a r ): In this way, instead of adaptively refining the grid in the transition zone, the time-step length is gradually increased. In practice, with the adopted discretization (Appendices B and C), a time-step increase with the same progression as the grid guarantees at any time that the discretization at the transition zone is fine enough to keep the truncation error under control. 10 2 10 1 10 0 dp and dp' (bar) dp and dp' (bar)  Spurious oscillations of the pressure derivative during the injection period: sensitivity to the initial time-step length (Dt 0 = 0.1, 3 s) with a time progression a t = 1.01; mobility ratio M = 1.19 and Corey exponent n w = n o = 2. The two derivatives overlap as time increased.
As it will be empirically demonstrated further on in this section, the only assumptions are that the well trajectory is vertical and that the numerical grid follows a geometrical progression, as expressed by equations (C4)-(C7). The values assigned to fluid properties, i.e. mobility and compressibility, and to reservoir properties, such as permeability, porosity, heterogeneity, anisotropy, do not jeopardize the validity of the relation. Furthermore, the well completion and test procedure, such as partial penetration/perforation of the well or the imposed rate, are not subject to any limitation as well.
It was verified that the time-step selection according to Equation (1) successfully prevented spurious oscillations to arise in all the considered scenarios. By way of example, the effectiveness of the proposed time-step selection (a t = a r ) is shown for the case most affected by oscillations, corresponding to a mobility ratio M = 0.25 and high Corey exponents (n w = n o = 3) (Fig. 9). It was noticed that the obtained time-step length corresponded to the oscillation semi-period (Fig. 9d).
The proposed time-step selection scheme is both computationally efficient and physically reliable. In fact, on the one hand the introduced numerical Time (h) dp and dp' (bar) dp and dp' (bar) dp and dp' (bar) dp and dp' (bar) Spurious oscillations on the simulated pressure derivative for the injection period. Sensitivity analyses to: a) reservoir absolute permeability (k = 10, 100 mD) for a fixed horizontal stabilization of the pressure derivative; b) pay thickness (h = 1, 10 m); c) porosity (/ = 0.1, 0.2); d) injection rate (q w = 95, 190, 380 m 3 /day) at a fixed mobility ratio M = 0.25 and Corey exponents n w = n o = 2.
formulation ensures accuracy and stability without requiring an excessive number of time-steps or an extremely fine mesh resolution. In fact, local grid refinement and adaptive time-step length are imposed via logarithmic progressions, thus efficiently reducing the overall computational cost with respect to a uniform time discretization with the same degree of detail. Furthermore, since the length of each time-step is univocally determined before the calculations over the time-step are performed, no additional computations are required; on the opposite, in methodologies based on constraints over variable variations, the time-steps need to be recalculated if the trial time-step length proves to be too long to comply with the imposed criterion. Besides preventing oscillations to occur, the proposed discretization does not significantly flattens the physically sharp saturation front, as it can be observed comparing the saturation profiles at the end of the injection period with the corresponding analytical Buckley-Leverett curves (Fig. 9b). Therefore, a reliable description of the solution is obtained even in the transition zone, with Sensitivity analyses to mobility ratio (M = 0.25, 1.19, 10) for a fixed pressure drop and for Corey exponents n w = n o = 2: a) simulated pressure response and b) pressure derivative for the injection period; c) water saturation profiles at the end of the injection period. Spurious oscillations are not visible on the pressure profile, but they are evident on the pressure derivative. a consequent correct description of the changing fluid mobilities and, in turn, of the pressure response.
It was also verified that the progression in time (1) does not involve any stability problem because the scheme is fully implicit. Furthermore, the required time resolution of the pressure response reduces logarithmically in time. In fact, as the diagnostic plot used for interpretation is log-log and the flow periods typically last a few hours, very short time-steps are needed at the beginning of each flow period (i.e. in the order of magnitude of seconds), while long time-steps are adequate at the end of the flow period (i.e. in the order of magnitude of hours). Therefore, a geometric progression in time would provide a good log-log representation. Sensitivity analyses to the Corey exponents (n w = n o = 1, 2, 3) for the mobility ratio M = 0.25: a) simulated pressure response and b) pressure derivative for the injection period; c) zoom of the oscillations during the horizontal stabilization of the pressure derivative; d) water saturation profile at the end of the injection period. Oscillations maintain exactly the same frequency, but grow in amplitude with the Corey exponents; moreover, a shift in time can be observed. Spurious oscillations on the pressure profile are also visible for high Corey exponents (n w = n o = 3).
Relation (1) is still valuable in heterogeneous cases. Exemplifying that, two heterogeneous cases were considered: a radial composite reservoir (Fig. 10a) and a layered reservoir (Fig. 11a). For a more general validation, the fluid properties that induced the most evident oscillations were considered, i.e. mobility ratio M = 0.25 and high exponent of relative permeability curves (n = 3). The pressure derivatives of the injection period are shown in Figure 10b and 11b, respectively. In both cases oscillations were successfully removed.
The introduced relation (1) is generally suitable to eliminate spurious oscillations also when a limited entry well is considered, even though the uniform grid discretization adopted in the vertical direction is coarser than that in the radial direction in the near wellbore zone (i.e. decametric vs centimetric). In fact, as injection begins, the vertical velocity is small compared to the radial velocity and an almost radial flow occurs in front of the perforated interval. In such cases diffusion in the vertical direction can be negligible (Lantz, 1971). After a short time (a few minutes in the presented case), when spherical flow begins, vertical diffusion is no longer negligible, but a time-step progression greater than or equal to the grid factor alpha is generally sufficient. Firstly, an oil-bearing homogeneous reservoir (M = 0.25) produced by a well with limited entry (perforations corresponded to 30% of the net pay) was considered (Fig. 12a). The effectiveness of the time discretization was also verified in potentially critical cases (even though  they have limited practical interest), when the vertical permeability is higher than the horizontal one. Three cases were addressed: -anisotropy ratio 0.1 and high order relative permeability curves (n = 3), -anisotropy ratio 10 and high order relative permeability curves (n = 3), -anisotropy ratio 10 and low order relative permeability curves (n = 2). For high order Corey functions (n = 3) oscillations were successfully removed for an anisotropy ratio of 0.1 (Fig. 12b) and significantly reduced for anisotropy ratio of 10 (Fig. 12c). Difficulties related to anisotropy ratio were not encountered for Corey exponent n = 2, as shown in Figure 12d.
Afterwards, an isotropic heterogeneous reservoir with different permeability layers (Fig. 13a) was investigated (M = 0.25, n = 3), where only the fraction of pay corresponding to the higher permeability layer was perforated. Again the oscillations were successfully removed (Fig. 13b).

Time (h)
k h = 100 mD h w = 3 m h = 10 m dp and dp' (bar) dp and dp' (bar) dp and dp' (bar)

CONCLUSIONS
A numerical solution is required for an adequate simulation of the injection/fall-off tests. Although applicable in theory, the analytical approach often implies excessive simplifications of the real system behavior, thus it fails to properly describe the reservoir characteristics and the phenomena occurring during the injection process. However, it was observed that numerical solutions can be affected by spurious oscillations, particularly evident on the pressure derivative calculated for the injection period at late time. This non-physical response arises in convection-dominated problems, which are characterized by the presence of a sharp saturation front. Sensitivity analyses proved that spurious oscillations depend on the parameters affecting the injected fluid velocity.
Oscillations can be mitigated, or even eliminated, by a reduction of the local truncation error, typically achieved by highly refined grids at the cost of remarkably increasing the computational effort. In this paper, a valuable alternative, based on an adaptive time-step calculation, was presented. It was observed that a time progression factor equal to the nodal grid progression keeps the truncation error efficiently under control. In fact, spurious oscillations were successfully avoided in all the considered scenarios, given that the only required assumptions are that the well is vertical and that the spatial discretization follows a geometric progression. The proposed timestep selection scheme is both physically reliable and computationally efficient, avoiding time consuming sensitivities for the identification of a suitable time-stepping. In fact, the presented time-step length which increases logarithmically provides a nice log-log representation of the pressure derivative, free from spurious oscillations also in convective-dominated problems. Furthermore, the computation of the saturation profile is not significantly affected by numerical errors, such as an artificial flattening of the front corresponding to an overestimation of the transition zone extension. Finally, accuracy and stability of the solution are ensured without requiring an excessive time or mesh resolution. Adaptive time-step length and local grid refinement are imposed via logarithmic progressions, significantly reducing the overall computational cost with respect to a uniform discretization with the same degree of detail. Opposite to methodologies based on fixed maximum variable variations, no additional computation is required because the length of each time-step is univocally determined before the calculations over the time-step are performed.
To some extent the proposed solution might appear pragmatic because it was derived empirically and it has a very simple form: it avoids oscillation by calculating the model solution with a time-step length equal to the oscillation semi-period. Nevertheless, it was shown that its validity is quite general in the context of injection test application. In fact, it was successfully applied to a variety of scenarios, from light oils to heavy oils, from fully perforated to limited entry wells, from homogeneous to heterogeneous cases (either radial composite or layered reservoirs). In addition, since injection rates and formation properties (such as permeability, thickness and porosity) do not significantly influence the oscillation phenomena, the applicability of the methodology is not limited to a range of parameters. For these reasons, the simplicity of the proposed methodology should be valued and not regarded as a limiting aspect.
where water pressure (p w ) and water saturation (S w ) were chosen as the unknowns. If the capillary pressures are zero,