Seismic correlated Mallik 3D gas hydrate distribution: Effect of geomechanics in non-homogeneous hydrate dissociation by depressurization

Journal of Natural Gas Science and Engineering 20 (2014) 250e270

Contents lists available at ScienceDirect

Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse

Seismic correlated Mallik 3D gas hydrate distribution: Effect of geomechanics in non-homogeneous hydrate dissociation by depressurization Mafiz Uddin a, *, Fred Wright b, Scott Dallimore b, Dennis Coombe c a b c

Alberta Innovates e Technology Futures, Edmonton, Alberta T6N 1E4, Canada Geological Survey of Canada, Natural Resources Canada, Sidney, British Columbia V8L 4B2, Canada Computer Modeling Group Ltd., Calgary, Alberta T2L 2A6, Canada

a r t i c l e i n f o

a b s t r a c t

Article history: Received 18 May 2014 Received in revised form 3 July 2014 Accepted 3 July 2014 Available online

The delineation of the Mallik gas hydrate field has utilized extensive well logging and substantial 3D seismic testing and interpretation. This study explores the use of seismic data to quantify the areally heterogeneous gas hydrate distribution. The available Mallik 3D seismic data was compiled and compared/contrasted with available well log data from two adjacent wells. Based on the seismic information, two areally variable (i.e. non-homogeneous) scenarios for gas hydrate distributions are considered: Scenario I having the same initial total hydrate amount as our earlier model areally uniform (homogeneous) distribution, and Scenario II with significantly less overall total hydrate, but honouring the same relative distribution. The scenarios of variable gas hydrate distributions are used in dynamic simulations of the lower Mallik zone. Simulations of each were conducted with and without the role of geomechanics. In Scenario I, we observed multiple gas production peaks (which quite similar to 6 days production behaviour) with higher localized pressure pulses occurred due to strong gas hydrate heterogeneity. In Scenario II, this drastic change in gas production rate was not observed (due to faster pressure evolution in the reservoir). In both Scenarios, the overall reservoir gas production peak is delayed compared to the homogeneous case. This is further delayed by the role of geomechanics. More interestingly, all simulation cases show a very similar overall production trend. This is probably a unique for the Mallik gas hydrate production using single vertical well, including a gas production peak but terminating in a stabilized period of lower but significant gas production. With geomechanics, gas production in general and the gas production peak is shifted and delayed. The geomechanics effect is not purely compaction drive (as in conventional reservoirs, gas production increases with geomechanics). The simulations utilized two set of geomechanical parameters obtained from logs (dynamic parameters) and rocks testing (static parameters). Geomechanical responses based on dynamic parameters were essentially equivalent to simulations ignoring geomechanical effects. The geomechanics simulations indicate an essentially elastic reservoir response (i.e. no plastic failure) assuming a cased vertical well. The Mallik upper zone A and middle zone B are closer to the permafrost and nearer to plasticity limits should be explored. Crown Copyright © 2014 Published by Elsevier B.V. All rights reserved.

Keywords: Gas hydrate Depressurization Hydrate dissociation Seismic Geomechanics Simulation

1. Introduction 1.1. The Mallik gas hydrate deposits The Mallik gas hydrate field is located in the Mackenzie Delta on the coast of the Beaufort Sea, Northwest Territories, Canada. * Corresponding author. Tel.: þ1 780 450 5047; fax: þ1 780 450 5242. E-mail address: mafi[email protected] (M. Uddin). http://dx.doi.org/10.1016/j.jngse.2014.07.002 1875-5100/Crown Copyright © 2014 Published by Elsevier B.V. All rights reserved.

In 1998, the JAPEX/JNOC/GSC Mallik 2L-38 gas hydrate research well was drilled as part of a CanadaeJapan collaboration to characterize a concentrated gas hydrate deposit and to evaluate its energy potential. A second multidisciplinary research and development programme was undertaken in the winter of 2001-02, and incorporated three new wells (JAPEX/JNOC/GSC 3L-38, 4L-38 and 5L-38) consisting of a production well (5L) and two monitoring wells (3L; 4L). Work included extensive well logging and 3D seismic testing and interpretation. Industry-acquired 3D seismic reflection

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

data were used to characterize the occurrence and spatial extend on the Mallik gas hydrate accumulation (Bellefleur et al., 2010). The Mallik 3D seismic coverage area is shown in the surface geology map, Fig. 1, including subsurface faults; undifferentiated moraine, lacustrine, marine, and estuarine deposits; fluviatile silt, sand, and gravel; Iperk Sequence and Kugmallit Sequence. More than 60 scientific and technical research papers documenting this work were published as Geological Survey of Canada Bulletin 585 (Dallimore and Collett, 2005). The Mallik gas hydrate resource is characteristic of a Class II hydrate deposit with hydrate co-existing with water. The gas hydrate intervals have high gas hydrate saturation values that, in some cases, exceed 80% of the pore volume (Dallimore et al., 1999; Dallimore and Collett, 2005). CanadaeJapan collaboration at Mallik culminated in March 2008 with the achievement of 6 days of continuous gas production via depressurization from a 12 m perforated interval, at an average gas rate of about 2000 m3/day, and a total produced gas volume of approximately 13,000 m3 (Kurihara et al., 2012; Uddin et al., 2012). Three drawdown stages were conducted, illustrating interesting combined kinetics and flow behaviour. In particular, after each pressure drawdown, rapid gas pulsing followed by flow stabilization was observed (Wright et al., 2011; Uddin et al., 2011). A summary of much of the recent Mallik research activities including simulations of the production test can be found in the Geological Survey of Canada's Bulletin 601 (Dallimore et al., 2012). 1.2. Gas hydrate production simulations As part of the Geological Survey of Canada Arctic Gas Hydrate Production Research Project, an extensive numerical modelling

251

component was pursued in several stages, based on the information gathered from the Mallik gas hydrate site. To match the observed 2008 production test, our earlier studies developed a multi-layer radial model (areally homogeneous, vertically heterogeneous) and a simple kinetics approach to model gas evolution and transport. Long term (8 years) production performance was also estimated with the model, assuming certain boundary conditions (Uddin et al., 2011, 2012). More recently, a simulation evaluation of the productivity of all three hydrate zones of the Mallik test well was conducted (Uddin et al., submitted for publication). 1.3. Objective and approach This study first incorporates both available seismic and well log information into a Cartesian static model for the Mallik hydrate deposit for lower zone C. A revised estimate of hydrate saturation is provided for this zone C. The Cartesian static model is then converted to an equivalent 3D radial model centred around well 2L-38. Consideration of the “tilted” and areally non-homogeneous nature of the hydrate deposit, observed from the seismic data, is incorporated into the model. Two limiting model scenarios of variable gas hydrate distributions are used in dynamic simulations of the lower Mallik zone. Simulations of long time production for each scenario were conducted without and with the role of geomechanics. More particularly, Section 2 first reviews and correlates well logs information for the wells 2L-38 and 5L-38 in the area of interest. Reasonable correlations between the two wells (separated by 50 m) are found for the lower two zones B and C. Section 3 explores available seismic data to quantify the areally heterogeneous gas hydrate distribution. Consistent with this

Fig. 1. Regional surface geology map showing major subsurface tectonic elements. Area of the Mallik 3D seismic coverage shown (red outline) within the Taglu fault zone (TFZ). Qu (dark yellow): undifferentiated moraine, lacustrine, marine, and estuarine deposits; Qf (light yellow): fluviatile silt, sand, and gravel; T1 and T6 (brown): Iperk Sequence and Kugmallit Sequence, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

252

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

