## Abstract

This paper describes the new concept of a ‘Statistical Reservoir Model’ to determine significant well-pair correlations. We solve this conceptual problem using a predictive error filter, combined with Bayesian methods that identify those well pairs that are related to each other with statistical significance, for the Gullfaks reservoir in the North Sea. Significant, long-range, correlations in the whole field are found at an optimal time lag of one month. The correlation function for significantly-correlated well pairs, after normalization for the distribution of available wells, shows a long-range power-law decay that is consistent with a critical-point response at the reservoir scale. A principal component analysis shows a strong correlation with the location and orientation of faults that intersect the main producing horizon. A predictive experiment shows that the model performs very well both in history matching and predictive mode on a time scale of about one month.

Hydrocarbon reservoirs comprise subsurface bodies of rock of suitable porosity and permeability to allow the storage and transmittal of oil and gas. Producer wells are sunk into a reservoir to allow hydrocarbons to be extracted, and injector wells to provide voidage replacement and/or to maintain a positive fluid overpressure that will facilitate an efficient sweep and extraction. Pressure changes associated with both injection and production in a poro-elastic material allow geomechanics to have a significant influence on hydrocarbon production rates through changes in the effective stress field (Segall 1989; Main *et al.* 1994; Maillot *et al.* 1999; Ngwenya *et al.* 2003; Zimmermann & Main 2004). Geomechanics not only predicts a reservoir response in the near field, but also at long range i.e. distances much greater than would be predicted by conventional drainage models (Heffer *et al.* 1995; Maillot *et al.* 1999). Direct evidence for the relevance of far-field stress and strain changes induced by reservoir pressure changes and poro-elasticity include surface subsidence and induced micro-seismicity (Healy *et al.* 1968; Segall 1989; Segall & Fitzgerald 1998; Grasso & Sornette 1998; Marsan & Bean 1999, 2003; Zoback & Zinke 2002; Rutledge *et al.* 1998, 2004). Changes in the far-field strain, and hence pressure, may occur through compliant faults and fractures (Segall 1989; Maillot *et al.* 1999; Rutledge *et al.* 2004; Zimmermann & Main 2004) or through matrix material with a stress-sensitive porosity or permeability (Main *et al.* 1994; Ngwenya *et al.* 2003). All these mechanisms of pressure, strain and permeability change can have a direct influence on production rates at long range. Such changes are particularly significant when the systems of faults and fractures and the stress field acting on them are in, or close to, a state of incipient failure or criticality (Main 1996). In such a state the hydraulic properties of the reservoir are at their most sensitive to changes in effective stress (Main *et al.* 1994; Ngwenya *et al.* 2003; Zimmermann & Main 2004).

## Long-range correlations in reservoir behaviour

One of the hallmarks of a critical point system is the presence of long-range correlations. For example natural earthquake sequences exhibit long-range triggering with a correlation length of 10 km to 20 km, and a maximum triggering distance of 150 km (Huc & Main 2003). Interestingly, oilfield well-rate data had also shown correlations at lag-time of less than a month over distances up to about 20 km, distances far too large to be due to Darcy flow between the relevant injector and producer wells in the time available (Heffer *et al.* 1995; Papasouliotis 2000; Huc & Main 2001).

The existence of long-range correlations is one more independent piece of evidence that the Earth's crust is in a near-critical state. Another characteristic property of critical systems is the sensitivity to what might otherwise be thought of as very small fluctuations. For example, 0.1 MPa is sufficient to trigger earthquakes (Grasso & Sornette 1998; Stein 1999). Other evidence includes the very low stress perturbations required to trigger micro-earthquakes in the far field (Healy *et al.* 1968; Segall 1989; Shapiro *et al.* 1997, 1999; Grasso & Sornette 1998; Segall & Fitzgerald 1998; Maillot *et al.* 1999; Zoback & Zinke 2002; Rothert & Shapiro 2003; Rothert *et al.* 2003), micro-seismicity along pre-existing, critically-stressed natural fractures or zones of weakness (Shapiro *et al.* 1997, 1999), preferential fluid flow from critically-stressed fractures in borehole data (Rothert & Shapiro 2003; Rothert *et al.* 2003), and micro-seismicity related to important faults (Simiyu 1999). Added to these evidences for criticality, at which all scales contribute to collective behaviour, is the scale-free population dynamics of faults, fractures and natural seismicity (Main 1996). These all testify to the ability of the Earth, through tectonic forcing and threshold dynamics, to tune itself to a global state of ‘self-organized criticality’ (Bak & Tang 1989) with long-range correlations.

