We evaluate the controls on production performance of the Wilrich tight gas sand play in West Central Alberta, and show that careful steering using 3D seismic to place the wellbore within the upper reservoir is the most important geophysical contribution to production outcomes. Geologic, geophysical, drilling, and production data from more than 20 wells are used in the analysis. The completion and production parameters within the study area are relatively invariant, creating a control experiment relative to other productivity factors. We thus isolate the effects of varying bottom hole pressures, porosity, wellbore length, number of stimulations, mud gas response, gamma ray measurements while drilling, mud weight, curvature, amplitude versus offset (AVO), amplitude versus azimuth (AVAz), velocity versus azimuth (VVAz), and position of the horizontal wellbore within the reservoir. These variables are treated separately and in a multivariate fashion to determine their relative and combined effect on the productivity of the wells. Several methods of statistical evaluation are used to test confidence in the results. The Wilrich sand is approximately 20-m thick, and it was expected that the multistage fracture stimulation would have minimized the importance of vertical permeability variations by adequately accessing the entire vertical reservoir section. Such is not the case; precise placement of the wellbore in the most permeable stratigraphy of the thin reservoir is of material importance. The pressure and porosity strongly affect the production performance, but to a lesser degree than vertical position within the reservoir. This suggests that stratigraphic concerns as they relate to permeability variation can be critical, even in thin fracture-stimulated reservoirs. Interesting relationships were observed between the AVAz and curvature measures, but neither they nor the AVO or VVAz attributes yielded statistically significant correlations to the production data.


Tight gas sand plays are an important subset of the resource plays available. Their permeability often plays an important role in productivity along with other factors such as mechanical properties and tectonic stress (Holditch, 2006). The Wilrich Member (Wilrich) tight gas sand play of West Central Alberta is a laterally extensive prograding marine shoreface complex of only moderate thickness and moderate-to-poor permeabilities (Zonneveld and Moslow, 2004). To achieve economic production rates with such permeabilities, horizontal drilling and multistage open-hole fracture stimulation is required. As with other resource plays, finding the Wilrich, and finding hydrocarbons is not in question (the Wilrich is gas-charged and locally pervasive). The goal, instead, is to maximize the productivity of each horizontal well to be drilled. Economic advantage in the resource industry fundamentally revolves around achieving the highest expected value from a given investment, which on a risked basis reduces to: producing for less cost, or to produce more for the same cost (Newendorp, 1975). The goal of the interpreting geophysicist in the Wilrich play is to determine the best placement of wells to maximize production. If this can be done in a sustained manner, it may represent a value-oriented strategic advantage (De Kluyver, 2000). It follows that two questions must be addressed: what controls the production of the Wilrich, and how can geophysical techniques be used to address those controls in an advantageous manner?

The geophysical relevance to oil and gas production is of natural interest and started with the earliest work on reservoir geophysics. Production starts at the reservoir, so reservoir geophysics must be the progenitor of production geophysics. Rummerfield (1954) first spoke about using the reflection quality to make indications about the reservoir. Hilterman (1975) summarized some of the early bright-spot era challenges and pitfalls of using amplitudes to interpret the reservoir. Their work has been followed by hosts of efforts to advance techniques that may be applied to understand the reservoir. Neff (Sheriff, 1992) illustrated the use of 3D seismic and synthetic modeling to map net pay. Amplitude versus offset (AVO) techniques (Ostrander, 1984; Shuey, 1985) were eventually used for hydrocarbon and lithology discrimination, and have been extended by many others in a variety of ways, including Goodway et al.’s (1997) AVO-derived lambda*rho and mu*rho (LMR) method. Avseth et al.’s (2005) work on quantitative interpretation illustrates how far the state of the art has come in regards to objectively relating rock properties, seismic, and the reservoir. The enthusiasm for quantitative interpretation and geomechanics expanded with the advent of resource plays such as shale gas and tight gas. Goodway et al. (2010) discusses the ability of LMR and AVO and amplitude versus azimuth (AVAz) to predict the geomechanical aspects of fracture stimulation as well as to reservoir quality. Sena et al. (2011) and many others have suggested similar methods of identifying so-called “sweet spots” where AVO analysis suggests the rock will be more geomechanically amenable to fracture stimulation. Close et al. (2012) quantitatively illustrates the relationship of LMR properties to the completion performance of wells in the Horn River basin. The strong correlation between these attributes and the instantaneous shut-in pressures of the multistage completion zones supports their use of these attributes to predict completion effectiveness. Close et al.’s (2012) work does not quantitatively illustrate the production performance of the wells being studied, which limits their argument to the realm of geomechanics.

The important step between either reservoir quality and production, or geomechanical properties (as they relate to completions) and production, is less well-studied. Da Silva and Marfurt (2012) suggest a multiattribute neural network approach to predicting estimated ultimate recovery (EUR) from LMR, AVAz, and other seismic data, but do not show actual results. Hunt et al. (2012b) do show the relationship of AVO, curvature, and AVAz attributes to completion data and production data in the Nordegg tight sand of West Central Alberta. This work shows correlations between seismic-derived porosity thickness (Phi-h) and production, as well as AVAz and production. The conclusiveness of Hunt et al.’s (2012b) work is limited by the small sample size of the validation data, which was the production log data from three horizontal wells. Alzate et al. (2012) also use production logging and 3D seismic to show an apparent relationship between LMR attributes and production. This experiment also had limited statistical data, but was similarly encouraging. Roth and Royer (2013) examine the controls on production in the Eagle Ford unconventional shale gas play of South Texas. This work uses an enormous publicly available database of 3500 wells to examine which variables affected production. They examine the effects of completion parameters as well as reservoir thickness, well length, proximity to faults, and reservoir depth, and employ nonlinear multivariate models in their analysis. Roth and Royer’s (2013) work shows that depth, or reservoir level build-up pressures, as well as reservoir thickness are key parameters that control production. This study does not lack for statistics, although the broad nature of the data limits its conclusions to being more regionally than locally applicable. There remains a paucity of studies of geophysical relationships to production rates for resource plays that have the necessary statistical size and detailed data to answer specific operational questions.