seismic information, three separate static models of Mallik zones A, B, and C were generated. Gas hydrate volumes were estimated for each of the three zones. Analysis showed that this data is missing expected aquifer zones however, probably because correlations with available well logs were not considered. A simple method is then applied to incorporate both available seismic and well log information into a Cartesian static model. Thereafter, only the lower gas hydrate zone C is considered further. A revised estimate of hydrate saturation is provided for this zone. Next, in Section 4, the Cartesian static model from Section 3 is converted to an equivalent 3D radial model centred around well 2L38. Consideration of the “tilted” and areally non-homogeneous nature of the hydrate deposit, observed from seismic, is incorporated into the model. Different scaling of the seismic data results in two model scenarios of variable gas hydrate distributions that are used in dynamic simulations of the lower Mallik zone. Simulations of long time production for each scenario are conducted without and with the role of geomechanics. Sections 5e7 present the simulations results for two model scenarios of variable gas hydrate distributions. Here, the effects of geomechanics in gas hydrate dissociations are simulated for both scenarios and compared with simulations without geomechanics. The choice of appropriate geomechanics parameters is explored in more detail in Appendix A. Here a simplified, essentially homogeneous, radial model is used to explore variety of geomechanical factors via a mechanistic sensitivity analysis. 2. Mallik hydrate bearing formation Geological data shows that the Mackenzie Bay Sequence extends from 350 m to 926 m below surface. It includes Upper Oligocene to Miocene sediments consisting of unconsolidated to weakly cemented sandstones and siltstones with rare shale interbeds. No hydrocarbon occurrences have been documented within the Mackenzie Bay Sequence. The Kugmallit formation was found at 926 m, and extends down to 1310 m. It consists of Oligocene sediments of unconsolidated to weakly cemented sandstones with interbeds of siltstone and shale. A number of regional hydrocarbon fields occur in Kugmallit strata, especially in offshore areas. At the original Mallik L-38 site, both solid gas hydrate and possible free gas were interpreted between 890 and 1106 m. Permafrost at Mallik L-38 extends to ~640 m. Ice bonding is quite variable through this sequence and talik or non-ice bonded sections are expected. Seismic and well logs data identified significant hydrate deposits in all three zones. The upper zone extending from 890 to 940 m occurs entirely within the Mackenzie Bay Sequence. Gas hydrate saturations range from 50 to 85%. The middle zone extending from 945 to 1035 m is a complex interbedded section comprising a series of 5e10 m thick gas hydrate bearing sand units separated by 0.5e1 m thick gas-hydrate-free silt layers. Gas hydrate saturations range from 40 to 80%. The lower zone (Fig. 2) extending from 1060 to 1120 m consists of two thick hydrate occurrences with many similarities to the upper gas hydrate bearing sand. Gas hydrate saturations are very high, ranging from 80% to 90%. The continuous well logs (2L-38 and 5L-38) for different variables plotted in a compacted form can be found in the GSC Bulletins 585, 601 (Dallimore and Collett, 2005; Dallimore et al., 2012). 2.1. Mallik lower hydrate bearing formation Continuous geophysical well logs were interpreted for the lower C zone hydrate bearing formation. Figs. 2 and 3 summarize some of the log variables in this C zone from wells 2L-38 and 5L-38. This data was extracted from more extensive logging information

provided by the Geological Survey of Canada. Comparing the lower C zone profiles between these two wells (Fig. 2), reasonable correlations between the two wells can be found, including the vertical profiles of high hydrate saturations and the separating coal layers. One should also note the clear presence of a bottom aquifer zone. More particularly, Fig. 3 illustrates the sonic well log data which is used subsequently for our geomechanical simulations. The 2007 and 2008 tests focused on production from the bottom zone only. This zone is characterized by an overlying clay layer, two regions of extensive hydrate deposits, and underlying free water zone at the base of the hydrate stability region. Indeed, the lower contact of gas hydrate zone, at a depth of 1107 m, marks an abrupt change from gas hydrate bearing sand to gas hydrate free sand. The pressure and temperature at the base of the hydrate stability zone are 12,962 kPa and 12.2  C, respectively. Wright et al. (2005) predicted that the pressure and temperature conditions for the base of the gas hydrate stability zone at 1107 m are consistent with in situ pore-water salinities of 45 ppt (parts per thousand). 3. Seismic based GH simulation cube This section presents seismic correlated 3D gas hydrate distribution in a selected simulation cube. Bellefleur et al. (2010) applied an acoustic impedance inversion technique to 3D true-amplitude seismic data acquired over the Mallik area to characterize gas hydrate occurrences. They constructed multi-dimensional heterogeneous models of petrophysical parameters (P- and S- wave velocities, and density) of hydrate-bearing sediments. Bellefleur et al. (2010) generated 3D gas hydrate concentration distribution in the vicinity of the Mallik gas hydrate wells by combining the modified BioteGassmann theory (Lee, 2002), statistical parameters estimated from well logs, and the horizontal correlation length ranged between 100 m and 500 m. The modified BioteGassmann theory of Lee (2002) relates density and velocities to hydrate concentration as VS ¼ Vpa(1  (1  c)∅). Here, VS, Vp are the formation (m) velocities; a is the matrix velocity ratio, V(m) S /Vp ; c is the gas hydrate concentration, and ∅ is the porosity. Fig. 4 shows the results of these calculations for an areal view of the inverted P-wave velocity for zone C. This view illustrates the heterogeneous nature of the hydrate distribution observed from the seismic processing, and as such, provides quantitative information useful in developing heterogeneous simulation models of the Mallik hydrate deposits. Our previous works (Uddin et al., 2011, 2012) focused on honouring near-well log information on hydrate distribution, but no information on the areal extent on the hydrate deposit was considered. In this work, we turn to an analysis of available seismic information in order to establish its utility in quantifying areal hydrate extent. Starting with seismic cube provided by Bellefleur et al. (2010), three independent simulation models for upper, middle and lower hydrate bearing formations are setup. The simulation area of 1 km2 (1000 m  1000 m) was chosen equivalent to our earlier 500 m radius radial grid model. The production well 2L-38 is located at the centre. The simulation grid sizes are 10 m, aerially and 1 m, vertically. 3.1. Setup simulation cube e without honouring well logs Fig. 4 shows the location of hydrate concentration cube (processed hydrate saturation data) received. This consisted of hydrate saturation values, the UTM coordinates areally, and a time value vertically. Data resolution was 21e22 m areally, and 0.002 s vertically. From this, a numerical simulation cube was setup, utilizing a simulation domain of (1000 m  1000 m), between (X ¼ 770,531e7,706,386; Y ¼ 512,840e513,911). Here the production well 2L-38 is located at the centre of the simulation domain.

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

253

Fig. 2. Comparing the well logs (2L-38 and 5L-38) for the Mallik lower gas hydrate bearing zone C e the well log 2L-38 showing total 11 log variables and the well log 5L-38 showing 2 key variables (formation materials and hydrate saturation). The wells 2L-38 and 5L-38 are approximately 50 m apart.

Fig. 5 illustrates the transformed simulation area. For this area, a simulation grid spacing of 10 m areally and a vertical grid spacing of 1.0 m was utilized. Three independent simulation models vertically were generated, consisting of upper (Z ¼ 890e940 m), middle (Z ¼ 940e1058 m), and lower (1058e1133 m) hydrate zones. This study presents the lower C zone gas hydrate production simulation. The hydrate saturation values for the simulation grid were generated from the original data utilizing simple interpolation, based on the areal coordinates of the original versus simulation grid locations. As illustrated in Fig. 4, a portion of the desired simulation grid (NE corner) had no equivalent original seismic data. These (missing seismic values) data was generated using a moving average approach. Thereafter, following the CMG STARS modelling approach, hydrate saturation was converted to hydrate (solid)

concentration values by multiplication with solid (hydrate) mole density. Next, void porosity distribution for the entire simulation field is generated, based on the log values at the centre of the grid, adjusted by the hydrate saturation areally for each layer (when Sh > 40% assume sand void porosity, when 20% > Sh > 40% assume sand/silt void porosity, when Sh < 20% assume silt void porosity). With hydrate saturation and void porosity defined, the fluid porosity values are generated by the simple formula ∅F ¼ ∅0(1  Sh). Here, Sh is the hydrate saturation and ∅0 is the void porosity. Computer programmes were written to perform all these seismic data manipulation steps. Fig. 5 illustrates the generated gas hydrate concentration and effective fluid porosity for this zone. Based on the earlier quoted