Correlation techniques from well rate fluctuation at pairs of wells have previously been used in oilfield analysis to demonstrate a relationship with local stress state and therefore the involvement of geomechanics in reservoir flow (Heffer *et al.* 1995, 1999; Heffer 2002). In contrast other workers have interpreted inter-well connectivity observed with this technique in terms of more local Darcian flow (Jansen & Kelkar 1996, 1997*a*, *b*; Refunjol & Lake 1997; Soeriawinata & Kelkar 1999; Albertoni & Lake 2002). In these studies, monthly well rates at each well are treated as time series, and then the Spearman rank correlation coefficient is calculated between rate series at well pairs. The rate correlations are considered to be indicators of inter-well communication through the reservoir, and therefore provide a potential tool for planning and managing waterfloods. The orientational trend found by Heffer *et al.* (1995), showing a strong alignment of the direction of the correlation in stacked data with the local maximum horizontal stress axis in all eight field cases tested, is similar to that of injection-induced micro-seismicity in a gas field (Rutledge *et al.* 2004). The main advantage of the Spearman rank correlation method over traditional least-squares regression or the Pearson correlation coefficient is that it does not rely on assuming an underlying frequency distribution to the rate fluctuations and is a definitive indicator of a true correlation, assessable at a chosen significance level. Its main disadvantage is that, in dealing with ranks, the phase information is lost, and hence it cannot be used for direct prediction of future rates.

This paper describes a new concept in the analysis of well rate data in the form of a truly predictive ‘Statistical Reservoir Model’. We demonstrate that the model may be used for calibrating faults and fractures, and to predict reservoir response to water flood through examining the flow rate statistics of injectors and producers. The concept behind the statistical reservoir model is simple. Given only a set of input flow rates at injectors and producers in the past, how can we predict the output flow rates at producers in the future? The model produces a matrix or array of correlation coefficients that have been identified as statistically significant at a given time lag. For a single time lag the Statistical Reservoir Model is a matrix that acts like a two-dimensional ‘transfer function’, for example converting an input set of production and injection data for a given month into a predicted output the next. For the reservoir examined here we used monthly flow rates for 106 wells over 11 years. The computational algorithm takes approximately 20 minutes for each producer well, implying that the full statistical reservoir model can be determined overnight. Larger oil fields with greater sampling rates or longer lifetimes would require more time for the computation.

Previous attempts to construct such statistical models have failed due to contamination of the output by chance correlations or amplified noise, but we have solved this problem using state-of-the-art Bayesian methods (Papasouliotis 2000). Here we test our new method on a field in the North Sea where there are sufficient data to establish the model, and validate the model through examining the principal components of the correlation matrix and by predicting the ‘future’ response using data not used to perform the history match. The first suite of results comparing the output of the model and its interpretation, in comparison with earthquake–earthquake correlations in natural seismicity, is presented separately (Main *et al.* 2006). Here we present more detail in the analysis of the results, in particular analysing the degree of predictability at different time lags. We also compare the results on real data with the output of a generic coupled flow– geomechanical model (Zhang *et al.* 2007) based on the poro-elastic mechanism operating in a near-critically-stressed crust, confirming such a geomechanical effect as a plausible mechanism for the observations.

## The Statistical Reservoir Model

The concept of a Statistical Reservoir Model is illustrated in Figure 1. The input data is an array of well pressure or flow rate signals, measured in the past up to the present, at different locations in the reservoirs. The output data is the pressure/rate response at a given time in the future. The array *C*_{ij} is sparse because chance correlations are deliberately screened out by our new technique.

