Numerical Modeling of Periodic Pumping Tests in Wells Penetrating a Heterogeneous Aquifer

This study investigated how to utilize multiple frequency components of pressure data from periodic pulse tests to estimate the intra-well permeability and compressibility distribution and also the presence of heterogeneities in a real field case. Periodic well testing is a technique in which injection or production pulses of a fluid are applied to a well in a periodic fashion. One of its main advantages is that ongoing operations do not have to be interrupted during the test as the superposed harmonic components can be identified using Fourier analysis. Further, modeling calculations are much faster than calculations in the time domain as no time-stepping is required and only the frequencies observed in the test need to be evaluated. We applied an earlier developed numerical code in the frequency domain to evaluate periodic-test results in a shallow aquifer and obtained a good match between data and calculations. The interpreted formation heterogeneity is in line with the local geology. Joints of various orientations constitute the main hydraulic conduits of the tested subsurface but they do not directly connect the wells. Thus communication between the wells has to be established through low-permeability features. The interwell periodic testing has corroborated the geological understanding of the aquifer and helped understanding the fluid flow pattern.


INTRODUCTION
Well testing is a valuable source of information on reservoir properties and a lot of work has been done and is being performed in this realm (Kruseman and Ridder, 1990;Horne, 1995;Fetter, 2001;Bourdet, 2002;Gringarten, 2008;Fokker et al., 2005;Beretta et al., 2007;Benetatos and Viberti, 2010;Viberti, 2010;Verga et al., 2011;Cancelliere and Verga, 2012). Two issues remain which need further attention in well test analysis: the inclusion of reservoir heterogeneity and the treatment of rate variations. Reservoir heterogeneity which cannot be handled with traditional analytical methods is usually approached using numerical techniques Copty and Findikakis, 2004;Ye et al., 2004;Houze et al., 2010); rate variations are approached with deconvolution techniques (Schroeter et al., 2004;Gringarten, 2010).
Periodic pulse testing in the way we presently apply it addresses these two issues. Periodic pulse testing is a form of testing in which the flow rate is varied a number of times with a fixed frequency and amplitude. Pulse testing was first proposed by Johnson et al. (1966); other authors have later elaborated on the idea and particularly on the benefits of periodic pulsing (Kuo, 1972;Black and Kipp, 1981;Marschall and Barczewski, 1989;Rosa and Horne, 1997;Hollaender et al., 2002;Renner and Messar, 2006;Rochon et al., 2008;Ahn and Horne, 2010;Fokker and Verga, 2011;Fokker et al., 2012). The periodic pulses consist in the repetition of two different rates at regular intervals. Pulses are created in one of the Science Publications AJES wells, called active well or pulser well and the pressure response can be measured both in the pulser well and in observation wells, also called responders. Because pulse tests do not require the observation wells to be shut-in, only limited disturbance is created to production or injection operations. The duration of the pulse and the flow rates are estimated a priori on the basis of the expected reservoir characteristics to ensure that a pressure variation can be detected at the responders. In all cases the amplitudes of the pulse test response measured at the observation wells are small, frequently less than 1 bar and sometimes less than 0.1 bar. The test duration may be longer than for a simple interference test because several periods are required for the frequency analysis, but the advantage is that the oscillating response due to repeating pulses is easier to identify in a noisy reservoir environment than a single pulse interference signal and that it is less affected by a possible drift of the pressure gauge.
The analysis of the distinct frequencies in the pressure response requires a harmonic approach to the pulse test interpretation (Renner and Messar, 2006;Ahn and Horne, 2010;Fokker and Verga, 2011). Fourier transformation unravels the signal by picking out the imposed frequencies and the corresponding response function. The identified amplitude and phase responses provide information about the reservoir properties in the volume of investigation. The amplitude response represents the signal strength in relation to the imposed signal; the phase shift represents its relative delay. In contrast to conventional well testing analysis, neither pressure derivative analysis nor deconvolution is required.
The interpretation of periodic pulse tests, based on the analytical solution of flow equations in the frequency domain, can be successfully applied to multiphase flow in anisotropic formations (Fokker and Verga, 2011). Ahn and Horne (2010) developed a solution for periodic testing that accounts for radial multi-composite reservoirs. To be able to obtain reliable results in more complicated heterogeneous systems we recently developed a numerical simulator in the Fourier space (Fokker et al., 2012). In this study we report on the application of this simulator to results of a field study (Renner and Messar, 2006). The original interpretation using analytical expressions for a homogeneous reservoir yielded highly varying results for the effective transmissibility and the effective storativity for the range of imposed pumping frequencies. This variability motivated us to use the same data in the present study in order to explore to what extent reservoir heterogeneity can explain all the measurements simultaneously and to what extent the characteristics of the heterogeneity can be resolved from such data sets.