254

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

Fig. 3. The Mallik lower gas hydrate bearing zone C e the well log 5L-38 showing the variations in hydrate saturation, geomechanical properties and seismic reading.

formula, these are seen to be correlated properties. Also the closest corner of these cubes is seen to be relatively less heterogeneous than the remaining regions of the cubes. This is a result of the moving average procedure used to generate the “missing” seismic data, as explained above.

3.2. Setup simulation cube e with honouring well logs This section applies a straightforward method to incorporate both seismic and well log information in a Cartesian model of the lower Mallik C zone gas hydrate reservoir. The same Cartesian grid

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

255

Fig. 4. 3D seismic cube by Bellefleur et al. (2010) and the numerical simulation cube, (a) maps of the inverted P-wave velocity for zone C (horizon map), zone C reveals a complex pattern with large variation in P-wave velocity over distances less than a few hundred metres, (b) numerical simulation cube of 1.0 km2 area equivalent to the earlier 500 m radius base case model (production well 2L-38 located at the centre of the simulation area).

utilized in the previous section is used to define the 3D model cube (1000 m  1000 m, with 10 m sized grid blocks). Vertically 75 m is simulated, via 1.0 m sized grid cells. As the model is centred on well 2L-38, the corresponding well log information (such as porosity and hydrate saturation) is employed to define the vertical heterogeneity of the model. The Mallik C zone vertically averaged seismic definition of hydrate saturation is used to generate hydrate areal heterogeneity. The seismic data is then used to rescale property values areally as follows. Grid gas hydrate saturation at any given grid blocks (i, j, k) is defined by the following expression: ðLÞ

Sh ði; j; kÞ ¼ Sh ðkÞ

ðSÞ Sh ði; jÞ ðSÞ Sh ð50; 50Þ

Here, S(L) h (k) is the hydrate saturation from well (2L-38) log at layer k; S(S) h (i, j) is the seismic gas hydrate saturation for C-zone; i, and j are the block indices along x, and y directions, and S(S) h (50, 50) is the

hydrate saturation value at the well log (2L-38) block. An additional check to restrict the rescaled hydrate saturation below 95% was also enforced. The grid block effective fluid porosity ∅F is defined as ∅F ¼ ∅0(1  Sh). Grid block effective permeability kE is then obtained using CarmaneKozeny type function as:

kE ð∅F Þ ¼ k0 ð∅F =∅0 Þck ½ð1  ∅0 Þ=ð1  ∅F Þ2 Following the above procedure, the resulting spatial distribution of hydrate concentration and fluid porosity are shown in 3D in Fig. 6. Here the top-most and bottom-most zones representing shale layers show a uniform areal distribution (based on zero hydrate saturations from the well logs), while the intermediate hydrate containing layers show the seismically corrected areal distributions. Analogous areal distributions in two selected layers containing hydrate are illustrated in Fig. 7. Here the profiles are seen to be correlated, which is a result of the scaling procedure with the average seismic distribution of zone C.

Fig. 5. 3D seismic simulation cube for the Mallik lower gas hydrate zone C (simulation area of 1 km2, gas hydrate production well 2L-38 is located at the centre), (a) gas hydrate distribution and (b) porosity distribution.

256

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

Fig. 6. Seismic and well logs correlated 3D static model showing the hydrate concentration distributions (data in resolution of 10 m  10 m  1 m grid spacing obtained based on moving average of 2 well logs (2L-38 and 5L-38) data with normalized seismic signature).

With this static model as a basis (honouring both seismic and well log data), the next section will turn to dynamic simulation gas hydrate production predictions. In an effort to minimize long run times and to compare with our earlier studies, these simulations will utilize a radial grid model. The conversion of the Cartesian model to an equivalent radial model will also be discussed. 4. Building 3D radial simulation model This section builds a 3D radial grid model for dynamic simulation of the gas hydrate production. This has several advantages. First connection to our previous based case 2D radial model (Uddin et al., 2012) becomes more apparent. As well, our initial 2D

geomechanics analysis presented in Appendix A can also be directly utilized and compared. The choice of radial gridding also has a significant numerical advantage. Run times are expected to be several orders of magnitude faster than an equivalently discretized Cartesian model. The main goal, however, is to study how an initial non-homogeneous spatial distribution of hydrate saturation will influence the predicted results. More particularly, we consider here two such possible spatial distributions, termed Scenario I and Scenario II. The former takes the original spatial seismic distributions generated and rescales them to achieve a total initial hydrate in-place consistent with the earlier 2D base case model. By so doing, we hope to generate a direct assessment of spatial inhomogeneity. This rescaling however

Fig. 7. Seismic and well logs correlated 3D static model showing hydrate and porosity distributions, areally, in layers 36 and 38, (a) hydrate concentration and (b) porosity distributions (data in resolution of 10 m  10 m  1 m grid spacing obtained based on moving average of 2 well logs (2L-38 and 5L-38) data with normalized seismic signature).

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

257

Fig. 8. 3D radial grid model, (a) plain view showing grid blocks in i- and j-directions, and (b) vertical cross-section showing grid blocks in i- and k-directions.

results in maximum hydrate saturations of 100% in some regions, which we feel is physically unreasonable. The latter case attempts to address this concern by rescaling all hydrate saturations downwards (but preserving the same relative distributions) such that maximum saturations are approximately 80%. This second scenario therefore has a significantly smaller total hydrate in-place value. Finally, before describing these simulations in detail, we should close by recognizing the limitations of the approach we are taking. The limited angular discretization results in discontinuous spatial distributions angularly. This therefore can have some physical and numerical consequences. The method is still robust enough to capture the major implications of a spatially non-uniform hydrate distribution, which is what we are attempting to establish. 4.1. Radial grid model setup A radial grid model (with outer radius 500 m) was designed to cover the same areal extent as the Cartesian model, and centred around well 2L-38. The uniform areal (x, y) discretization of 10 m in

the Cartesian model was replaced by a varying grid size radially and 8 angular subdivisions. The conversion of properties areally between the two grids (Cartesian to radial) requires simple interpolation, due to the varying grid sizes. The vertical discretization was maintained in both models at 1.0 m thickness. However, in order to compare with our earlier models, the lower 13 m of the original Cartesian model was ignored. The net result is grid dimensions of 101  101  75 in the original Cartesian grid and 287  8  62 in the radial model. Fig. 8 illustrates the generated radial grid model. Finally the varying elevation of Mallik zone C gas hydrate deposit, also found from the seismic data (Bellefleur et al., 2010) was approximately incorporated into the radial model, utilizing the choice of angular subdivisions in the model. Relative to the centre point of the radial grid (well 2L-38), the elevations at the centre points of each angular subdivision at 500 m were established at the base of the hydrate stability zone. This was done by analyzing one vertical seismic section across 1000 m, plus the average areal seismic definition of hydrate saturation for the 1000 m  1000 m

Fig. 9. 3D dynamics radial grid model showing grid blocks hydrate concentrations in layers 16 and 36 (3D static model in Fig. 6 is converted into 3D dynamics radial grid model).

258

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

Fig. 10. 3D dynamics radial grid model showing the grid blocks hydrate concentrations in 4 vertical cross sections, 1e5, 2e6, 3e7 and 4e8 (3D static model in Fig. 6 is converted into 3D dynamics radial grid model).