The Statistical Reservoir Model itself is illustrated in Figure 2. Each element in the array is the product of a real regression slope *W*_{ij} and a binary filter *N*_{ij} to retain only those correlations that are significant (*N*_{ij}=1). The result is a parsimonious array representing the response of the *j*'th producer to the *i*'th injector or producer at different time lags *k*. To determine the optimal regression model we minimize the prediction error with respect to the model parameters, resulting in a standard regression model (Draper & Smith 1998). Bayesian methods (West & Harrison 1997; Leonard & Hsu 1999; Papasouliotis 2000) are then used to determine the significance matrix. The method of solution and the Bayesian techniques used are outlined in the Appendix and presented in more detail in UK patent application 0524134.4.

To calibrate the response of faults and fractures to injection and production, the significance matrix *N*_{ij}, for time lags of 0 or 1 month (Fig. 2) is decomposed to obtain its principal components, each being a set of values across well locations. The first few principal components encapsulate the major proportion of the variance in pressure/rate fluctuations. Heffer *et al.* (1999) and Heffer & King (2006) describe one method of interpolating values between wells for individual components to emphasize spatial trends which are consistent with geomechanical principles. Finally we examine the performance of the model in both history-matching and predictive mode.

Although the Statistical Reservoir Model is capable of predicting pressures/rates, its validity is dependent upon constancy of the physical processes in the reservoir; any changes to process can only be handled with a physically-based simulator. A combination of the two approaches is therefore the most advantageous in any practical situation.

## Results

The concept of the Statistical Reservoir Model presented above has been applied to flow-rate data from the Gullfaks field in the North Sea. The main Gullfaks field lies in Block 34/10 in the northern part of the Norwegian North Sea, and has been developed from three production platforms. The Gullfaks A platform started production on 22 December 1986, Gullfaks B platform on 29 February 1988 and the C platform on 4 November 1989. The recorded production data are composed of monthly measurements of flow rate (volume in standard cubic metres, Sm^{3}) per day, averaged over 1 month intervals) taken over a period of 133 months starting from December 1986. A total of 106 platform wells, 27 injectors and 79 producers were recorded in the dataset during the time period December 1986 to December 1997. This field has been under water flood since January 1988.

In initial examination of the raw flow rate data (treated as time series), we observed that there were some periods of zero flow rate in individual wells (e.g. see Fig. 3) due to well downtime. Further, a number of wells were operated as producers for some time prior to conversion to injectors (injection of water and/or gas). Such wells are possible sources of error or bias in examining the correlations using the Spearman rank technique, since they could provide spurious correlations and they have been treated separately by Heffer *et al.* (1995) and Albertoni & Lake (2002). However, such arbitrary treatment could lead to further systematic biases in determining inter-well correlations. Our new method for the Bayesian analysis includes a phase of pre-processing to handle missing data. In the regression stage, the missing values are omitted. However, in the Bayesian analysis stage, the missing values are tackled in two ways. For predictor variables, each missing value is replaced by a weighted average from neighbours. For response variables, a missing value is imputed from a Gaussian distribution, conditional on the observed values of predictor variables (Little & Rubin 1987).

The flow rate at wells produce time series that are both ‘noisy’ and ‘patchy’ compared to other applications of time series analysis (Fig. 3). Injection data include many large and sudden systematic fluctuations, but fluctuations in production data seem to be less discontinuous. This illustrates the ‘damping’ effect of the hydraulic, mechanical and structural responses of a complex reservoir.

Only those wells in operation for more than 24 months were selected for this test. For these data the optimal regression model (*W*_{ij}) for each producer was first established via a history match of past data, and then the Bayesian methods were used to identify those wells that are significantly related to each other (*N*_{ij}). Finally the Statistical Reservoir Models (*C*_{ij}) were formed for periods of 0–67, 0–77, 0–85, 0–103, 0–115 and 0–133 months via the product of a real regression array and a significance matrix. In establishing the Gullfaks models, an optimal lag time of one-month was found for the most significant correlations. For example, the effect of time lag between pairs of wells on the cross-correlation and autocorrelation methods is shown in Figure 4. The one-month lag correlations have more significance than the zero-lag correlations, implying a causal relationship, but also a synchronicity in response at zero lag. The optimal time lag (one-month) and scale lengths (several km) involved confirm that the correlations are not solely caused by Darcy flow.