We show a specific case study for the Wilrich tight gas sand wherein sufficient production data exists to address the play specific controls on production with statistical confidence. We examine a variety of seismic attributes such as AVO, AVAz, and curvature, but also consider engineering and geologic attributes including well length, number of stimulations, bottom hole pressure, gamma ray log, mud log, and mapped phi-h. More than 20 horizontal wells were used in the study. All of these wells were drilled and completed in a similar manner. We collected production data in the same manner for each well. This means that fracture type, rates, size (75 tons for each frac), proppant properties, and flow conditions were all held essentially constant for the well sample. In addition, we employed the same well site geologic team and geophysical steering team throughout the project. The relative position of the horizontal wellbore within the Wilrich was regularly and finely estimated by this team, enabling a stratigraphic position code to be defined and treated as a ranked variable (Sokal and Rohlf, 1995) for all wells. The integrated geologic and geophysical approach was useful in unforeseen ways: geologic data from numerous vertical wells are used to estimate Phi-h when AVO methods proved to be compromised by the effects of variable overlying coals. Our results are unbiased toward technical discipline, and suggest that bottom hole pressure and phi-h are important to the productivity of the wells, but that the relative stratigraphic position of the well in the 20-m-thick reservoir is the most important single factor in production. It has been shown that much thicker sandstones may require multiple fracture stimulation to access the entire vertical section. Green et al. (1997) demonstrates that multistage stimulation yielded improved production of vertical wells in the Eagle Sand in Montana; however, the Eagle Sand is greater than 120 m thick, more than six times the thickness of the Wilrich. Our Wilrich study differs from the Eagle Sand in that it is pursued with horizontal wells, and at 20 m thick, fracture stimulations were expected to have height growth enough to encompass the entire vertical section. Previous sweet spot mapping for horizontal operations in the literature (Goodway et al., 2010; Sena et al., 2011; Alzate et al., 2012; Close et al., 2012) have been primarily concerned with optimizing the lateral position of wells using LMR methods, but we have examined the vertical position of the well within a relatively thin reservoir. The results show a surprising importance to the vertical position of the wellbore given the significant size of the multistage fracture stimulations used and the limited thickness of the reservoir. We examined the results in a variety of ways: as single variable regressions, as multivariate regressions, and as stepwise regressions (Davis, 1986). We have tested the confidence of the results with p-tests (Fisher, 1935; Kalkomey, 1997), with cross-validation tests (Hampson et al., 2001), and by examining the physical consistency of the multilinear weights. We have also used different length of production data to further examine our results. The conclusion that the vertical wellbore position was the dominant control on production in the Wilrich is confirmed by each of these tests. This work emphasizes the importance of careful steering of horizontal wells in a vertical (stratigraphic) as well as a lateral sense. It also suggests the importance of near well-bore permeability despite the use of multiple, significantly sized, fracture stimulations.

The Wilrich tight gas sand

The sequence stratigraphic framework, depositional origin, and lateral facies continuity of the Wilrich (or Falher G Member) was defined in outcrop and adjacent subsurface and was interpreted as a series of prograding shoreface complexes by Zonneveld and Moslow (2004). Figure 1 is a stratigraphic column showing the Wilrich as the oldest of the Falher Members within the Spirit River Formation. Due principally to its lack of conglomerate facies and associated high permeability, the Wilrich has been ignored until recently as a reservoir unit with potential economic significance for natural gas production. In Figure 1, the Wilrich is represented only by the deeper water shale element of the shoreface complex, in keeping with the legacy of pessimism over the zone. However, the advent of horizontal drilling, multistage completions technology, and a unique set of reservoir characteristics have combined to make the shallower water shoreface sandstone one of the most prolific tight gas resource plays in the Western Canadian Sedimentary Basin (Moslow and Zonneveld, 2012).

The Wilrich is characterized from core analysis as being of low-to-moderate porosity (6%–8%) and very low permeability (0.05 to 0.1 mD), which is consistent with other tight gas sandstone plays. Figure 2 shows a petrophysical type log for a typical Wilrich shoreface sequence. A coal section overlies the Wilrich reservoir, which has variable thickness, and very low velocity and density characteristics that have a strong impact on the reflectivity of the seismic. Core porosities and permeabilities are given by black dots. There is a good agreement between the log effective porosities and the core effective porosity, as seen in Figure 2b. Figure 2c illustrates the increasing dolomitic cement content lower in the Wilrich. Detailed petrographic analyses and core data shows a series of three reservoir units: a lower, middle, and an upper unit. The different units are separated by a dashed blue line in Figure 2.