area. Then via interpolation, the elevations at every 5 m radially at each angular subdivision were established, thus forming a definition of the complete bottom surface of the base of the hydrate zone. These calculations were performed in Excel resulting in the creation of an “index file for bottom hydrate elevations”. This surface was then extended vertically to the top of the hydrate bearing zone, based on the values at the well (innermost grid cell) to complete the description of the “tilted” hydrate deposit thickness. Then, based on the defined radial/angular/vertical grid cell dimensions, an array of hydrate saturation values for each grid cell was created as an input file for the simulations. A separate computer programme was written to read the Excel index file and to generate the array of “tilted” hydrate saturations for the simulations. For completeness, it is emphasized that this method could be readily extended to any future Cartesian model incorporating “tilted” hydrate deposits. In this case, the Excel index file defining the bottom hydrate surface elevation would consist of elevations every 5 m in both x, y directions (i.e. 100  100 values), and the computer programme would read this file and generate the array of “tilted” hydrate saturation values for every Cartesian grid cell in the model. Fig. 9 illustrates the results of these calculations for two layers (16 and 36) of the radial model. Because of the (limited) angular discretization and the overall tilting of the hydrate deposit, the generated hydrate saturations have a “blocky” appearance. As such, this represents a trade-off between physical information and Table 1 Initial material in place for two hydrate deposition scenarios in 3D radial grid model. GH reservoir

Material in-place Formation pore volume (m3) Hydrate volume (m3) Water volume (m3, STD) PeT conditions Average pressure (kPa) Average temperature ( C)

2D model (areally uniform gas hydrate distribution)

3D model Scenario I

3D model Scenario II

1.493  107

1.569  107

1.581  107

5.399  106 9.33  106

5.216  106 1.047  107

4.511  106 1.130  107

11,005 11.38

11,005 11.38

11,005 11.38

efficient numerical simulations of the dynamic depressurization process. This representation does however capture the heterogeneous nature of the hydrate deposit and incorporates consistently the available information provided by both seismic and well logs. Fig. 10 shows various cross sectional views of the generated “tilted” hydrate deposit (radial cross sections 1e5, 2e6, 3e7 and 4e8). Here the views are seen to be much more continuous and confirm the basic correctness of the methods used to generate them.

Table 2 Basic model parameters for the Mallik lower zone C (static and dynamics simulation parameters). Model basic properties

Values

Field log data e well 2L-38 Lower GH bearing zone (1060e1120 m) Perforation interval (1093e1105 m) Pressure at the mid perforation (kPa) Temperature at the mid perforation ( C) Absolute permeability for GH layers (mD) Absolute permeability for bottom aquifer (mD) Void porosity for GH layers (%) Void porosity for bottom aquifer (%) GH saturation (%) Water saturation (%) Salinity at the base of GH zone (ppt)

60 12 11,005 11.378 0.10e850 10e1800 15e45 20e40 10e85 100 (Class II GH) 45

Numerical parameter e radial grid model Modelling radius (m) Vertical layers (e) Number of areal grid blocks (e) Areally layer permeability (permeability variance, s2) Permeability anisotropy (kv/kh) Grid fluid porosity, fF (%) a Grid effective permeability, kE (mD) Layer GH saturation, horizontally Critical gas saturation for GH layers Capillary and relative permeability curves

500 57 287 Mild variation (s2 ¼ 0.10) 0.20 fF ¼ f0(1  SH) CeK function Varied (seismic based) 0.05e0.12 Van Genuchten (1980) and Parker et al. (1987) (calibrated with production data history match)

a Effective fluid permeability any given time and grid block was updated via CarmaneKozeny function (CeK function).

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

The resultant material-in-place values for the lower C gas hydrate zone are summarized in Table 1. Two scenarios combining seismic and log information have been considered. Also reported for comparison are the equivalent values for the radial base case 2D model (i.e. “angularly homogeneous”) model utilized in our earlier modelling work. More particularly, Scenario II utilizes directly the areal seismic data to rescale the areally uniform vertical layer data derived from the well logs in the angularly homogeneous model. As this procedure does not guarantee correct log values at the inner radial grid cell, Scenario I further rescales the seismic data to consistently match the well log data at the inner radial grid cell. This is analogous to the method described in Section 3.2. It is noted that this rescaling results in a total hydrate volume for Scenario I that essentially matches the total volume in place of the base case areally-uniform model. Thus a comparison of the Base case and Scenario I will illustrate the roles of a non-uniform areal hydrate distribution, while comparison with Scenario II illustrates the additional effect of a lower overall hydrate volume. Scenario II can be considered a lower limit estimate of the Mallik lower zone hydrate deposit. The basic model parameters for this Mallik lower gas hydrate zone C are given in Table 2.

259

of hydrate (2D and 3D models, respectively). The effects of geomechanics in 3D on the gas and water production behaviour are also considered, although the geomechanics spatial distributions will not be discussed here to shorten the presentation. Rather, more specific and detailed effects of geomechanics will be presented for Scenario II (see Section 6). 5.1. Gas and water productions Fig. 11(a) contrasts the 2D and 3D gas production rates and cumulative gas predicted by the two models up to 4.5 years of production. The 3D gas rates are somewhat higher than the 2D rates up to the maximum time simulated. In principle, either earlier or later rates might be expected, and observed behaviour is a consequence of the specific distribution profile utilized. Also compared are the simulations with and without geomechanics in 3D (the 2D simulation ignores geomechanics). Initially geomechanics is shown to increase gas production. However it is possible that at later times, the case without geomechanics may actually overtake the run with geomechanics and produce higher gas rates. This comment is based on our understanding of the

5. GH production results e Scenario I This section compares the predicted production behaviour of the areally-homogeneous and areally-heterogeneous distributions

Fig. 11. The effect of spatial gas hydrate distributions and the geomechanics on (a) gas and (b) water productions over the length of 8 years depressurization (numerical simulations are, 2D homogeneous (Uddin et al., 2012), 3D non-homogeneous (Scenario I in Table 1) with and without geomechanics).

Fig. 12. The effect of spatial gas hydrate distributions and the geomechanics in overall reservoir volumes in place, i) hydrate-pore volumes and ii) gasewater volumes over the length of 8 years depressurization (numerical simulations are, 2D homogeneous (Uddin et al., 2012), 3D non-homogeneous (Scenario I in Table 1) with and without geomechanics).

260

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

competing roles of geomechanics for freely flowing phases (hydrodynamics), and its effects on the kinetics of gas release. (The 3D simulations without geomechanics have not been run to full completion because of excessive run times.) Section 6 discusses the equivalent simulations for Scenario II where the non-geomechanics 3D run does overtake the run with geomechanics at much later times. The numerical simulation also shows frequent gas production peaks in 3D non-homogeneous gas hydrate distribution case (highlighted in Fig. 11(a)). Here two regions of rapid gas production for case (c) and one region for case (b) are observed. We emphasize these are a physical results as we repeated the calculations with different time step choices and obtained reproducible behaviour. The origin is an inhomogeneous trapping of large amounts of free gas which builds up until it is ultimately released. This is reminiscent of what was observed more locally several times during the Mallik six days field test. As well, a recent field test in Alaska exhibited a similar effect (unpublished). Occasional numerical oscillations also occur in these simulations but are instantaneous (one time step) behaviours having no effect on cumulative production levels. Fig. 11(b) analogously contrasts the 2D and 3D water production rates and cumulative water predicted by all three models. The 3D water rates are again somewhat higher than the 2D case up to the maximum time simulated. It is also seen that the water rates are less sensitive to geomechanics. In principle, for all three cases either higher or lower rates might be expected and observed behaviour is once again a consequence of the specific distribution profile utilized. Observed (but lower) water production spikes are tied to gas production behaviour.

5.2. Variations in material in place The changes in hydrate and pore volumes are shown in Fig. 12(i). It should be noted that hydrate volumes and pore volumes are related quantities since hydrate decomposition increases available fluid pore volumes. Overall there is a similar trend in 2D and 3D behaviour without geomechanics (the scales have been adjusted to emphasize differences). However, with geomechanics the time behaviour has a distinctly different pattern. Longer run times are needed before further conclusions may be drawn. Water and gas volume evolution are shown in Fig. 12(ii). Again, at early times, the slopes with time for 2D and 3D cases without geomechanics are very similar, while the gas evolution with geomechanics is quite different. However at later times, all curves are tending to be more similar. 5.3. Hydrate and gas saturation distributions Fig. 13(a) compares the hydrate distributions for the 2D model, the 3D model without geomechanics (Sections 4e8), and the 3D model with geomechanics (Sections 4e8) at simulation time of 780 days. This illustrates the evolution from a localized dissociation to a more extended one. The 2D initial distribution is by definition homogeneous areally, while the 3D initial distribution has a significant high saturation zone at the left side of the shown sections in some vertical layers. This is significant as this zone is much more resistant to dissociation, causing the 3D evolutions to be noticeably asymmetric relative to the 2D sections. The 3D case with geomechanics shows a slightly smaller dissociation region at all times, relative to the 3D case without geomechanics.