Examining the correlation matrix *C*_{ij}, we observe that each producer has only a few wells, between about 7 and 25, that are significantly correlated to it, confirming the parsimonious nature of *C*_{ij}. As an example, Figure 5 shows the spatio-temporal distribution of significant correlated well pairs, (i.e. for wells *j*, where *N*_{1j}=1, rather than 0, in the significance matrix of Fig. 2) that are correlated to the producer well 1 over a period of 85 months. Long-range correlations are observed up to a distance of about 4.5 km and are consistent over time, in line with the results of Heffer *et al.* (1995), Papasouliotis (2000) and Huc & Main (2001). Such long-range correlations are similar to the correlation length-scales in hydrocarbon production-induced seismicity (Healy *et al.* 1968; Segall 1989; Grasso & Sornette 1998; Segall & Fitzgerald 1998; Rutledge *et al.* 1998, 2004; Marsan & Bean 1999, 2003; Zoback & Zinke 2002) and in the natural earthquake event correlations (Huc & Main 2003). All of these are likely to be due to geomechanical feedbacks.

The correlated pairs of wells are not uniformly distributed in the reservoir, and, in this case, take up (broadly) an approximately north–south alignment. This alignment is parallel to the orientation of the major faults in this field (see the fault map of the Gullfaks Field, Fossen & Hesthammer 1998), similar to the alignment of hydrocarbon production-induced seismicity with faults (Simiyu 1999; Rutledge *et al.* 2004), both implying a significant mechanical response. The most likely explanation of these effects is that the Gullfaks field is in, or close to, a state of criticality, such that the poro-elastic stress disturbances caused by fluctuations in flow rates induce far-field poro-elastic stress changes via large-scale inelastic deformation involving shear, possibly along faults (Main *et al.* 2006; Zhang *et al.* 2007). We shall see that this relationship with faults is complementary to a relationship with stress state.

The evolution of the correlated well pairs with time for this oilfield (Fig. 6) is also consistent with that observed in induced micro-seismicity elsewhere. Figure 6 shows the evolution of the raw number of significant correlated well pairs at different distance away for six different durations for the Gullfaks oilfield data. We note that:

most of the correlated wells occur at distances between 1 km and 6 km;

a peak in the number of correlated wells is reached at around 2 km;

the height of the peak decrease and the curve becomes flatter with time; and

no effect is detected for distances large than 10 km at this example, i.e. where no wells are available.

These results are of course strongly conditioned by the prior distribution of well location, but when corrected for the background distribution of well positions, result in a power-law correlation function and anomalously slow diffusion of the form 〈*x*〉 ∼ *t*^{H}, where <*x*> is the mean correlation length, *t* is the duration of the record, and the exponent *H*=0.33 (Main *et al.* 2006). The slow diffusion observed has a direct analogue in the triggering of natural seismicity (Marsan & Bean 1999, 2003; Huc & Main 2003). This result confirms that the long-range correlation is not caused by a solitary wave, which would have *H* =1 or Fickian pressure diffusion (*H*=0.5; Main *et al.* 2006).

It is well known that the power-law behaviour of earthquakes frequency distribution has been associated with self-organized criticality (Main 1996; Grasso & Sornette 1998; Turcotte 1999). Main *et al.* (2006) confirmed that the correlation function of Figure 6, corrected for the well locations, is also a power‐law for the period 0–85 months in the Gullfaks data. Figure 7 shows the frequency of correlations in the histogram of Figure 6 as a function of distance for a variety of other time periods, this time normalized for the pre-existing distribution of well locations. The fall-off of correlations at the three additional durations confirms the results of Main *et al.* (2006), that the best fitting distribution is the power-law, as seen in earthquakes, therefore consistent with self-organized critical behaviour (Main 1996; Grasso & Sornette 1998; Turcotte 1999; Huc & Main 2003). The long-range nature (Fig. 7) observed is also in good agreement with the slowness of the fall-off in a generic coupled geomechanics–flow model but only when the pre-existing tectonic stress state is near-critical (see cases 1 and 2 in the accompanying paper of Zhang *et al.* 2007).