The lower unit is comprised of very fine-grained sandstone with quartz and larger amounts of dolomite as cement and detrital grains. Intergranular porosity is common with low permeabilty. This rock acts as a storage system in this play. The middle unit is comprised of fine-grained sandstone with equal amounts of quartz and dolomite. The unit has lower porosity and permeability and acts as a baffle between the upper and lower unit. The upper unit is comprised of fine- to medium-grained sandstone with high quartz content. Intergranular porosity is common and permeabilites are higher. Core permeabilities in the upper unit, as shown in Figure 2d, are as high as 0.2 md, which is 2 to 10 times higher than in the middle or lower units. Figure 2e shows thin section representations of each unit. The larger grain size in the upper unit and lack of dolomitic cements drive its superior permeability. The upper unit is therefore the best reservoir and was chosen early in the play development as the horizontal target interval. As a consequence of the potential importance of the upper reservoir, we acted early in the project to evaluate our relative stratigraphic position within the Wilrich reservoir during drilling. All horizontal wells were monitored in real time with gamma ray and mud gas logs, with the seismic time-depth estimate, and by sample collection and evaluation. This integration of well-site and seismic information is used to interpret the position of each wellbore within the reservoir at regular intervals (5 m).

Correlation coefficients and linear regression

We use linear regression as our primary tool to predict production from the hypothesized engineering, geologic, and geophysical variables, and the study of the correlation coefficient as our primary means of evaluating the significance and potential usefulness of the predictions. Linear regression and the Pearson correlation coefficient are together perhaps the most important statistical techniques used in experimental science (Rodgers and Nicewander, 1988). These techniques are used pervasively in biometry and psychometrics, as well as in earth science. Linear regression and the correlation coefficient are tools to describe or define the relationship between two observed variables, X, and Y, which are a typically a sample statistic rather than a population parameter. The estimate of the sample means for X and Y are X¯ and Y¯, respectively, and there are n samples. The estimate of the standard deviations of the two variables are denoted as sx and sy, and the variances can be denoted as sx2 and sy2. The correlation coefficient is often written in this form 
which can also be written 
There are many ways of interpreting the correlation coefficient (Rodgers and Nicewander [1988] define 13 ways and claim there are probably more), but formula 2 has a straightforward interpretation. The correlation coefficient is the standardized covariance. The covariance of two variables is highest when they vary together: a high value of X yields a high value of Y, or the opposite. When the variables change in a connected fashion, they covary; the covariance is high. Division by the product of the standard deviations normalizes the covariance so that the size of the values does not affect the measure of interdependence. The values of the correlation coefficient are between ±1. Correlation coefficients do not require normal sample distributions (Rodgers and Nicewander, 1988).
Linear regression involving one or more variables is well-described in many statistical texts, and we briefly describe it so it can be related to the correlation coefficient. In linear regression, a dependent variable Y is being predicted by an independent variable X, and their functional relationship can be written (Sokal and Rohlf, 1995) 
where Y^ is the estimate for Y. The weights of the b terms bo and b1 are determined from the least squares normal equations, and have the objective of minimizing the error in estimating Y (Davis, 1986) 
Equation 3 can be expanded to include many variables, in the following way 
The multilinear weights can be solved using an expanded version of the normal equations used in single variable regression (Hampson et al., 2001).
A quality or goodness of fit measure of the line to the points Yi is defined by (Davis, 1986) 
where R2 is mathematically equivalent to the square of the correlation coefficient term r in equation 1, and is given the name coefficient of determination (Rodgers and Nicewander, 1988). The relationship between the minimization of error in linear regression and the correlation coefficient is thus demonstrated, and are useful tools in examining the association between variables.

Tests for significance: Observation and p-values

Several techniques should be used when deciding that a correlation between variables is significant. Some of the most important tests can be performed without mathematics. The first of these is that the data is related by a hypothesis rather than only a correlation coefficient. Our scientific method is, by its nature, an exercise in inductive confirmation, and as such a hypothesis must come before data. Hughes et al. (2010) emphasize that data alone do not an argument make. Observational studies notwithstanding, data comparisons should have some justification. The possibility of a correlation by chance or noise in the data becomes higher when we compare data with no reasonable claim of physical or causal relationship. The second technique that does not necessarily require mathematics is observation of the scatter plots of the data being related. Observation alone may suggest that the data has a nonlinear relationship, which should be addressed, or that the correlation is being driven by a very few outliers (Davis, 1986). Davis (1986) demonstrates statistical methods of doing this, called tests of model, but such inquiry can often be done by simple observation.

Fisher (1935) introduced significance testing of correlation coefficients. His approach was to test the null hypothesis that the correlation coefficient could have been found by chance. Fisher’s null hypothesis, applied to the correlation coefficient, is stated 
Statistical work is performed to evaluate the strength of the evidence against the null hypothesis in the form of a probability. The probability being measured is called a p-value, which represents the probability of observing a correlation coefficient as far from zero as the one observed, under the assumption that the null hypothesis is true. If the p-value is achieved, statisticians (and we) suggest that the null hypothesis should be rejected (Davis, 1986). Fisher suggests a value of p-value of less than 5% as his significance level, but this is completely arbitrary. Stern and Smith (2001) suggest that much smaller p-values should be used; however, they are referring to enormous observational medical studies with hundreds of variables. Choosing too high a p-value introduces the possibility of a chance correlation being interpreted as meaningful, while too low of a p-value can lead to the false dismissal of too many potentially important relationships. In this study, we consider the 5% p-value, because the multivariate nature of the production question led us to expect smaller single variable correlation coefficients, and it was felt to better to be more inclusive than exclusive in the study. There are other more sophisticated hypothesis tests, such as the Neyman-Pearson method (Stern and Smith, 2001) which involve alternative hypothesis testing. We did not feel this more sophisticated approach was likely to be helpful or practical at this stage of inquiry. The relationship between correlation coefficients and p-values depends also on the sample size. The larger the sample size, the more certain the correlation coefficients become, and the lower the correlation coefficient must be to pass a given p-test. We define the p-values with a two-tailed sample size dependent test statistic called the Student’s t-distribution, which looks a lot like a normal distribution, except that it has a wider tail (Davis, 1986). The smaller the sample size, the broader the t-distribution becomes. The so-called t-test is very commonly used in evaluating correlation coefficients, and is what Kalkomey (1997) used in her discussion on correlation coefficients for geophysical inquiry. The t-test for significance of the bivariate correlation coefficient is 
The t-test p-values can be looked up in tables, such as in Davis (1986), although many statistical programs will also calculate them. Multivariate t-tests for significance can also be found in Kalkomey (1997). The 5% p-value for the 21 samples we used in most of our work occurs at: r=0.433, and a correlation coefficient of 0.549 is required to pass the 1% p-test.