Fig. 13. (a) Hydrate concentrations and (b) gas saturations at the end of 780 days (2D radial model of areally homogeneous hydrate, 3D radial model with non-homogeneous hydrate). (All vertical length scales 60 m; all radial length scales 500 m from central well.)

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

261

Fig. 14. (a) Pressure, and (b) temperature distributions at the end of 780 days (2D radial model of areally homogeneous hydrate, 3D radial model with non-homogeneous hydrate). (All vertical length scales 60 m; all radial length scales 500 m from central well.)

Fig. 13(b) compares the gas saturation distributions for the 2D model, the 3D model without geomechanics (Sections 4e8), and the 3D model with geomechanics (Sections 4e8) at time of 780 days. This again illustrates the evolution from a localized dissociation to a more extended one. Initial gas distributions for all cases appear very similar and symmetric. At later times, the profiles extend laterally more and more, but such that the 3D case with geomechanics is somewhat less extended than the other two. This probably implies that the peak production peak in gas rate observed in Fig. 11(a) for the 2D case around 2300 days will occur much later in the 3D case with geomechanics, and possibly somewhat earlier than the 2D case for the 3D case without geomechanics. We are currently extending these simulations to confirm this speculation. 5.4. Pressure and temperature distributions Fig. 14(a) compares the pressure distributions for the 2D model, the 3D model without geomechanics (Sections 4e8), and the 3D model with geomechanics (Sections 4e8) at time of 780 days. This illustrates the pressure evolution to final state with a defined very low pressure region around the well surrounded by a vertically stratified distribution. The time evolution of these 3 cases leading to a similar final state is somewhat different, however. In particular, the 3D case with geomechanics (i.e. a “softer reservoir”) shows only the evolution of a localized low pressure region near the well. The “harder reservoir” cases, in 2D and 3D, behave similarly with an evolving high pressure region above the open zone of the producing well at earlier times, possibly until 1½ years approximately. This may be connected to the tendency to trap free gas in the reservoir, as illustrated in Fig. 13(b).

Fig. 14(a) compares the temperature distribution profiles for the 2D model, the 3D model without geomechanics (Sections 4e8), and the 3D model with geomechanics (Sections 4e8) at time of 780 days. This illustrates the temperature evolution to final state with a defined low temperature region around the well caused by the endothermic nature of the hydrate dissociation process. The time evolution of these 2 cases without geomechanics is remarkably similar, and includes a lower temperature region away from the well at later times in the region containing free gas. This effect is not seen in the 3D case with geomechanics, possibly due to the less extended free gas region. Also the role of the bottom aquifer in maintaining a higher temperature near the bottom of the well is observable from the temperature distributions in all 3 cases. Finally, we studied the average field pressure and temperature evolution for each simulation cases. An overall, there is decline in both pressure and temperature. The early increase in pressure for cases without geomechanics is a result of a buildup of free but trapped gas until it is able to flow. This is a local effect relatively near the well. This is only observable as an “average” pressure because the remaining regions of pressure are unchanged from their initial values. At later times, additional equivalent local pressure effects are expected to occur, but are unobservable as an average effect. Geomechanics is seen to dampen such pressure responses. The average temperature decline is a result of the endothermic nature of hydrate dissolution. 6. GH production results e Scenario II This section compares the predicted production behaviour of the 3D areally-heterogeneous distributions of hydrate with and

262

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

without geomechanics for Scenario II. Here we will also extend the time of simulation to probe expected long-time areal boundary effects. As shown in Table 1, this case has significantly less initial gas hydrate than Scenario I. From the simulations presented in Appendix A, the no-geomechanics case, or rigid geo-model, is approximately equivalent to a geomechanics run with mechanical parameters chosen from sonic well logs (i.e. dynamic parameters). The geomechanics parameters chosen in this 3D simulation correspond to the static parameters used in Appendix A (i.e. from laboratory testing). 6.1. Gas and water productions Fig. 15(a) contrasts gas production rates and cumulative gas predicted by the two 3D models, with and without geomechanics, this time up to 10.5 years of production. The gas rates curves cross at times 700 days and again at 2100 days. This variation is probably due to the overall competing field effects of kinetics and hydrodynamics behaviour, each affected differently by geomechanics. The “softer” geomechanical reservoir somewhat inhibits pressure changes and thus delays the kinetic gas dissociation rate It is probable that the gas rate peak occurs earlier for 3D without geomechanics. This was observed in 2D geomechanics simulation in Appendix A. We have confirmed that the gas spike at 250 days is numerical artefact, and could be removed by decreasing the time step size during this period. Note that the gas production rate peaks

Fig. 15. The effect of geomechanics on (a) gas and (b) water productions over the length of 8 years depressurization (numerical simulations are: 3D non-homogeneous (Scenario II in Table 1) with and without geomechanics).

at 68,000 m3/day around 3500 days and then plateaus at approximately 64,000 m3/day for the case without geomechanics. A somewhat delayed peak gas rate is indicated for the case with geomechanics. As discussed in our earlier work, this peak in gas production is an outer boundary effect (i.e. a limited areal extension of hydrate). Fig. 15(b) contrasts water production rates and cumulative water predicted by the two models. Here, both cases show very similar production trend with a crossover in rates after 3000 days. The deviation gradually increases with time. Again in all cases, we observed a very similar water production trend for both cases. The geomechanics water rates are noticeably lower. 6.2. Variations in material in place Fig. 16(i) contrasts the changes total in place hydrate volumes and pore volumes. Note that hydrate volumes and pore volumes are related quantities since hydrate decomposition increases available fluid pore volumes. Overall there is a similar trend in the 3D behaviour with and without geomechanics (the scales have been adjusted to emphasize differences). The variation in total hydrates volume (and available pore volume) for the first 500 days are different, but after 500 days, there is no significant variation.

Fig. 16. The effect of geomechanics in overall reservoir volumes in place i) hydratepore volumes and ii) gasewater volumes over the length of 8 years depressurization (numerical simulations are: 3D non-homogeneous (Scenario II in Table 1) with and without geomechanics).

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

263

Fig. 17. Gas hydrate distributions at the end of 780 and 2880 days, 3D model (seismic based areally non-uniform GH distribution) simulation, (a) & (c) without geomechanics and (b) & (d), with geomechanics. (All vertical length scales 60 m; all radial length scales 500 m from central well.)

Water and gas volume evolution are shown in Fig. 16(ii). Again, at early times, gas and water evolution with geomechanics is quite different. However at later times, all curves are tending to similar slopes. 6.3. Hydrate and gas saturations distributions Figs. 17 and 18 compare the hydrate and gas saturations distributions at two simulation times for the 3D geomechanics model (Sections 4e8). The 3D case with geomechanics shows a slightly smaller dissociation region at both times relative to the

3D case without geomechanics. When compared to Scenario I results at 780 days (Fig. 13), the Scenario II case shows a somewhat less extended profile in the layers containing the open well production zone but otherwise are very similar (note the colour scales, in the web version, for Scenario I and II are somewhat different). 6.4. Pressure and temperature Figs. 19 and 20 compare the pressure and temperature distributions at two times for the 3D geomechanics model (Sections 4e8). The

Fig. 18. Gas saturation distributions at the end of 780 and 2880 days, 3D model (seismic based areally non-uniform GH distribution) simulation (a) & (c) without geomechanics and (b) & (d) with geomechanics. (All vertical length scales 60 m; all radial length scales 500 m from central well.)