To calibrate the response of faults and fractures to injection and production, the significance matrix *N*_{ij} (in Fig. 2) from the established Gullfaks Statistical Reservoir Model was decomposed into its principal components (or eigenvectors) for each of the available cumulative times of analysis: 0–67; 0–85; 0–103; 0–115; 0–133 months. Each principal component had a value at each of the wells in the field. Interpolation between wells was performed using the strain interpretation technique of Heffer *et al.* (1999) and Heffer & King (2006). In their studies, it is assumed that the rate fluctuations are direct indicators of geomechanical deformation. In particular each value is assumed to be equivalent to a function of the strain tensor at the well concerned. The interpolated strain maps of the first principal component (presented in Main *et al.* 2006) indicate strong linear features associated with well rate correlations in this field.

The mapped first component for the early stages of field development (0–67 months), overlain on the map of faults at top Etive Formation (Fig. 8) suggests association of two of those linearities with two major north–south faults that bound a rotated fault block: supplementary NE–SW trending linearities are seen at later times (Main *et al.* 2006). Linearities are less clear on the interpolated strain map of the eigenvector corresponding to the second highest eigenvalue in the earlier months of development (0–67 months, Fig. 9). However, a strong north–south linearity does develop by 103 months that persists until 115 months; this linearity follows the trace of the major north–south fault that provides the boundary between the ‘domino’ and ‘horst’ regions of the structure (Fig. 9). After 133 months the emphasis of the second eigenvector has shifted further to the east. This is presumably associated with later development of the Cook and Statfjord reservoirs to the east of the initial development of the Brent reservoir. Our interpretation of these principal components is that they reveal hydraulically reactive features that align in position and orientation with the predominant faults in the area revealed by seismic data (Main *et al.* 2006). The correlation of aspects of the Statistical Reservoir Model with structural properties is important because the model was not conditioned a priori by this independent data. This therefore represents an independent validation of the method.

In summary, the application of the Statistical Reservoir Model to the Gullfaks field has yielded inter-well dependencies that are (a) long-range up to the scale of the field domain and (b) fault-related, with principal components of the correlation matrix that have trends along the major north–south faults. These characteristics imply that the mechanism behind the inter-well signals is not hydraulic fracturing. Although hydraulic fractures, facilitated by thermally-induced stress reductions, are quite likely to have been created as local features at each of the injector wells, it is highly unlikely that these would have caused connections between wells at distances up to field scale for two reasons: (i) injection pressures rapidly dissipate and fall below the magnitude of the minimum principal stress, which itself increases beyond the region of convective cooling from injection; (ii) long continuous hydraulic fractures would be obvious to the field operator through pressure transient tests or extremely early water breakthroughs. Characteristic (b) adds weight to the alternative explanation of shearing on or around pre-existing faults or fractures caused mainly by the poro-elastic stress changes accompanying rate and pressure changes. Regions of dilatation and compression result from shearing which are likely to be associated with increases and decreases in local permeability respectively. Coherent shearing on a complex of faults under the over-riding influence of the *in-situ* stress state, whether compressive or extensional, is also consistent with characteristic (a). Additionally the north–south ‘domino-style’ faults in the western part of the field, which are shallowly dipping at about 30° to the east, are also quite likely to be activated in sections, possibly contemporaneously with strike-slip shearing on other faults. It must be added that the interpolation technique outlined above will preferentially tend to link high values of a principal component that are spatially close. Therefore the north–south trends along faults in Figures 8 and 9, although revealing true coherences in rate fluctuation of some degree, are not necessarily reflecting the highest correlations in rate. A better way of mapping multiple-well correlations is to be sought. The relationship (b) of principal components with faults, if taken in isolation, might be explained in terms of conventional Darcy flow, focused by the faults either as barriers or as natural conduits, with no geomechanical changes involved, but this is inconsistent with (a).