Stepwise regression

We expect that production is a multivariate problem involving pressure, permeability, and other variables. Moreover, we expect that some of the key production-controlling variables themselves can only be predicted indirectly through the use of several variables. We must therefore eventually employ multivariate analysis. We employ simple multilinear regression as defined by equation 5. The key question is how do we decide which variables to use in our multilinear regression model? We employed a method called forward stepwise regression (Davis, 1986). The method is described as follows:

  1. Test all variables for correlation coefficients with single variable regression.

  2. Choose the single variable with the highest correlation coefficient. This variable is called the first variable, and its regression formula defines the best single variable production prediction model.

  3. Perform a p-test. If the p-test is passed, continue.

  4. Perform two variable regression, using the single variable chosen in step 2 and, in succession, each of the remaining variables. The variable, that when used with the first variable, gives the highest correlation coefficient is chosen as the second variable. These two variables and their regression formula define the best two variable production prediction model.

  5. Perform a p-test. If the p-test is passed, continue to the best three-variable test.

This process continues until the p-test is failed.

Cross-validation testing

Another method of estimating the correct number of attributes to use in a multivariate study is called cross validation. The method differs from p-testing in that it evaluates the actual predictive capabilities of the regression formulas. The idea is to evaluate the effect of overtraining and outliers on the predictive solution by dividing the sample into two groups, one in which the multilinear weights are calculated, and another, hidden group upon which the weights are tested in a predictive fashion. The predictive error is calculated for the hidden group. The prediction error we calculate is the total squared error in prediction from equation 4 divided by the number of samples, making it an root-mean-square (rms) average of the error 
We express the results in terms of a percentage error from the sample mean of production. Similarly to Hampson et al. (2001), our method is a leave-one-out cross validation, which we employ in conjunction with stepwise regression. The method follows:

  1. Take the best single variable.

  2. Choose one well from the sample as the hidden well.

  3. Calculate the regression weights using all wells except for the hidden well.

  4. Use the regression formula to predict production for the hidden well.

  5. Record the predicted value and the prediction error for the hidden well.

  6. The hidden well goes back into the regular sample group and another well is chosen as the hidden well.

  7. Repeat the calculation of regression weights with the new sample and apply them to predict the production of the new hidden well.

  8. Record the predicted value and the prediction error for the new hidden well.

  9. Repeat this process until all wells have been the hidden well one time.

  10. Calculate the rms average prediction error for the single variable cross-validation test.

  11. Perform the same process for the best two variables from the stepwise regression analysis.

  12. Calculate the rms average prediction error for the two variable cross-validation test.

  13. Continue this testing process for the best three-variable solution, and so on.

  14. Plot the rms average error versus the number of variables in the solution. When the error bottoms out, or its downward slope is significantly reduced, the optimum number of variables has potentially been found.

The predictive error calculated in formula 9 always decreases with increasing number of variables when the same sample is used for regression analysis and for prediction. With cross validation, this is not the case. Prediction error calculated with our method of cross validation will typically bottom out or reach a minimums when the number of variables in the solution is optimum (Hampson et al., 2001). Additionally, the regression weights are calculated n times, and can be examined for consistency.

Darcy’s law as hypothetical justification

We stated earlier that the variables that we test should be supported by a hypothesis. Darcy’s law provides broad justification for the variables we evaluate. Darcy’s law suggests that production will vary as a product of pressure differential, permeability, the inverse of viscosity, and an area term (Holditch, 2006). We ignore the viscosity term because gas composition is relatively invariant locally, and there is no water production. We include pressure in our analysis, as pressure does change in the area. The horizontal drill and multistage fracture-stimulation method is predicated on increasing the exposed area of the reservoir, which supports our recognition of horizontal length and the number of fracture stimulations performed. Permeability is very difficult to discuss with tight gas reservoirs (Holditch, 2006), and may have additive components such as matrix and fracture permeability terms (Nelson, 2001). Of additional difficulty relative to permeability is that nothing has been done to directly measure it. Many of the variables we test are evaluated because of their potential to infer some element of permeability. Phi-h (or porosity) may or may not infer permeability correctly, as Figure 2 demonstrated, whereas position alone in the reservoir may also be inadequate. Any AVO attributes that predicted porosity could be justified in this test. Gamma ray or mud log values may have a porosity-permeability indicative relationship and may therefore be tested. The fracture inferring seismic parameters of AVAz, velocity versus azimuth (VVAz), and curvature may infer fracture presence, although there is no calibration data in the area to validate this or to transform any of those measures to permeability. For the sake of completeness, we chose to test all of the variables mentioned, including the AVO measures, on the basis of Darcy’s law. One particular measure, λ/(λ+2μ), may have a relationship to fracability and therefore area, as well as a potential relationship to porosity (Close et al., 2012). Although Darcy’s law is multiplicative, the many variables to be considered might be treated in a multilinear fashion in recognition of the additive components of permeability. The multilinear approach respects the empirical nature of this work, and is also easy to use even if it is not perfect. Encouraging results from this method could be used to inform a more sophisticated linearization of Darcy’s law at a later date.