264

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

Fig. 19. Pressure distributions at the end of 780 and 2880 days, 3D model (seismic based areally non-uniform GH distribution) simulation (a) & (c) without geomechanics and (b) & (d) with geomechanics. (All vertical length scales 60 m; all radial length scales 500 m from central well.)

3D case without geomechanics shows the pressure and temperature pulses expanded in a considerably larger region relative to the 3D case without geomechanics. When compared to Scenario I results, Fig. 14, at 780 days, the Scenario II case shows a somewhat greater extended profile in the layers containing the open well production zone but otherwise are very similar (note the colour scales, in the web version, for Scenario I and II are somewhat different). The localized pressure pulse spreading is considerably higher in this Scenario II due to lower initial hydrate deposit and higher permeability. The average pressure and temperature evolution for these cases were studied. Overall, there is decline in both pressure and temperature. The early increase in pressure for cases without geomechanics is a result of a buildup of free but trapped gas until it is able to flow. This is a local effect relatively near the well. This is only observable as an “average” pressure because the remaining regions of pressure are unchanged from their initial values. At later times, additional equivalent local pressure effects are expected to occur, but are unobservable as an average effect. Geomechanics is seen to

dampen such pressure responses. The average temperature decline is a result of the endothermic nature of hydrate dissolution. 7. Geomechanical results e GH distribution Scenario II 7.1. Stress Cross-sectional views of the effective radial and vertical stresses changes over time are illustrated next. Fig. 21 shows the stress distributions predicted for the static moduli case, with the largest magnitudes in effective stress occurring near the open well intervals. The stress distributions predicted for largest times simulated, showing a wide regional change in effective mean stress by these times. 7.2. Volumetric strain and vertical subsidence Analogously, cross-sectional plots of volumetric strain distributions are considered next. Fig. 22 shows an earlier (t ¼ 780 days)

Fig. 20. Temperature distributions at the end of 780 and 2880 days, 3D model (seismic based areally non-uniform GH distribution) simulation (a) & (c) without geomechanics and (b) & (d) with geomechanics. (All vertical length scales 60 m; all radial length scales 500 m from central well.)

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

265

Fig. 21. The Changes in effective stress in the radial and vertical directions at the end of 780 and 2880 days, 3D model (seismic based areally non-uniform GH distribution) simulation with geomechanics. (All vertical length scales 60 m; all radial length scales 500 m from central well.)

and a later (t ¼ 2880 days) times behaviours. Note the effects of non-uniform (layer-by-layer) distribution of initial hydrate saturation on near-well volumetric strain changes are clear in these plots. Again the initial hydrate saturation layering continues to manifest its effects on volumetric strain. Cross-sectional plots of the vertical subsidence of considered next. As illustrated in the same Fig. 22, subsidence first occurs directly above the well. At the longest times, subsidence is seen to have a significant areal expansion throughout the reservoir (although maximum levels are still found directly above the well). At later times, substantial subsidence begins to occur further in out into the reservoir as the zone of hydrate melting expands.

8. Conclusions and future exploitation of Mallik gas hydrate resources 8.1. Conclusions The paper investigated two topics both individually and in combination. First we explored the use of well logs and seismic data to quantify the areally heterogeneous gas hydrate distribution. Reasonable correlations between the two wells of 2L-38 and 5L-38 (separated by 50 m) were found, especially in the two lower zones B and C. Seismic analysis showed that data was missing expected aquifer zones, however, probably because correlations with available well logs were not considered. In this paper, a method was presented to incorporate both available seismic and well log information into a Cartesian static model for the Mallik hydrate deposit for lower zone C. A revised estimate of hydrate saturation was provided for this zone C. The Cartesian static model was then converted to an approximately equivalent 3D radial model centred around well 2L-38. Consideration of the “tilted” and areally non-homogeneous nature of the hydrate deposit, observed from seismic, was incorporated into the model. The result was two limiting model scenarios of variable gas hydrate distributions that were used in dynamic simulations of the lower Mallik zone.

The second topic explored the role of geomechanics on long term hydrate production of the Mallik deposit, based on a mechanistic sensitivity study, described in Appendix A. A simplified radial model formed the basis of our treatment which allowed a comprehensive study of a variety of geomechanical factors via a sensitivity analysis. The model setup followed the radial grid choice utilized earlier for the non-geomechanics models, but with a coarser grid vertically, and with uniform properties provided in 5 zonal regions (a permafrost region overlying the upper/middle hydrate zones, and underlain by a hydrate free sand and aquifer). Geometric parameters obtained from logs (dynamic parameters) and rocks testing (static parameters) were both utilized in these sensitivity studies. The choice of dynamic mechanical parameters was equivalent to ignoring the geomechanical behaviour. Based on our simulations, the lower Mallik zone C can be expected to not reach plasticity limits, assuming reasonable mechanical parameters and a cased production well. This is contrary to earlier simulations in a vertically higher and less compressed reservoir produced via a horizontal well (Rutqvist et al., 2009). Ultimately, simulations of long time production for Scenario I and II were conducted without and with the role of geomechanics. A kinetic limitation for gas evolution limits the expected compaction drive effect to very long times. A vertical subsidence of approximately 0.17 m can be expected after 8 years of production. The conclusions from the present study confirm and extend our previous modelling studies. We now confidently reestablish the Mallik gas and water production trends simulated earlier. Gas production is divided into 4 distinct phases, which we feel is characteristic of gas hydrate dissociation occurring in a layered reservoir produced from a single vertical well. First, there is a rapidly rising rate period followed by an extended slowly rising rate. This is followed eventually by a significant period of very high rate gas production, finally terminating in another stabilized period of lower but significant gas production. With geomechanics, gas production in general and the gas production peak is shifted and delayed. The geomechanics effect is not purely compaction drive (as in conventional reservoirs, gas production increases with geomechanics). Based on the available seismic information, two areally

266

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

Fig. 22. Volumetric strains and vertical displacements at the end of 780 and 2880 days, 3D model (seismic based areally non-uniform GH distribution) simulation with geomechanics. (All vertical length scales 60 m; all radial length scales 500 m from central well.)

variable (i.e. non-homogeneous) scenarios for gas hydrate distributions were considered e Scenario I having the same initial total hydrate amount as our earlier (areally) homogeneous models, and Scenario II with significantly less overall total hydrate, but honouring the same relative distribution. In Scenario I, we observed multiple gas production peaks (which quite similar to 6 days production behaviour) with higher localized pressure pulses occurred due to strong reservoir gas hydrate heterogeneity. In Scenario II, this drastic change in gas production rate was not observed (due to faster pressure evolution in the reservoir). In both Scenarios I & II, the overall reservoir gas production peak is delayed compared to the homogeneous case. This is further delayed by the role of geomechanics. In all cases (2D homogeneous, Scenario I and II, with or without geomechanics), we observed a very similar overall production trend. This is probably a unique for the Mallik gas hydrate production using single vertical well. 8.2. Future work Future field simulations, especially considering the full Cartesian geometry, would benefit from intensive discussions between the seismic and production groups. Geomechanics responses in upper and middle zones (closer to the permafrost and nearer to plasticity limits) should be explored. Also the role of horizontal or slanted wells needs further assessments. Because future simulation involves larger areas of investigation and possible repeating well pattern development, a large volume of unprocessed available 3D seismic data, well logs from nearby hydrate wells (L-38, 3L-38, 4L-38) could be employed to get a more complete statistical representation of the interwell vertical stratigraphic variation in reservoir properties. This will allow quantification of gas hydrate accumulations and porosity variations. This could also provide a protocol in the use of 4D seismic methods to monitor gas hydrate dissociation spatially during long term production, analogous to that used to quantify steam zone evolution in heavy oil reservoirs (Zou et al., 2006). Further multi-scale analysis (laboratory work and MD simulation) of the kinetic model of gas exsolution seems warranted,