The Statistical Reservoir Model has also been validated independently in predictive mode, as illustrated in Figures 10 and 11. In this experiment the data was split into two sections labelled ‘past’ and ‘present’. The past history was first used to construct a predictive statistical reservoir model. The model was then left unchanged, and used recursively to predict the future response for the later part of the time series (which importantly was not used to construct the model), one month at a time. Apart from a few outlier ‘spikes’, the prediction of flow rate (month by month) is successful within 95% confidence limits that are determined not by the future data, but from calibration on past data (Fig. 10). For the aggregated rate history of a group of 28 wells, both the history match and the predicted rates are in good agreement with the observed flow rates, as shown in the summed time series of Figure 11. Some outliers survive despite the aggregation, but the prediction of flow rate again performs very well both in history matching and predictive mode, indicating a statistical convergence through averaging. This confirms that the statistical reservoir model can in principle be used to predict the response of the reservoir to given injection scenarios at timescales of one month or so. In future tests will be necessary to validate the model in truly predictive mode using true ‘future’ data.

## Discussion

The results presented here show long-range correlations in flow rate that are consistent with a critical point response to stress perturbations induced by pressure change. It is difficult to suggest an alternative explanation for a power-law correlation function up to a distance of 10 km at up to one month lag. The similarity of the characteristics of correlations from Gullfaks with those from several other fields previously analysed by Heffer *et al.* (1995), and from generic coupled geomechanics and flow modelling by Zhang *et al.* (2007) both lend weight to the geomechanical explanation. However, the results at shorter range may in principle reflect a more conventional reservoir response based on geological architecture. Similarly the principal component analysis shows a pattern that follows the pre-existing normal faults revealed by seismic data which might or might not be geomechanically active in the present-day stress field. A principal component analysis of itself does not reveal the nature of what the principal component means. Here it may simply be reflecting the long-term sealing or transmissibility characteristics that may be modelled by conventional reservoir description. Alternatively there may be a component of incipient reactivation, consistent with the geomechanical simulations of Zhang *et al.* (2007), or some combination of the two. In future work it will be important to calibrate the method against independent data from tracer tests, fluid pressure monitoring or 4D seismic observations, and compare the results with conventional reservoir simulation.

The term ‘critical’ is used in two senses in the literature. In the physics literature a critical point system is one with a correlation length comparable to the system size (Bruce & Wallace 1989). This is consistent with the long-range power-law correlation function, up to at least 10 km on the normalized plot of Figure 7. In the geological literature, ‘critically-stressed’ usually means a system that is *locally* near failure (Barton *et al.* 1995). ‘Self-organized criticality’ is a *global* concept (Bak & Tang 1989) where a large correlation length can be maintained while deformation is concentrated on larger fault systems. Simulations in anti-plane geometry show how fault populations can evolve to the critical-point system size by repeated earthquakes, resulting in ‘stress shadows’ in between the large faults, implying the system is not uniformly near failure everywhere in a local sense (Cowie e*t al.* 1993).

The statistical analysis presented above implies that the system, if not stationary, is generally evolving only slowly in time (see Fig. 7a). In other reservoirs there may be more dramatic systematic interventions that will introduce significant non-stationarity in the raw data (time series), implying that the method would be applicable only within time intervals that could be regarded as quasi-stationary.

The results presented here can be interpreted solely in terms of the subsurface response to flow rate change. This implies that surface operational changes, such as control of surface injection pressure, or the application of chokes on producers, do not ‘stack’ as coherently as the geological or geomechanical responses over time. Such effects may have been screened out by the significance tests, since they may be more likely to add a random component to the time series. Nevertheless, in future it will be important to use documented data on such interventions to see if any systematic effects may be observable, and hence to devise schemes to minimize their effect on the statistical reservoir model.

Here we have shown that the statistical reservoir model can in principle improve reservoir description by highlighting faults and fractures at different scale lengths, and improve short term recovery by directly forecasting production rates. Other applications may also be possible, including screening for geomechanical effects (by looking for long range correlations); subsurface monitoring of CO_{2} injection (looking for evidence of reactivated faults in time-lapse mode); or in continuous testing of incremental predictions of decline curve analysis.