3D seismic initial observations and vertical wells

Three-dimensional surface seismic covers most of the six-township Wilrich play area, with numerous well ties. A total of 151 vertical wells penetrate the Wilrich and tie the 3D data; these wells were subsequently used, along with the 3D seismic, for time-to-depth conversion. The top of the Wilrich directly underlies the coal sequence illustrated in Figure 2, and is close to the trough-peak zero crossing on seismic. The Wilrich is typically picked as the peak directly underlying the coal sequence, and we call that peak the Wilrich peak or the sand peak. The Wilrich reservoir velocity is nominally 4700m/s, and with the 37.5 Hz dominant frequency in the area, is approximately one-sixth of a wavelength in thickness. This thickness is midway between Widess’ (1973) thin bed thickness of one-eighth of a wavelength and the tuning thickness of one-quarter of a wavelength. Figure 3 shows a synthetic seismic image from the same well used for the petrographic work performed earlier. The wavelet used is an SEG standard zero phase 10/1560/70Hz wavelet. Interpretive windowing of the Wilrich zone has been executed in recognition of the thickness of the zone and the synthetic image of Figure 3, and typically includes the coal trough above the Wilrich, as well as the Wilrich peak. The trough just below the Wilrich, which we call the sand trough, is also considered in the windowing schemes in an effort to avoid the effects of tuning from the overlying coal. Figure 3 clearly shows that rock property variations in the overlying coal and silt sections have an enormous magnitude, and may dominate the seismic response. Windowing schemes must vary for different seismic properties in recognition of their different physical meaning and calculation parameters. Figure 4 illustrates the different windowing methods we used. The reflectivity and migrated stack data use horizon referenced windows of sample width, as illustrated in Figure 4a by the arrows. Attributes with this sample method have the name “AMP” associated with them. A pick is made for the coal trough, the sand peak, and the sand trough. Figure 4b and 4c show the VVAz Vani and K1 curvature attributes, which are defined by a single pick close to the Wilrich peak pick. The VVAz method is horizon-interval-based, so any sample within the horizon interval is sufficient, whereas the curvature calculation method uses long-time windows and is somewhat insensitive to the pick site. Figure 4d illustrates that the AVO LMR attributes use a large window that spans the coal trough to the sand trough. Arithmetic averages, maximum, and minimum values are calculated within these windows. A smaller windowing for AVO attributes that excludes the coal trough is not included in this paper because it failed to yield better results. Figure 4e illustrates that the AVAz Bani attribute uses a larger rms average window, defined according to Hunt et al. (2012a). This windowing scheme is used for the horizontal well extractions as well as these vertical well extractions.

Most of the 3D data was processed with a 5D minimum weighted norm interpolation algorithm prior to prestack time migration imaging and AVO analysis as in Hunt et al. (2010). AVO gather preparation and measurement used time variant gather trimming, adaptive super gathering (Xu and Chopra, 2007a), and NMO-stretching corrections (Xu and Chopra, 2007b). The AVO attributes are handled isotropically. LMR properties were subsequently estimated, as were AVAz anisotropic gradient (Bani) (Ruger, 1996), VVAz velocity anisotropy (Vani), and curvature attributes including the K1 principal curvature using the methods of Chopra and Marfurt (2007, 2010). The attributes were correlated to geologic control data from the vertical wells with the same method as Hunt et al. (2010). A total of 114 vertical wells were drilled within the AVO-processed data area. Table 1 shows the correlation coefficients for the well to seismic comparison.

The correlation coefficients given in Table 1 show primarily that all of the attributes, with the exception of curvature, enjoy a statistically significant correlation with the coal isopach. Every attribute fails the 1% p-test for statistical significance relative to the Wilrich sand. Only the AVAz Bani attribute manages to pass the 5% p-test for statistical significance relative to the Wilrich reservoir. The AVO LMR attributes did not show a materially improved independence from the effects of the overlying coals despite the fact that they are estimated from inversion (Goodway et al., 1997). These results suggest that Phi-h mapping of the Wilrich based upon the seismic data alone should be treated with skepticism.

The poor performance of the AVO attributes in this test could have several causes, with the first being the dominance of rock property variation of the overlying coal section. Other contributing factors could include the data quality itself, problems in the processing, or a lack of variation in the Wilrich sand rock properties, which would be a practical companion of the large variations in the coal sequence. Figure 5 illustrates the effect of varying coals and varying Wilrich sand. The original analysis well is used, and is given the name “Well A.” Another well with much poorer Wilrich sand and the same logging program was chosen, and given the name “Well B.” Well B has a thicker sand, but much of it is tighter and because those tight elements do not meet the 3% porosity cut-off, Well B has a poor Phi-h value compared to Well A. The variation in the coal sequence is typical for wells in the area. A variation of Well A (called “replaced Well A”) was also produced wherein only the Wilrich sand from Well B replaces the sand section in Well A. The Spirit River Formation isopach was kept constant for this comparison, which is supported by the Well B Spirit River isopach. This replacement well is placed in the middle of the model of Figure 5a. These three wells therefore represent a control experiment relative to changes in the Wilrich sand as well as changes in the coal sequence. Figure 5b illustrates the normal incidence synthetic model of the cross section. The Well A seismic character has barely detectable differences with the replaced Well A character, but the replaced Well A character is quite different from the Well B character. It is clear that the coal is the dominant effect on the seismic amplitudes. The problem appears to be obvious: when the coals change, compressional and shear velocities as well as density change significantly, whereas when the sand changes, these rock properties have very small relative changes. With the coal directly overlaying such a thin sand interval, it would appear that this problem is physically challenged with data of this frequency content.