especially considering its impact on longer term geomechanics effects of compaction drive. Large scale laboratory models should also be investigated as a bridge between fundamental studies and field level experience. Acknowledgements The authors acknowledge the funding support from the Geological Survey of Canada (GSC) and the Alberta Innovates e Technology Futures (Project Number e 6658, HOOS 13255-2010). We thank G. Bellefleur and M. Riedel from Geological Survey of Canada for their support in collecting seismic data. Appendix A. Geomechanical sensitivity studies e gas hydrate dissociation by depressurization This appendix A investigates the impact of geomechanics on gas hydrate dissociation and production at the field scale, especially at longer production times. Here a simplified radial model forms the basis of our treatment which allows a comprehensive study of a variety of factors via a sensitivity analysis. The conclusions obtained here is applied in this study to a more targeted analysis of the role of geomechanics on gas hydrate production in a well-defined model of the Mallik lower (Zone C) gas hydrate deposit. A.1. Setup geomechanical simulation domain The model setup follows the radial grid choice utilized earlier for the non-geomechanic models, with an outer radius of 500 m. Vertically, however, a coarser grid was selected, with uniform properties provided in 5 zonal regions. From the top, a permafrost region down to 640 m, overlying the upper/middle hydrate zones down to 1158 m, below which the lower gas hydrate zone, consisting of 28 m of hydrate-free sand, 20 m of hydrate deposit, and underlain by 14 m of hydrate free sand and aquifer. These regions are illustrated schematically in Fig. A.1. A subset of this model, eliminating the permafrost and upper/middle hydrate zones was also utilized for many of the simulations to speed up run times. As

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

267

such, attention must be paid to the choice of geomechanics boundary condition at the top of the model to represent the effects

Fig. A.1. A schematic view showing simulation domain for the Mallik lower GH zone to study the geomechanical responses induced by GH dissociation process. Total 12 grid blocks were chosen near the wellbore to evaluate the failure mode during dissociation process.

of the eliminated overlying zones. This latter model will be reminiscent of the model chosen in this study to represent the detailed model of the lower Mallik GH region. The choice of layer-by-layer properties for this idealized geomechanics model results in someone different volumes in place from the base case models used in our previous studies of the Mallik lower GH zone (Uddin et al., 2012). These differences are summarized in Table A.1. A.2. Geomechanical data Following Rutqvist et al. (2009) two geomechanics parameter sets are chosen for our investigations, based on mechanical properties from laboratory geomechanics testing (static properties), versus properties derived from sonic well logs (dynamics properties). The presence or absence of hydrate can also be expected to affect the mechanical properties. The chosen values are listed in Table A.2. Essentially, the dynamic moduli represent “harder” material than their static counterparts, and the presence of hydrate adds further “hardness” in each case. In the present geomechanical simulation, the initial values of the stresses on finite element must be defined. Fig. A.2 shows the convention of stresses on a block. Here, sx, sy, sz are the normal stresses perpendicular to the x-, y-, and z axis, respectively and sxy, szx, szy are shear stresses on the surfaces perpendicular to the x-, y-, and z axis, respectively. In addition

Fig. A.2. Schematic of a local grid block showing the convention of the normal and shear stresses.

shear stresses are symmetry, therefore, sxy ¼ syx; szx ¼ sxz;

szy ¼ syz.

Finally, the initial stress state is somewhat uncertain, so sensitivity to this initial stress state was conducted utilizing several reasonable choices, listed in Table A.3. It is expected that smaller initial stress state ratios might bring stress paths closer to plasticity limits so this will be investigated as part of these simulations.

268

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

Table A.1 Initial material in place. GH reservoir

2D base case model (Uddin et al. 2012)

2D idealized model (geomechanical sensitivity studies)

Material in-place Formation pore volume (m3) Hydrate volume (m3) Water volume (m3, STD)

1.493  107 5.399  106 9.33  106

8.0  106 3.692  106 4.307  106

PeT conditions Average pressure (kPa) Average temperature ( C)

11,005 11.38

11,005 11.38

Table A.2 Basic geomechanical parameters. Geomechanical properties

Static properties (data set 1)

Dynamic properties (data set 2)

Silt/clay layers Young's E. modulus (kPa) Unload bulk modulus (kPa) Poisson ratio (e) Cohesion (kPa) Gamma (kPa) Exponent, n (e)

0.50e6 0.0 0.15 6.9  105 4.22  106 0.40

5.0e6 0.0 0.40 69.0  105 4.22  106 0.40

Sandy layers (Sh ¼ 0%) Young's E. modulus (kPa) Unload bulk modulus (kPa) Poisson ratio (e) Cohesion (kPa) Gamma (kPa) Exponent, n (e) Friction angle ( ) Dilation angle ( )

0.5  106 3.16  104 0.15 0.5  103 4.22  106 0.40 30 10

5.0  106 3.16  104 0.40 5.0  1013 4.22  106 0.40 30 10

Sandy layers (Sh ¼ 100%) Young's E. modulus (kPa) Unload bulk modulus (kPa) Poisson ratio (e) Cohesion (kPa) Gamma (kPa) Exponent, n (e) Friction angle ( ) Dilation angle ( )

1.8  106 3.16  104 0.15 2.0  103 4.22  106 0.40 30 10

18  106 3.16  104 0.40 20.0  1013 4.22  106 0.40 30 10

Fig. A.3. The effect of geomechanics on (a) gas production (rate and cumulative), and (b) overall field gas hydrate deposit and pore volumes e depressurization simulations over the length of 8 years prediction with geomechanics (based on the static and dynamics parameters) and without geomechanics simulations. No noticeable differences between the dynamic geomechanics and without geomechanics simulations.

Table A.3 Initial stress at the first grid cell (1,1,1) for geomechanical simulation runs (stress gradient option used to define initial stresses for all other grid cells). Initial stress (kPa)

szz sxx syy sxy syz sxz

Static parameter (set 1)

Dynamic parameter (set 2)

Run 1

Run 2

Run 3

Run 4

Run 5

Run 6

Run 7

Run 8

20,817.0

20,817.0 0.50szz 0.50szz 0.0 0.0 0.0

20,817.0 0.67szz 0.67szz 0.0 0.0 0.0

20,817.0 0.75szz 0.75szz 0.0 0.0 0.0

20,817.0

20,817.0 0.50szz 0.50szz 0.0 0.0 0.0

20,817.0 0.67szz 0.67szz 0.0 0.0 0.0

20,817.0 0.75szz 0.75szz 0.0 0.0 0.0

szz szz 0.0 0.0 0.0

A.3. Effect of geomechanics on GH production First, the role of geomechanics on gas hydrate production is examined, comparing the roles of dynamic versus static geomechanics parameters. These simulations utilize an initial stress state ratio of 0.67, although results basically turn out to be independent of this choice.

szz szz 0.0 0.0 0.0

Fig. A.3 shows gas rates and cumulative gas predicted for the static versus dynamic mechanical properties. This shows that the “softer material” (static properties) results in delayed gas production response although cumulative gas eventually reaches equivalent levels. Fig. A.3 also shows the variation in hydrate saturation and generated pore volume as the hydrate melting occurs. The harder (dynamic) material melts more quickly.

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

A.4. Geomechanical responses The spatial distributions of stress, volumetric strain, and vertical subsidence over the length of 8 year simulations were plotted and studied (Uddin, 2013). The results can be summarized as: 1) The stress distributions predicted for the static moduli case, with the largest magnitudes in effective stress occurring near the open well intervals. The stress distributions predicted for largest times simulated, showing a wide regional change in effective mean stress by these times. 2) Analogously, cross-sectional plots of volumetric strain distributions are studied next. The effects of non-uniform (layer-bylayer) distribution of initial hydrate saturation on near-well volumetric strain changes are clearly noticeable in the crosssectional plots. The initial hydrate saturation layering continues to manifest its effects on volumetric strain on later time distributions. 3) Subsidence first occurs directly above the well. At the longest times, subsidence is seen to have a significant areal expansion throughout the reservoir (although maximum levels are still found directly above the well). At later times, substantial subsidence begins to occur further in out into the reservoir as the zone of hydrate melting expands. A.5. Material failure assessment The sensitivity studies demonstrate that the lower Mallik zone C can be expected to not reach plasticity limits, assuming reasonable mechanical parameters and a cased production well. The failure criterion is discussed in this section.