## Conclusions

We have developed a new concept, the Statistical Reservoir Model, for calibrating hydraulically reactive faults and fractures, and predicting reservoir response to water flood solely through examining the flow rate statistics of injectors and producers. The results of Main *et al.* (2006) and those presented here confirm that long-range correlations between well pairs are associated with reservoir faults, as observed in some previous studies based on simpler techniques. Long range correlations are found over an optimal time lag of one month that cannot be explained by Darcy flow, but that are consistent with a geomechanical origin based on the poro-elastic mechanism operating in a near-critically stressed crust. The principal components from the correlation matrix reveal features that align in position and orientation with known faults where they intersect the main producing horizon. The existence of long-range correlations in hydrocarbon production is one more independent piece of evidence that the Earth's crust is in a near-critical state, and underlines the importance of geomechanics in modelling reservoir response to water flood.

The method has the potential to identify the most compliant structures that may be reactivated through geomechanical effects. This picture of the reservoir can therefore answer a very simple question: namely, when is it necessary to carry out a full geomechanical simulation to explain long-range effects and when is a normal drainage model (Darcy flow) sufficient? It can also be used alongside structural models developed from seismic and borehole data, notably to identify the most significant hydraulically reactive faults and fractures.

The model represents a new way of examining structurally complex reservoirs, giving independent information not contained in conventional deterministic reservoir models. More work is required to validate the model in different environments, but in principle the method could be used to validate deterministic models; to interpret active faults; to predict short-term production; to optimize offtake and injection rates; and to highlight where there is missing information, particularly in longer-range processes.

## Acknowledgments

This work began as part of NERC CONNECT grant GR3/C0022 with matching funding from BP, and was completed as part of the COFFERS project, under the Industry Technology Facilitator Complex Reservoirs Programme. We are grateful to sponsors: Hess, BG Group, BP, ConocoPhillips, DTI, Kerr-McGee, Statoil, Shell and Total, who supported this project by providing funds, data or robust and constructive criticism. We thank Statoil for providing the dataset used, but do not imply that Statoil is necessarily in agreement with any of our interpretations for the Gullfaks field.

- © The Geological Society of London 2007

## Appendix

The statistical method used in the Statistical Reservoir Model (Fig. 2) is explained here.

To determine the response of the reservoir to perturbations at well sites we minimize the mean prediction error
1
between the observed fluid flow rate *y*_{i, t} at the *i*'th producer for times *t*=2, *T* and that predicted, *ŷ _{i,t}*, by multivariate regression on a vector

*X*

_{t-k}of elements comprising the fluid flow rates

*x*

_{j,t-k}at all

*N*producers and

*M*injectors at time

*t*−

*k*, where

*k*is a lag time. The solution to (1) for all

*y*

_{i, t}is the Statistical Reservoir Model 2 where

*Ŷ*is a vector of predicted flow rates at all

_{t}*M*producers and

*R*_{k}is a matrix of the regression parameters. For more than one time lag

*R*_{k}would be a three-dimensional array with elements

*r*

_{i,j,k}:

*i*= 1,

*N*;

*j*= 1,

*N*+

*M*;

*k*=1,

*K*.

The inversion for the optimal Statistical Reservoir Model is done in two steps. First the well pairs that are significantly correlated at different lag times are identified using a Bayesian Information Criterion (BIC. Papasouliotis (2000), modified after Leonard & Hsu (1999), eqn 1.1.6). This removes well pairs that do not significantly contribute information. Pragmatically the search is stopped for a given producer when the multivariate regression coefficient is *R*^{2}=0.99 or a given number of iterations is completed. Second, Bayesian Dynamic Linear Modelling (Papasouliotis (2000), based on models presented in Leonard & Hsu (1999), sections 5.5 and 5.6) is used to eliminate a lower number of pairs whose optimal regression slope is not significantly different from zero. These two steps together define a binary significance matrix, *N*_{ij}, where most elements are zero, resulting in a parsimonious model. Typically only 5–25 out of the 106 wells in the field are needed to achieve *R*^{2}=0.99 for a given producer.