We also investigated the AVO response of the change in sand by creating a Shuey (1985) model for the replaced Well A log with exactly the same parameters as we did for Well A earlier. Figure 6 shows a comparison of the AVO gathers of the original Well A and the replaced Well A. The differences are barely detectable to the eye. Figure 7 shows a crossplot of the AVO amplitudes from the center of the Wilrich peak in the models. Because only the sand is changed, we compare to the AVO with the identical scale bars. The difference in intercept is a barely perceptible 2.4% relative to Well A, and the gradient has less difference. We therefore conclude that the problem with the AVO attribute-Phi-h prediction is likely a fundamental problem with the physics and geology.

The poor correlations of AVAz, VVAz, and curvature to the Wilrich reservoir Phi-h were not discouraging to further examination of these attributes. Their use for anisotropic identification, particularly for fractures (Lynn et al., 1996; Ruger, 1996; Chopra and Marfurt, 2007; Hunt et al., 2012b), remained unevaluated from the Phi-h test. Moslow and Zonneveld (2012) suggested the Wilrich would be commonly fractured in this area, so numerous maps were created with these attributes at the Wilrich level. Figure 8 illustrates a key map from the evaluation of the anisotropic attributes: the study area corendering of K1 curvature and AVAz Bani at the Wilrich level. The use of opacity allows the higher values of each attribute to be seen simultaneously. Note the dominant trends are in the northwest to southeast direction, which is the direction of maximum horizontal stress in the area. These two attributes have a low but statistically significant correlation with each other: 0.285 on 51,477 samples. The correlation is qualitatively apparent in Figure 8, and appears to be related to extensional draping over underlying Devonian carbonate (reef) features, which are present and have been drilled in the area. The correlation of these attributes appears encouraging; however, their meaning relative to production is not assured because they are correlated with each other. The observed relationship between the anisotropic attributes suggests that they might be of interest to production, motivating their use later in the production evaluation.

Data-handling method for horizontal wells

Twenty-four horizontal wells were drilled. The same geologic, geophysical, and engineering team worked on all the wells. Three-dimensional seismic-derived time-depth estimates were used in conjunction with wellsite data to steer the wells. All wells were completed with the same type of fracture (frac) stimulation, including very similar frac rates, fluids, proppant types, and frac size or tonnage. The nominal frac size was 75 tons. The wells flowed into a proprietary gathering system. Line pressures were known and were held essentially constant over the study area. Initial test rates, which are short preproduction flow tests, and actual initial production (IP) rates were collected with the same protocols. Table 2 summarizes that production data, excluding liquids. The Wilrich within the study area produces no water, and consistently low ratios of condensate in the 10 standard barrels per million cubic feet range.

Longer production data is preferable for analysis. All of the data was used, although the 60-day average initial production (IP 60) was used the most heavily as it had the best trade-off of statistics (there are 21 wells with IP 60 data) versus length of test. The same production data was also divided by number of fracture stimulations to give an equivalent (IP) per frac set of measures. Some might argue that the analysis method should only use IP per frac data, instead of unscaled IP data. Both measures are used to avoid prejudicing the results by the assumption that production scales according to the number of stimulations. The wells are open hole completed, so production may not be entirely stimulation driven. Additionally, by using unscaled IP data, the dominance of all variables can be compared, including the number of stimulations. A statistical description of the IP 60 production data is given by the histogram in Figure 9. The data has an average of 4.3mmcf/d, the median is 4.4mmcf/d, and a range of 7.9mmcf/d. The central tendencies of the data are clear, although it is unclear if the IP 60 distribution is Gaussian; the expectation from Holditch (2006) is that tight gas reservoir permeability and production distributions should be log-normal. Of immediate concern to the operator is the fact that the lower production values (less than 2mmcf/d) experienced by some of the wells would not be economic.

Wellsite geologic data such as gamma ray, rate of penetration (ROP), and mud gas were collected at submeter samples. Actual rock samples were collected every 5 m, and the relative stratigraphic position was estimated every 5 m along the horizontal wellbore. Stratigraphic position, or position, is defined by a code with one of three values at each 5 m sample point: one if the wellbore is in the upper reservoir unit, two if the wellbore is in the middle reservoir unit, and three if the wellbore is in the lower reservoir unit. Geophysical data from the 3D were gridded and mapped every 20 m to their correct position along each wellbore. The gridded and mapped data were collected and converted to log format as in Hunt et al. (2012a) and displayed with the geologic wellsite data. Phi-h data came from geologic gridding of the 151 vertical wells in the study area. Geologic log-based Phi-h was used as opposed to geophysical estimates of Phi-h due to the poor seismic to Phi-h relationship. The vertical well density was sufficient to map the porosity trends reasonably. This Phi-h data was mapped to the horizontal wellbores, collected, and displayed as log data along with all of the other information.