Numerical solution in the Fourier space
A comprehensive description of the numerical solution scheme is provided in Fokker et al. (2012). Here we only give the main line of thought and the extension presently applied. We consider a heterogeneous reservoir with space-dependent permeability k containing a slightly compressible fluid. The fluid pressure p obeys a diffusion Equation (1): where, Φ, c and µ denote porosity, fluid compressibility and fluid viscosity, respectively. Boundary conditions at the reservoir edge are a constant pressure for an open reservoir and no flow for a closed reservoir.
For water, compressibility and viscosity are nearly independent of pressure in the range of hydrological applications and thus the diffusivity equation is linear. The linearity allows for treating the combined pressure solution of many frequency components as a superposition of the individual contributions and each harmonic component p ω (x, t) = P ω (x)e iωt can be considered independently. We evaluated the time derivative and incorporated the source term q ω (t) = Q ω e iωt to obtain an equation in the frequency domain without explicit time dependence. In addition to our earlier work we also incorporated a wellbore storage coefficient C w = c f V well in the equation for the cell containing the pulser. We then have: Equation (2) can be discretized by the finite difference method so that a numerical solution in the spatial domain can be obtained independently for each frequency. We have used a 2D implementation, but the approach is fully applicable to 3D.
Simulation results can be compared with real data in terms of amplitude ratio and phase shift. The amplitude ratio A ω is defined as the ratio between the pressure amplitude and the rate amplitude for the corresponding frequencies; the phase shift ϕ ω is defined as the normalized time delay between the pressure and the rate for the corresponding frequencies (Equation 3):

Aquifer Test Details
We re-evaluate previously presented experiments in three shallow boreholes penetrating a jointed sandstone formation (Renner and Messar, 2006). Pumping tests had been performed in a triangle of three water wells, called BK1, BK3 and BKs at distances of 20.6, 38.1 and 51.8 m. The western and southern sides of the reservoir are bordered by a quarry wall; the south-eastern side is bounded by a fresh-water reservoir. The boreholes penetrate steeply-dipping layers of sandstone and claystone that cannot be correlated between the wells (Fig. 1). A series of injectivity and interference tests for puckered intervals demonstrates that joints constitute the dominant hydraulic conduits in the subsurface. While the wells hydraulically communicate the communication apparently has a strong spatial restriction. For example, communication between wells BK1 and BK3 was shown to be strongest between the bottom of BK3 and a shallow Science Publications AJES zone in BK1 in some agreement with the traces of joints projected from borehole logging information (Fig. 1).
Periodic pulses had been applied both in well BK1 and in well BK3 and the response was measured in the pulser and in the observer wells. Cycle periods varied between 10 s and 5400 s providing the opportunity to test a large range of frequency responses both in the pulser and in the observer wells. The imposed injection and production rates had been corrected for wellbore storage effects.

RESULTS
We re-analyzed the original pressure tracks. A fast Fourier transform was first operated on the corrected flow rate and on the detrended pressure records (example given in Fig. 2) yielding amplitude spectra (Fig. 3). In the rate spectrum of every experiment, only those frequencies were selected that gave at least three times larger amplitude response than the background. For these frequencies the pressure responses were evaluated: only if the amplitude was significantly larger than the background signal and if the phase difference with the corresponding signal in the rate was in line with the general trend they were used in the analysis. In this way, all frequencies for which a good response could be identified rather than restricting the analysis to the leading frequencies were used in the analysis. Indeed, as it was already pointed out, the fact that the records contain more frequency components than the leading ones, is one of the advantages of periodic testing with square or non-sinusoidal pulses (Fokker and Verga, 2011). The signals used in the further analysis were the amplitude ratio and the phase difference between pressure and rate signal responses in the Fourier domain.
The identified frequency spectra were combined to a single set for every test setup. Five setups delivered useful results: the observed pressures in BK1, BK3 and BKs for pulses applied in BK1 and the observed pressures in BK3 and BK1 for pulses applied in BK3. The amplitude ratio and phase difference values are considered to be independent of the test details like injection rate, but they are influenced by the local permeability distribution. As a first attempt an homogenous permeability field was assumed; matching results of the tests where BK1 is the pulser are shown in Fig. 4 and 5 in terms of permeability and total compressibility, respectively.
At a later stage we used our 2D numerical tool to arrive at a single heterogeneous permeability field matching all the data. For the reservoir thickness we took 5 m, which was deemed a representative number for the high-permeability streaks identified with the well logs and we chose constantpressure boundary conditions because of the reservoir characteristics. The matches obtained from each single-well response suggested the presence of two homogeneous zones of different permeability.
The presence of these zones is in line with the characteristics of the system as observed from the well logs, showing uncorrelated large-permeability streaks. We completed the geological setting by introducing permeability boundaries between the wells, which ensured the reproduction of the response in the observer wells. Pressure is thought to communicate from well to well through lower-permeability layers of considerable extent. A map of the different zones is presented in Fig. 6. We matched the test responses using a geometry with four simplified geological zones. The matching procedure was quite ill-conditioned and largely depended on the weighting of the different sources of information. Most areas mainly influenced one of the 5 test scenarios; only the BK1-BK3 barrier critically affects both the pulser-observer pairs BK1-BK3 and BK3-BK1. To cope with the different dependencies an iterative matching procedure was employed starting with the most sensitive parameters and then fine-tuning the model response by adjusting all the parameters together. After the general geometry had been determined using a broad sensitivity analysis we subsequently applied an optimization algorithm to estimate permeabilities in different zones: firstly around the pulser wells using the responses in those; then between the pulser and observer wells using the observer responses; and, finally, accounting for all the measurements simultaneously. Results are reported in Table 1.
In addition to the reservoir properties, we accounted for wellbore storage effects. In fact, the wellbore acts as a storage tank when the rates are changed and causes the flow rate at the interface between well and reservoir to approach the requested surface flow rate as a function of time (Ramey, 1965). Results are shown in Fig. 7 and 8.