269

The Coulomb failure criterion in a porous material is given (Bratli and Risnes, 1981; Risnes et al., 1982) by: s1  p ¼ 2S0 tan a þ (s3  p) tan2 a, here, p is the fluid pressure, s1 is the greatest principal stress, s3 is the smallest principal stress, S0 is the shear strength and a is the failure angle(p/4 þ ∅/2), here, ∅ is the internal friction angle. At the borehole surface, the criterion becomes, s1  p ¼ 2S0 tan a. In poorly consolidated sand, S0 will be less than 100 kPa and a will have values of about 60e65 . The right hand side of the above equation thus will have values less than 430 kPa. Under normal conditions, the left side will be much greater than this, regardless of which stress components in cylindrical coordinates (sr, s, sz) is the greatest. This means that critical stress values are reached near the borehole surface and the elastic solutions will not be valid in this region. To further illustrate the geomechanical behaviour occurring during the depressurization process, stress path plots (maximum effective stress versus minimum effective stress) are shown for selected locations in the reservoir. Fig. A.4 compares behaviour in layer 7 (near the top of the production zone) and layer 19 (near the middle of the production zone) for the static mechanical property choice. These paths are characterized by a generally increasing value of effective stress, first primarily in minimum effective stress but then followed by a distinct change to increasing maximum stress direction. This corresponds to the melting of the hydrate. Fig. A.4 also shows the analogous stress path results in the same two layers, but employing the dynamic mechanical properties. Here the shift in stress path direction occurs much earlier, and the overall behaviour is a relatively equal increase in minimum and maximum effective stress throughout the depressurization. Note these curves for different locations in one layer are often superimposable, as the time behaviour is not explicitly shown. Locations closer to the well should show stress path behaviour changes earlier in the simulation.

Fig. A.4. The geomechanical simulation with the i) static and ii) dynamic parameters showing the stress paths at 4 grid cells in two layers near the wellbore, (a) stress paths at layer 7 and (b) stress path at layer 19. The Coulomb failure lines are also shown in these figures.

270

M. Uddin et al. / Journal of Natural Gas Science and Engineering 20 (2014) 250e270

References Bellefleur, G., Riedel, M., Huang, J., Milkereit, B., Brent, T., 2010, Oct. 19e21. Towards Seismic Detection and Characterization of Gas Hydrate Accumulations in Permafrost Environment: an Example from the Mallik Gas Hydrate Field, NWT, Canada. CSUG/SPE 137499, presented at the Canadian Unconventional Resources and International Petroleum Conference, Calgary, AB, Canada. Bratli, R.K., Risnes, R., 1981. Stability and failure of sand arches. Soc. Pet. Eng. J., 236e248. Dallimore, S.R., Collett, T.S., 2005. In: Dallimore, S.R., Collett, T.S. (Eds.), Scientific Results from the Mallik 2002 Gas Hydrate Production Research Well Program, Mackenzie Delta, Northwest Territories, Canada, Geological Survey of Canada, Bulletin 585. Dallimore, S.R., Collett, T.S., Uchida, T., 1999. In: Dallimore, S.R., Uchida, T., Collett, T.S. (Eds.), Scientific Results from JAPEX/JNOG/GSC Mallik 2L-38 Gas Hydrate Research Well, Mackenzie Delta, Northwest Territories, Canada, Geological Survey of Canada, Bulletin 544. Dallimore, S.R., Yamamoto, K., Wright, J.F., Bellefleur, G., 2012. In: Dallimore, S.R., Yamamoto, K., Wright, J.F., Bellefleur, G. (Eds.), Scientific Results from JOGMEC/ NRCan/Aurora Mallik 2007e2008 Gas Hydrate Production Research Well Program, Mackenzie Delta, Northwest Territories, Canada, Geological Survey of Canada, Bulletin 601. Kurihara, M., Sato, A., Funatsu, K., Ouchi, H., Yamamoto, K., Fujii, T., Numasawa, M., Masuda, Y., Narita, H., Dallimore, S.R., Wright, J.F., Ashford, D.I., 2012. Analysis of 2007 and 2008 gas hydrate production tests on the Aurora/JOGMEC/NRCan Mallik 2L-38 well through numerical simulation. In: Dallimore, S.R., Yamamoto, K., Wright, J.F., Bellefleur, G. (Eds.), Scientific Results from JOGMEC/ NRCan/Aurora Mallik 2007e2008 Gas Hydrate Production Research Well Program, Mackenzie Delta, Northwest Territories, Canada, Geological Survey of Canada, Bulletin 601, pp. 261e289. Lee, M.W., 2002. BioteGassmann theory for velocities of gas hydrate bearing sediments. Geophysics 67 (6), 1711. http://dx.doi.org/10.1190/1.1527072. Parker, J.C., Lenhard, R.J., Kuppusamy, T., 1987. A parametric model for constitutive properties governing multiphase flow in porous media. Water Resour. Res. 23 (4), 618e624.

Risnes, R., Bratli, R.K., Horsrud, P., 1982. Sand stresses around a wellbore. Soc. Pet. Eng. J., 883e898. Rutqvist, J., Moridis, G.J., Grover, T., Collett, T., 2009. Geomechanical response of permafrost-associated hydrate deposits to depressurization-induced gas production. J. Pet. Sci. Eng. vol. 67, 1e12. Uddin, M., Wright, F., Dallimore, S.R., Coombe, D., 2012. Gas hydrate production from the Mallik reservoir: numerical history matching and long-term production forecasting. In: Dallimore, S.R., Yamamoto, K., Wright, J.F., Bellefleur, G. (Eds.), Scientific Results from JOGMEC/NRCan/Aurora Mallik 2007e2008 Gas Hydrate Production Research Well Program, Mackenzie Delta, Northwest Territories, Canada, Geological Survey of Canada, Bulletin 601, pp. 261e289. Uddin, M., Wright, F., Coombe, D., 2011. Numerical study of gas evolution and transport behaviours in natural gas-hydrate reservoirs. J. Can. Pet. Technol. 50 (1). Uddin, M., Wright, F., Dallimore, S.R., Coombe, D. Gas hydrate dissociations in Mallik hydrate bearing zones A, B, and C by depressurization: effect of hydration number (cavity occupancy) in hydrate dissociation. Journal of Natural Gas Science and Engineering (submitted for publication). Uddin, M., 2013. Numerical Studies of Methane Production from Arctic Gas Hydrate Reservoirs: IV. GH Dissociation Simulation with Geomechanics. Report prepared for the Geological Survey of Canada, Natural Resources Canada. Van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892e898. Wright, F., Uddin, M., Dallimore, S.R., Coombe, D., July 17e21, 2011. Mechanisms of gas evolution and transport in a producing gas hydrate reservoir: an unconventional basis for successful history matching of observed production flow data. In: Proceedings of the 7th International Conference on Gas Hydrates (ICGH 2011). Edinburgh, Scotland, United Kingdom. Wright, J.F., Dallimore, S.R., Nixon, F.M., Duchesne, C., 2005. In situ stability of gas hydrate in reservoir sediments of the JAPEX/JNOC/GSC et al. Mallik 5L-38 gas hydrate production research well. In: Dallimore, S.R., Collett, T.S. (Eds.), Scientific Results from the Mallik Gas Hydrate Production Research Well Program, Mackenzie Delta, Northwest Territories, Canada, Geological Survey of Canada, Bulletin 585. Zou, Y., Bentley, L., Lines, L., Coombe, D., June 2006. Integration of seismic methods with reservoir simulation, Pikes Peak heavy-oil field, Saskatchewan. Lead. Edge, 674.