A variety of other geologic and engineering data is similarly gridded, and mapped to the horizontal wells. This data includes bottom hole pressures for the Wilrich, which were taken from the vertical well control data. Pressure is expected to be a key variable, so its inclusion is critical. Given that the gathering system line pressures are essentially invariant, bottom hole pressures are used as the pressure term in the production evaluation. The bottom hole data from each vertical well was gridded throughout the area. The pressure grid values were mapped every 20 m along the horizontal wells to produce a log, in a similar fashion to the other gridded variables. This method should represent the pressure accurately because the Wilrich has such limited thickness. Coal thickness and Phi-h from overlying Falher channels were also assigned with the gridding method. We evaluated the potential effects of these overlying zones under the hypothesis that the fracture stimulations could potentially access them. A variety of curvature attributes were used, including a very long wavelength parameterization (Chopra and Marfurt, 2007), that was dubbed “S-long.” Euler curvature (Chopra and Marfurt, 2012) was also used. AVO LMR attributes were collected with the same average, minimum, and maximum, windowing scheme from the Phi-h evaluation described earlier.

Final step: Wellbore averages

The final step in the method is to bring the position, engineering, geologic, and geophysical data together with the production data. Because all of the production data is in the form of single IP numbers for each well, a single wellbore average for every well must be made for each of the other data measures, too. This means that the submeter sampled gamma ray, mud gas, and ROP data, the 5 m sampled position data, and the 20 m interval mapped values of pressure, Phi-h, AVAz Bani, VVAz, curvature, and AVO attributes, are each averaged to give mean values for each horizontal well. If production logs existed for every well, the log data display format could be maintained, but there is no production log data for these wells. This situation is common at present. This wellbore averaging approach is not only reasonable, but likely to be inevitable for anyone attempting to use a large statistical sample for production analysis in a local area. With this done, all of the variables are compared to the IP data in a single and a multivariate fashion.


The first comparison is a simple correlation matrix, and is illustrated in Table 3. This matrix gives the correlation coefficients for each of the potential production controlling variables with the actual IP data. Because the IP 60 data is to be used heavily in further analysis, the correlation coefficients for it with the other IP data are also shown. We use the 5% p-test as our measure of statistical significance. IP 60 has a statistically significant correlation with all of the other production measures. The other test variables are rank-ordered according to the magnitude of their correlation coefficient with IP 60 per frac, but are separated into two groups: the geologic, engineering, and wellsite variables in the first group, and the seismic variables in the second group. Very few of the variables have statistically significant correlation coefficients. The only variable with consistently significant correlation coefficients to the various IP data is the position variable. The position correlation to the IP data always improved, if only slightly, when the per frac scaled version of each IP period was used.

The petrographic analysis, permeability, and grain size observations led to the hypothesis that position might be important, but it was not expected to be the most dominant control on production. Further investigation of the data relationships was performed by creating scatter plots of the position variable with the IP data. Figure 10 illustrates of these plots. Figure 10a shows the position versus IP 60 scatter plot. The negative correlation is consistent with the prediction of greater permeability higher in the reservoir, which corresponds to a lower position value. Figure 10b shows the same comparison except that IP 60 data is now scaled by the number of fracture stimulations. The frac-scaled scatter plot is clearly more natural-looking and better correlated. Figure 10c illustrates the position variable plotted with the 24 wells with test data scaled by frac and also has a believable qualitative appearance. Figure 10d compares position against the one-year frac-scaled average production. This plot has the most impressive and obvious correlation, and passes the 1% p-test for statistical significance. These scatter plots support the correlation coefficient data in suggesting that the position of the wellbore is related to production of the wells. The fact that this position relationship holds for long and short periods of production has significant economic ramifications to steering within the upper reservoir.

Considerations of Darcy’s law suggest that several other variables should be considered for further analysis, despite having failed the 5% single variable p-test. Figure 11 shows scatter plots from four of those variables. Figure 11a shows the number of frac stages versus IP 60. The positive correlation is in accordance with the expectation that a greater number of stimulations will yield better production rates. Vertical striping in the data suggests that there might be a lack of range and distribution of the number of fracs, but also that the number of fracs is not likely to be the dominant variable relative to production. Figure 11b shows the well length versus IP 60. This scatter plot does not support any argument for longer well bores, and its negative trend is against expectations. Figure 11c is geologic Phi-h versus IP 60. The positive slope of the data trend is in accord with expectations of a positive porosity permeability trend. Figure 11d is the bottom hole pressure versus IP 60. The positive sloping trend is congruent with Darcy’s law. This scatter plot is the best distributed of the four shown in this illustration. The correlation coefficient for pressure almost (but not quite) passes the 5% p-test. Failing the 5% p-test does not prove that these variables are not related to production; however, it does warn of the possibility that the apparent correlations are a random rather than an informative event.

Given that we know production has multivariate controls, it is reasonable to expect that single variable correlations may look poor. We therefore undertake the stepwise multilinear approach of evaluation. We allow all variables to be considered and ranked, including the geophysical variables. Figure 12 shows the correllation coefficient analysis of the stepwise multilinear regression for IP 60. P-tests are also performed for each variable. Variables added after the three-variable solution failed the p-test. A vertical (purple) dashed line is drawn in Figure 12 to illustrate the p-test failure point. The three-variables that this analysis supports as production controls are: position, Phi-h, and pressure. Stepwise regression with the IP 30 and IP 7 production data yields the exact same choice of the first five variables in the same order. The one-year average IP data does not have enough data points to perform stepwise regression. The seismic LMR, AVAz, VVAz, and curvature variables do not pass the 5% p-test in this evaluation.