DISCUSSION
By assuming a homogeneous field permeability the responses in the pulser could be matched when wellbore storage and skin were employed and when constantpressure boundaries were used. However, as expected from the earlier work by Renner and Messar (2006), it was impossible to simultaneously match the responses in wells BK1, BK3 and BKs upon pulsing in BK1 with a homogeneous model, as explicitly demonstrated in Fig. 4 and Fig. 5. Adjusting the parameters (i.e., k and c f ) to honor all the recorded responses yielded poor results. We also evaluated an analytical model and obtained the same unsatisfactory results. The model was a homogeneous, radially-symmetric scheme with finite-extent and open or closed boundaries and contained two-layers: the lowerpermeability layer fed a higher-permeability layer which connected to the wellbores. These analyses corroborated that a heterogeneous model was needed.
Considering the heterogeneous permeability configuration with four simplified zones (Fig. 6) and using an optimization procedure the test responses were successfully matched. Such a process can be performed by using any standard optimizer. In this study the active set method for constrained programming was used within a Sequential Quadratic Programming algorithm (Gills et al., 1981). This iterative method calculates a numerical approximation of the Hessian of the Lagrangian function to generate a Quadratic Programming subproblem whose solution is used to form a search direction for a line search procedure. The active set strategy was adopted to constrain the optimization parameters to physical ranges.
The optimized permeabilities in the barriers are much smaller than in the reservoir parts connected to the wells ( Table 1). Compressibilities did not vary as consistently and showed smaller differences. It should be noted here that the resulting values for both are correlated with the layer thickness of the different zone, which was chosen to be 5 m for all. The wellbore storage numbers amounted to C w = 0.1 m 3 /bar for BK1 and C w = 0.03 m 3 /bar for well BK3. A good final match was obtained between the measured and calculated responses from both the pumping in BK1 and in BK3 ( Fig. 7 and 8, respectively).

CONCLUSION
We applied a numerical simulator in the frequency domain to a field case of periodic pulse testing. The simulator performs the calculations for a limited number of distinct frequencies rather than for proceeding time steps. As a result, calculations can be performed much faster. For periodic pulse testing, field operations can also be simpler than for conventional well testing because production does not have to be interrupted but can be continued during testing. In fact, the frequency signal can be isolated from the overall signal using Fourier analysis.
Our field example is a well triplet penetrating a shallow aquifer. Imposed pulsing periods ranged from 10 to 5400 seconds. Numerical interpretation of the periodic tests facilitated a realistic and simultaneous match of all the frequency components of the analyzed well pairs, in contrast to the earlier analysis with analytical modeling. An effective 2D description of the aquifer geological features was obtained in agreement with the well log data. The zones of decreased permeability between the wells are in line with the geological setting of a very laminated reservoir with steeply dipping layers. As a result the layers are not continuous between the three wells and communication has to take place through the intermediate layers. The interwell periodic testing thus has corroborated the geological understanding of the reservoir.
Prior geological and operational information was required to constrain the ill-posedness of the inverse problem. In particular the geometry of the various regions was input in our optimization scheme. A further improvement can presumably be obtained by employing Science Publications AJES a global optimization scheme which also updates the geometrical layout. Ensemble-based schemes like the Ensemble Kalman Filter (Evenson, 2009) are promising in this respect as they combine a robust optimization with a probabilistic approach and result.