Multivariate analysis, such as the exhaustive stepwise regression approach used here, evaluates enough combinations of the data variables that it can overfit the solution (Kalkomey, 1997). We therefore employ the leave-one-out cross-validation method we described earlier. Figure 13 illustrates the cross-validation prediction error analysis. The full data prediction error solution is depicted in red, and it decreases as more variables are added. The cross-validation prediction error is shown in purple. It is larger than the full variable cross-validation error and drops until it comes to a local minimum, as Hampson et al. (2001) also describe. In the case of the Wilrich production data, the local minimums occurs at four variables, which is when the 5% p-test failed. This is the point at which skepticism should be applied to the notion of using additional variables. The variables of position, Phi-h, and pressure remain the only three that this testing supports.

There are other fundamental methods of evaluating the multivariable solutions. The simplest remaining method is by inspection of the solutions. The scatter plots for position, Phi-h, and pressure are examined in Figure 10 and Figure 11. None of these scatterplots reveal a dominant outlier. As a single variable solution, only position passes the 5% p-test. Multilinear estimates of IP 60 are shown in Figure 14. Figure 14a shows the multilinear estimate of IP 60 as estimated using position and pressure plotted against the actual IP 60 data. The correlation coefficient for this solution is 0.578. Figure 14b shows the multilinear estimate of IP 60 as estimated using position and Phi-h plotted against the actual IP 60 data. The correlation coefficient for this solution is 0.640. Figure 14c shows the multilinear estimate of IP 60 as estimated using position, Phi-h, and pressure, plotted against the actual IP 60 data. The correlation coefficient for this solution is 0.753. Each of these scatter plots have a natural appearance. They are spread well across the plotted space, and there are few remarkable outliers. The three variable solution is clearly the best in appearance. The sign of the multilinear weights for each solution is another element of information that can be examined. From our understanding of Darcy’s law and the Wilrich reservoir, it should be expected that the position weight should be negative, whereas the weights for each of Phi-h and pressure should be positive. This is checked for each of the multivariable solutions: All of the weights have the expected signs.


The analysis of the wells and their production data suggests that the position of the wellbore within the Wilrich reservoir is the most dominant control on production. It is not surprising that the position in the reservoir would affect production, but it is surprising that this variable is dominant within the control group of wells. There are several potential contributing factors to this result, the first of which is that the position variable was well sampled. Many wells were drilled in the upper, middle, and lower parts of the reservoir. The converse to this point is that some of the other key variables had limited range. The number of fracture stimulations, for instance, only varied between 7 and 13. A limited range reduces the perceived value of the production control. Similarly, the relationship of AVAz, VVAz, and curvature to production might be limited by a lack of actual fractures variability in the study area. Although it was expected that there would be fractures in the Wilrich, it is possible that there were not many, or that most of them had limited permeability. The correlations between these variables is interesting, particularly AVAz anisotropic gradient and curvature, but does not appear strongly relevant to production in this data. Another reason for the poor performance of some of the variables could be that they are too inaccurate. The correlation coefficients of stack amplitudes and LMR attributes to Wilrich Phi-h is strongly indicative that these variables would fare no better as indicators of production. If the seismic cannot accurately predict Phi-h or mechanical properties, there is no reason to expect it to predict production. This problem is specific to the local Wilrich with its thick, variable, overlying coal section. The geologic estimate of Phi-h that was used may have had reduced predictive use due to its lack of accuracy between the vertical well control. Improvements in seismic acquisition and processing that would overcome these accuracy problems have the potential to make a material difference in the prediction of production.

Value of precise steering

An economic assessment of the potential value of precise steering within the Wilrich can be made using the analysis performed in this work. The economic assessment is in the form of a simple model for production, using the three-variable multilinear estimate for IP 60 production. To isolate the value of steering, we hold the Phi-h and pressure variables to be constant at the 21-well average. Estimates of IP 60 production can then be made using position values of one, two, and three, and are shown in Table 4. These position values span the possible range. The differences in expected IP 60 are material. Although it is quite likely that wells drilled in the upper reservoir will be economic at 5.52mmcf/d, it is less likely to be so for wells drilled in the middle reservoir. Wells drilled in the lower reservoir are not expected to be economic.


The analysis of the Wilrich production data shows that the variables of position, Phi-h, and pressure have the clearest affect on production. Position is the only variable that consistently passes the many tests for statistical significance with the various lengths of production data. The strong correlation of position to production for short and long production periods has clear economic ramifications for the value of precise steering in the upper Wilrich reservoir. Multivariate analysis confirms that the variables of position, Phi-h, and pressure, have plausible weighted interactions in creating a superior prediction of production. Other variables such as the number of fracture treatments or well length were expected to have a positive and significant correlation with production, but their affect may have been reduced because of the limited range of their variation in the well group. AVO and the other seismic methods likely failed because of the effects of coals and the weak AVO variability of the Wilrich.

The numerous statistical tests undertaken in this study yield strong conclusions that are specific to the study area and methods of drilling and completion. More general conclusions can only be made by considering the expected reasons for the observed behavior. Geoscientists should consider what in situ variations exist, and if they have a relationship to flow through Darcy’s law, or to fracture-stimulation behavior. If permeability varies significantly, as in this case, properties and measurements that infer permeability and may direct drilling toward it, may be useful for maximizing production. In this case, seismic was used to aid in positioning the well. In other cases where stresses or fracturing are more relevent, other approaches may have more value.


The authors would like to thank Pulse Seismic Inc. for permission to display images from 3D seismic data licensed from them. We would also like to thank Santonia Energy for its support of this work.

Biographies and photographs of authors are not available.

Freely available online through the SEG open-access option.