3D Temperature distribution and numerical modeling of heat transfers in an active fault zone: Eugene Island 330, Offshore Louisiana.


Any reference should be made as "Guerin, Gilles 2000. Acoustic and Thermal Characterization of Oil Migration, Gas Hydrates Formation and Silica Diagenesis, PhD. Thesis, Columbia University"
 

4.1 Abstract

    At the center of an active growth fault system underlain by high-relief salt diapirs, the Eugene Island 330 field is a case study for the dynamics of active oil fields. Here, Plio- Pleistocene sandstone reservoirs are supplied with mature hydrocarbons by migration of fluids from overpressured shales upwards along the active fault system. The heat carried by the ascending fluids and the presence of highly conductive salt diapirs create strong temperature anomalies affecting the entire mini-basin. We present a method to interpret such complex temperature distributions by using some of the most commonly recorded data in modern oil field exploration: Bottom Hole Temperatures (BHT) and 3D seismic interpretation. More than 600 BHT from 200 wells allowed us to determine the present-day temperature field in a 15 ´ 15 ´ 6 km datacube. Using a lithological model built from the main seismic reflectors interpreted from 3D surveys, we perform a 3D numerical simulation of the conductive regime in the mini-basin. The results identify a 5-8° C temperature anomaly over the main salt diapirs. After subtraction of these results from the observed temperature field, we isolate a residual 5-10° C temperature anomaly potentially created by the upward migration of warm fluids along the multiple fault systems. Finally, using the fault surfaces interpreted from the 3D seismic survey to define the migration pathways, we build a 3D numerical model of the heat and fluid transfers required to reproduce these anomalies. Our simulations show that they could be created by periods of fluid expulsions lasting 1,000 to 5,000 years within the fault system. Such fluid expulsion activity is likely to have filled the most productive reservoirs of the EI330 oil field.

4.2 Introduction

    The Eugene Island 330 (EI330) oil field in the Northern Gulf of Mexico (Figure 4.1) has been the focus of an unprecedented interdisciplinary effort to understand the dynamics of an active oil field. The largest Pleistocene oil field in the world [Anderson et al., 1994], EI330 is at the core of an active growth fault system. Over the 30 years since discovery, multiple evidences have shown that hydrocarbons have been actively migrating along the fault during the recent history of the field and are still flowing: 1) the presence of hydrocarbon seeps at the sea floor, 2) discrepancies between oil maturity in the multiple reservoirs, and 3) temporal changes in the composition of oils produced over the last 20 years [Anderson et al., 1994, Holland et al, 1990]. Identifying the mechanisms and pathways of the hydrocarbon migration could indicate where to drill to take advantage of natural expulsion processes.

Figure 4.1

    To study this oil field, a 15 km ´ 15 km ´ 6 km `data cube' was assembled, combining multiple 3D seismic surveys, fluid samples, core, log and production data from over 600 wells. Region growing techniques were applied to High seismic Amplitude Events (HAE) in the 3D seismic data to identify the trails left by hydrocarbons moving along migration pathways [Anderson et al., 1994]. However, seismic imaging alone cannot be used to determine the volumes, duration and timing of hydrocarbon migrations. The comparison of successive 3D surveys indicates changes in seismic attributes within reservoirs under production, but this 4D analysis detects changes in seismic attributes due to forced fluid substitution and pressure depletion during the production of reservoirs [He, 1996, Chapter 3]. It is mostly sensitive to the movement of strong seismic events such as oil/water and gas/oil contacts on the scale of a reservoir [Anderson et al., 1994]. It cannot be used for the much slower migrations within the fault that feed the reservoirs over thousands of years.

    Unlike seismic attributes, formation temperature can be an indicator of long term fluid circulation. The upward migration of large volumes of fluid coming from deeper and warmer formations along localized conduits can create significant thermal anomalies [Hopper, 1990, Borner and Sharp, 1988]. The amplitude and the extent of the anomaly is dependant on both velocity and duration of the fluid migration, and it can remain detectable well after the end of the migration because of the low thermal diffusivity of most sediments. However, in addition to possible fluid migration along the fault, the presence of heterogeneous structures such as salt diapirs common in the Gulf of Mexico can also create strong thermal anomalies. Therefore, the interpretation of the temperature distribution requires to identify and quantify the various components of the thermal regime in order to isolate the specific contribution of fluid migrations.

    More than 600 Bottom Hole Temperatures (BHT) have been recorded in 200 wells penetrating the EI330 field. After correcting the data for drilling disturbances, we present the current 3-D temperature distribution in EI330, which displays several anomalies indicative of a complex thermal regime. To identify its various components, we numerically simulate first the 3D conductive regime to determine the temperature anomalies produced by the salt diapirs underneath the mini-basin. Assuming that the differences between this model and the measured temperatures result from advective heat flow through the fault, we reproduce these differences by a complete 3D numerical modeling of heat and fluid flows in the fault. Our results indicate the fluid migration history and mechanisms that could be at the origin of the observed temperature distribution.

4.3 Geological Setting

    The Eugene Island Block 330 field lies under 85 m of water, approximately 270 km southwest of New Orleans, near the southern edge of the Louisiana Outer Continental Shelf (Figure 4.1). The extremely rapid influx of Upper Cretaceous and Lower Tertiary terrigeneous sediments from the ancestral Mississippi delta has mobilized the Jurassic Louann salt, which left a trail of isolated diapirs across the shelf, and formed semi-continuous salt uplifts under the slope of the Gulf Coast margin [Woodbury et al., 1973]. The study area is made of nine property blocks located at the transition between these two salt provinces [Holland et al., 1990], with two salt diapirs underneath the South-East and the North-West corners of the study area (Figure 4.2).

Figure 4.2 

    The Tertiary and Quaternary sediments overlying the salt define three main facies representative of the normal evolution of a deltaic system prograding across a continental margin [Selley, 1988; Holland et al., 1990, Alexander and Flemings, 1995]: 1) directly over the salt, massive shales and turbidites were deposited in a prodelta environment. 2) They are overlain by a sequence of proximal deltaic sands and transgressive shales deposited during sea-level fluctuations when the delta slope was located nearby. The most productive reservoirs of the field are in this interval [Holland et al., 1990]. 3) The uppermost section was deposited after the delta had prograded southward, and is composed of fluvial massive sands. The transitions between the different sedimentation phases were identified in the 3D seismic data by reflectors corresponding to transgressive episodes. Biostratigraphic markers Cristellaria "S" (Cris S) and Small Gephyrocapsa (2) (Sm. Gep(2)) [Alexander and Flemings, 1995] define the bottom and the top of the proximal deltaic phase, respectively (Figure 4.3).

Figure 4.3

    Structurally, the EI330 field is a classic salt-withdrawal shelf minibasin [Alexander and Flemings, 1995]. It is bounded by two main fault zones (see Figure 4.3). The northern and eastern boundaries are defined by the listric, normal, "Red" Fault, while the south and west borders are defined by the counter-regional "Blue" Fault that developed as extensional compensation during the salt withdrawal to the south. Oil and gas reservoirs are trapped under two rollover anticlines on the downthrown side of the Red Fault.

4.4 The Temperature distribution in Eugene Island 330

4.4.1 Data Correction

    As one of the most productive oil fields in the Gulf of Mexico, EI330 has been extensively studied and sampled. Over 600 Bottom Hole Temperatures were collected from about 200 wells in the 9 blocks study area. During the drilling of a well, large quantities of "mud" are circulated in the borehole to facilitate the drilling, evacuate the cuttings and stabilize the hole. The influence of this circulation on the borehole temperature was first discussed by Bullard [1947]. Typically, the temperature of the mud flowing down the drill-pipe is close to sea surface temperature and lower than the formation temperature at depth. Therefore the deepest part of the well gets cooled, while the mud ascending outside the pipe carries some of the heat from the deep formation and increases the temperature in the upper section of the well. The time elapsed between the end of the mud circulations and the BHT measurement is usually not long enough for the borehole temperature to return to the formation temperature [Bullard, 1947]. It is however possible to estimate the influence of mud circulation on the borehole temperature and to determine the actual formation temperature from these data.

    Brigaud [1989] reviewed various methods to estimate the equilibrium temperature from BHT's. The most commonly used, often referred as the `Horner Plot', was originally theorized by Bullard [1947], and later developed by Horner [1951] and Lachenbruch and Brewer [1959]. It assumes that the effect of drilling and mud circulation is equivalent to the effect of a linear heat source in the formation, whose strength (S) is independent of depth and time. It also assumes that the formation is homogeneous and infinite. With these hypotheses, the formation temperature (Teq) can be expressed as a function of the measured temperature (TBHT), of the time elapsed since the end of the circulation (t), and of the duration of the circulation (tc):

  (1)

where Q is the heat influx per unit time and unit length of borehole wall during drilling and circulation and K the thermal diffusivity of the formation. In this form, Eq. (1) is difficult to use because the strength of the source (S = Q/4pK) depends on poorly constrained parameters, including pumping rate, mud temperature and thermal properties of the mud and formation. Two temperatures recorded in a well at the same depth but at different times with no thermal disturbances in between can help overcome this difficulty, by providing two equations for the two unknowns S and T. By reworking (1) for two BHT's T1 and T2 measured at times t1 and t2, the formation temperature at thermal equilibrium can be estimated by:

.    (2)

Even if the assumptions made are extreme simplifications of the influence of the drilling process, Lachenbruch and Brewer [1959] show that the maximum relative error is less than tc/[6(tc + t)], which for typical values of tc = 2 hours and t = 6 hours is about 4%.
    This method was used to correct about 25% of the measurements. For the other data where only one measurement was available at each depth, the impossibility to determine S in (1) required an other approach. We used a method described by Gable [1977, 1986], who shows that within a confined area, the single-measurement BHT's can be corrected by a polynomial function of depth Tcorr(z) defined from the actual Horner corrections made in this area:

Teq = TBHT + Tcorr(z). (3)
    To determine Tcorr(z) for the study area, we calculated a 2nd order polynomial fit for the data already corrected by the Horner Plot method (t(z) = 26.94 + 0.0249z - 1.748´ 10-07z2), and for the entire set of uncorrected BHT's (tBHT(z) = 29.49 + 0.0171z - 3.714´ 10-08z2). The correction law Tcorr(z) is the difference between these two functions:
Tcorr(z) = t(z) - tBHT(z) = -3.45 + 0.0078z - 1.377´ 10-07z2. (4)
The negative constant coefficient (-3.45) illustrates the warming of the upper section of the hole because of the mud circulation.

    This method does not take into account differences between separate wells in the duration of the mud circulation, or in the time elapsed between the end of the circulation and the measurement. It implicitly assumes that the drilling conditions, the timing of operations and the drilling mud were comparable for all the wells in the area, and that the thermal diffusivity of the formations surrounding the different wells were also similar. These assumptions are difficult to verify but reasonable since all the wells in this small area encountered the same type of lithology and were drilled with the same techniques.

Figure 4.4

    Figure 4.4 shows the entire set of BHT's before and after correction and the two polynomial fits t(z) and tBHT(z). Two temperature logs measured in separate wells (A-20 and C-19) in this area are also shown for comparison and illustrate the evolution of the temperature in a borehole. Both logs were recorded a short time after the end of circulations. The log in C-19 was recorded after a shorter time than in A-20, and displays lower temperatures at depth, but both logs are close to the polynomial fit for the uncorrected BHT's and within the range of these data.

Figure 4.5

    In the interval with the highest density of data (1750m-3000m, see Figure 4.4), BHT's were available at several depths in most wells. We calculated the temperature at specific depths in these wells (i.e. 1750, 2000, 2500 and 2750 m) by linear interpolation or extrapolation between the two most closely located corrected measurements. The calculated temperature distribution is shown in Figure 4.5. Using a regular mesh of 50 ´ 50 elements for consistency with the following modeling, these surfaces were produced by a gaussian kernel smoothing of bandwidth 5.

4.4.2 Heat transfer in sedimentary basins

    To interpret the temperature distribution in Figure 4.5, we will summarize the constituents of the thermal regime in sedimentary basins, before considering the specific features of EI330. The thermal regime in sedimentary formations can be described by a general expression of the conservation of energy in porous media :

(5)

where V and S are respectively the volume and the external surface of an arbitrary unit of sediments and  is the outbond normal to S. The left side is the variation of internal energy, a function of the bulk density (r) and specific heat (C) of the materials and of the temperature (T). The first term on the right represents the heat transfer by conduction defined by Fourier's law and proportional to the thermal conductivity l of the material and to the temperature gradient () across its surface S. This term is usually the dominant heat transfer component in sedimentary basins [Rybach, 1981] and can generate strong anomalies in the presence of geological bodies of high thermal conductivity, such as salt diapirs [Mello et al., 1994] or ore deposits [Castoviejo et al., 1997]. The second term represents the heat carried by fluids of specific heat Cf and density rf flowing with Darcy velocity . We assume here that the pore space is fluid-saturated, and that the volumetric flow rate can be expressed by Darcy's law for single phase fluid flows = -k/m[ - rf], where is the pressure gradient, k the permeability of the sediments, m the dynamic viscosity of the fluid and the acceleration of gravity. The last term in Eq. (5) is the heat generated by eventual sources of volumetric strength H.

    For the numerical formulation, assuming that the temperature distribution is continuous and varies smoothly on the scale of the grid blocks, it is more convenient to rewrite (5) as a local partial differential equation:

    (6)

4.4.3 Heat sources

    In most sedimentary formations, possible heat sources are the radiogenic heat produced by the decay of radioactive isotopes of Potassium, Thorium and Uranium, and the heat produced by chemical reactions such as diagenesis or hydrocarbon generation.

4.4.3.1 Radiogenic heat

    Values collected by Rybach [1981] show that in clastic sediments the volumetric heat production due to radioactive decay is about 1.0 mW/m3. If the specific heat and the bulk density of sediments are respectively 2000 Jkg-1K-1 and 2500 kg/m3, the time required to increase the temperature of a unit mass of sediments by 1fC would be 150,000 years, which is much longer than the periods of activity in the Red Fault, lasting a few thousand years [Roberts et al., 1996, Anderson et al, 1994]. Moreover the globally uniform lithology of EI330 makes radiogenic heat a practically uniform process over the entire area, that cannot be associated with significant localized anomalies such as seen in Figure 4.5. Because of its numerical simplicity, this term is however included in the following simulations to remain as close as possible to the actual thermal regime.

4.4.3.2 Heat generation from oil maturation

    Most thermal regime studies do not consider the generation of hydrocarbons as a possible source of heat [Rybach, 1981, Bachu et al, 1995]. Studies using BHTs to constrain temperature distributions are usually devoted to larger areas and longer histories, and use more sparsely distributed data sets [Bodner and Sharp, 1988, Bachu et al., 1995, Armstrong et al., 1996]. This limit their scope to the large scale dominant components: conduction and advection, as well as radiogenic heat for large time scales. By contrast, our goal was a comprehensive understanding of the temperature distribution in a relatively smaller area using a higher density temperature data set so that we could study short-lived time-dependent processes. Therefore, before discarding any possible component to simplify the formulation, we have to estimate the strength of most possible components, including the heat generated by the transformation of kerogen, which is mostly exothermic. While a complete account of the thermochemistry within the source rocks is beyond the scope of this study, we can estimate an order of magnitude of the heat released in the source rocks. Since the reactions are multiple and complex, we use two different estimations of this process: an average value over time, and its intensity at its peak productivity.

England et al. [1987] estimates that a typical source rock will yield about 20g of hydrocarbon per kg of sediments or 48 kg/m3 for an average sediment density of 2400 kg/m3. This release occurs mostly between 120 and 150° C [England et al., 1987], which in reasonable heating conditions (1 to 10° C/ Ma) corresponds to an average generation rate of 5.0 10-11 to 5.0 10-10 g/m3/s. Assuming an average molecular weight of 150 g/mol (~ C15H32, , which is the median component of the chromatogram measured by Wheelan et al, [1994]) , and an enthalpy of reaction of 200 kJ/mol (based on Mallard and Linstrom, [1998]) the heat generation rate is 6.710-2 to 6.7 10-1mW/m3, lower than the radiogenic heat that we already discarded.

    This value represents an average over several million years, underestimating the heat generated at the peak of hydrocarbon productivity [Tissot et al., 1987]. From the description of oil generation by Burnham [1995], it is possible to use chemical kinetics to estimate the heat generated at the peak of kerogen transformation. Oil generation from kerogen obeys an Arrhenius-type of law:

      (7)

where C is the oil concentration in the source rock, T the temperature (in K), R the perfect gas constant, A is the frequency factor in s-1, and E the activation energy of the reaction, in J/mol. Burnham [1995] gives values of 3.0 1013 s-1 for A and 200 kJ/mol for the activation energy for type II oils (marine sediments). For the same concentration of organic matter as used previously (20g/kg of rock, or 320 mol/m3 using the same values for the density and molecular weights as before), and at temperatures of maximum generation (T = 130° C, from Tissot et al., [1987]), the rate of kerogen transformation is dC/dt= 7.010-11 mol/s. With the same enthalpy of reaction used earlier (200 kJ/mol) the heat generation rate is of the order of 14mW/m3 at its peak.

    Both estimates are of the order of the radiogenic heat (1 mW/m3), already discussed and discarded. Moreover, the most productive source rocks are located in strata at 120 to 150° C, which are the deepest layers of our model. We will therefore ignore the influence of oil generation on the shallower anomalies of Figure 4.5.

4.4.4 The effect of sedimentation

    In the following interpretation and modeling, we will assume that the conductive regime is permanent. This simplification neglects the effect of sedimentation and salt movement on the temperature distribution. While the salt movement beneath EI330 stopped about 1.5 Ma ago, sedimentation is an active process occurring at an approximate rate of 1.0 mm/year, or 1km/Ma [Alexander and Flemings, 1995]. We have to estimate the magnitude of its influence on the thermal regime before any further simplification. A crude comparison of the sedimentation and heat diffusion time scales can be derived from Eq. (6). The time Dt required for heat to diffuse through a thickness D of newly deposited sediments can be estimated by the following approximation:

    (8)
where k (= l/rC) is called the thermal diffusivity of the sediments. For 1 km of newly deposited sediments of typical diffusivity 10-6-10-7 m2/s, the time required for heat diffusion through the layer is 1012 -1013 s or 0.03 - 0.3 Ma, which is comparable with the deposition time (1 Ma at 1mm/year). Because sediments are deposited continuously over time, this estimation exaggerates the slowness of heat diffusion, whose reach increases only with the square root of time ( D a (Dt)1/2). While one meter of sediments gets deposited in 1,000 years, Eq. (8) shows that it takes less than a year for heat to diffuse across this distance. These values mostly underline how each process has to be weighted in terms of time scale and distance for a correct interpretation.

To estimate more accurately the effect of sedimentation on the conductive regime we can solve the 1D equation of heat conduction in a homogeneous medium moving with a velocity U equal to the sedimentation rate [Jaeger, 1965]:

.     (9)
Carslaw and Jaeger [1959] and Jaeger [1965] use Laplace transform to solve (9). In an infinite medium with a constant sedimentation rate, a linear temperature profile at the origin and a constant temperature at the surface, the solution is:

(10)

where erfc is the complementary error function: erfc(x) = 1 - (2/^¼)òox e-u2du, go is the initial temperature gradient and To the temperature at the surface.

Figure 4.6 

    Figure 4.6a shows the evolution of the temperature profile vs. depth (dotted lines) through time over 6 Million years, equivalent to the deposition of our data cube, and Figure 4.6b shows the evolution of the temperature at the same depths as the surfaces of Figure 4.5. Both figures indicate a steady decrease of temperature with time, at a rate of about 2.0 °C/Ma at the end, corresponding to the present day. This analytical solution overestimates the cooling effect of sedimentation because the infinite medium hypothesis does not allow the heat support provided by the regional basement heat flow. As a result, at any depth, (10) converges to the surface temperature T0, which is much more dramatic than the actual effect of sedimentation. However, Figure 4.6b shows that even in such extreme model, the temperature change over the last 100,000 years is less than 0.2°C. Without an exact analytical solution, a likely more reasonable estimate of the 1-D effect of sedimentation can be calculated by numerical method. A Forward-time Centered-Space discretization of (9) is:

     (11)
where Tin indicates the temperature at time step n at the node i of the discretized vertical column. Dt and Dz are the constant time step and uniform node spacing, respectively. An implicit Crank-Nicholson formulation of (11) expresses the temperature increment between successive time steps as the average of the increments at the two time steps. If (11) can be written by Tin+1 = Tin + f(Tin), this algorithm can be expressed as Tin+1 = 1/2.[f(Tin)+f(Tin+1)]. If we note a = UDt/2Dz and b = kDt/Dz2, the matrix form of this formulation is:
(12)
which comes to invert a tri-diagonal system of equations while imposing the following boundary conditions: constant seafloor temperature at the top and constant temperature gradient at the base of the sediment column. The results are also shown in Figure 4.6a and 4.6b (continuous lines) for a 6 km-thick sedimentary sequence deposited at 1mm/year with a 30°C/km linear temperature profile at the origin. They show in particular that the temperature is almost in steady state after about 1.5 Ma, and practically constant at the present time. The equation of the steady state profile can be derived by solving (9) without the time dependent component:
k2T/?z2 = U?T/?z     (13)
    The solution of this linear differential equation is:
    (14)
where D is the thickness of the sedimentary sequence and gb the temperature gradient at this depth, assuming a constant purely conductive heat flux. This solution, shown in green in Figure 4.6a with a 1 mm/year sedimentation rate and 30 °C/km gradient at depth, is barely distinguishable from the Crank-Nicholson solution after 1.5 Ma.

    These 1-D estimations of the influence of sedimentation can not be representative of the actual processes in EI 330 where sedimentation rates have varied in space and time and where sediment properties are not uniform, but they show that for periods of time of the order of 10,000 to 100,000 years the present conductive regime can be considered stationary.

4.4.5 Analysis of the present-day temperature distribution

    The most distinct features present at all depths in Figure 4.5 are the lower temperatures on the downthrown block in the South-West (Block 337), and the high temperature zone extending from the South-East corner (Block 339) to the North-West (Blocks 331, 314), across Block 330.

    In Figure 4.6a we also display the result of Eq. (14) at two different sedimentation rates, 1mm/year (green) and 2 mm/year (red) showing that faster sedimentation rates can decrease considerably temperature at depth. Therefore, the faster sedimentation rate in the down-thrown block [Alexander and Flemings, 1995] can contribute to the lower temperature observed in the SW in Figure 4.5. An other reason for these lower temperatures is the difference between the relative thickness of the sand-rich and shale-rich formations across the fault (see Figure 4.3). The proximal deltaic and fluvial sections are sand-rich and have a higher thermal conductivity than the underlying shale. Assuming a uniform vertical heat flow from the basement, the thick sands maintain a lower temperature gradient over a greater depth in the downthrown block (~4000 mbsf) than in the upthrown section where shale becomes dominant below only 2000 m.

    The higher temperatures along the SE-NW trend can be visually correlated in Figures 4.2 and 4.3 with the salt diapir system and with the general pattern of the Red Fault system. The high thermal conductivity of the salt compared to the surrounding sediments (4-6 Wm-1K-1 vs. 1-3 Wm-1K-1, [Mello et al., 1994]) and the high amplitude of salt topography, described by Anderson et al. [1994] as the most rugged topography on earth (see Figure 4.2), can generate significant temperature anomalies [Mello et al., 1994]. So could warm fluid ascending along the fault [Roberts et al., 1996]. The respective contributions of these two regimes have to be differentiated.

    Coelho et al. [1996] show that the presence of oil and gas decreases the thermal conductivity of reservoir sediments and generates a thermal blanketing that increases temperature underneath the reservoirs. They argue that this influence of hydrocarbons could create a positive temperature anomaly as high as 15° C below the main reservoirs and could explain by itself the observed anomalies. However, using their equations, we find that such blanketing would require about 60% of gas distributed uniformly in the pore space with an average porosity of 40 % between 1,000 and 2,000 mbsf. This is hardly representative of the stack of mostly oil-saturated low-porosity reservoirs constituting the Block 330 field, and could not explain the shallower anomalies. This blanketing must contribute to some of the observed anomalies, particularly in the deeper layers, but we choose to take it into account in the following models only by attributing lower conductivity to the sands of the downthrown block.

4.4.6 One dimensional analysis of heat advection

    The dominant terms in Eq. (5) are the conductive and the advective components, whose relative strength define the nature of the thermal system [Rybach, 1981]. The differentiation between the strength of these two terms can be performed on a first order by calculating the Peclet number of the system, a dimensionless number that can be calculated from (6) by

Pe = rfCfVD/l (15)
where D is a characteristic length for fluid migrations. For large values of Pe (Pe >>1) heat advection dominates, for low values (Pe <<1), conduction is preponderant, and for intermediate values the two modes have comparable influences. To determine the strength of migrations in the Red fault that could produce anomalies comparable to the influence of the salt, we can calculate the velocity required for the Peclet number to be equal to 1. A value for D can be found in Holland et al. [1991] who observe that migration must have occurred over more than 3650 m vertically in the recent history of the EI330 field. For this comparison, the strength of the advective regime has to be compared to the strongest conductive regime in the area, within the salt. Reasonable values for the salt conductivity and for the density and specific heat of the fluids are, respectively l = 6 Wm-1K-1, rf=1000 kg/m3 and Cf = 4.185.103 JK-1kg-3. Therefore, the velocity of the fluid migration necessary for heat advection to compete with the heat drainage through the salt is V= l/rfCfD ? 4.10-10 m/s, or ? 1.2 cm/year. Even though this is a crude one-dimensional comparison of complex processes, this value is not exceptional in an area of significant geopressures and high pressure gradients, and is considerably smaller than the maximum velocities (up to 100 m/yr) considered by Roberts et al. [1996], who modeled the expulsion mechanisms in the same area. It shows that the influence of the salt and of ascending fluids have to be both considered in the interpretation of Figure 4.5.

4.5 Numerical modeling of the conductive regime

4.5.1 Model description

    Assuming that sedimentation is uniform and at a constant rate, it is possible to model the combined effects of sedimentation and heat conduction in permanent regime. Using control volumes fixed with regard to the seafloor (the top of the sediments), sedimentation can be considered as a particular form of heat advection. In permanent regime, this can be written:

     (16)
where represents the sedimentation, assumed vertical and downward, and r and C are the density and specific heat of the sediments. For simplification, we consider here that all sediments have the same heat capacity (rC). We use a finite difference formulation to rewrite this equation in finite difference within a regular rectangular grid. For an element n of dimensions Dx,Dy,Dz, it can be rewritten

    (17)

where the subscripts refer to the six elements connected: (T)op, (B)ottom, (W)est, (E)ast , (N)orth, (S)outh, and dm and Sm are respectively the grid dimension along and the surface orthogonal to the m direction. The thermal conductivity between two adjacent elements is the harmonic "in series" average of the thermal conductivity of n and m :
    (18)
    The centered space over-relaxation method [Press, 1994] used to solve (17) is an iterative algorithm where an initial value for the temperature is attributed to all the nodes at the initial step (T0), and the temperature Ti+1 at step (i+1) for any node is calculated from the values of the previous step (i) at the adjacent nodes, until reaching a convergence threshold. In our rectangular uniform grid, where the node dimensions are h ´ h ´ v the recurrent formula is :

Tni+1 = (1 - a)Tni +
b[h2(TTikT + TBkB + g(TTi - TBi)) + v2(TNikN + TSikS + TWikW + TEikE)]    (19)
with 

where a is the over-relaxation coefficient, between 1.0 and 2.0 for unconditional convergence [Press, 1994].

    If the heat capacity of the sediments is uniform, the construction of this model requires only to specify the thermal conductivity. We use four different types of formation with specific attributes in our model, corresponding to the main facies associated with the development of the minibasin: the salt (1) is overlain by massive shales (2), underlying interbedded sands and shales (3), and the uppermost section is composed of massive sands (4). The salt topography and the transgressive surfaces Cris S and Sm. Gep (see Figure 4.3) used to delineate these units were directly imported and gridded from the 3D seismic interpretation. Reproducing each individual sand between these two surfaces would require a much finer grid and would not change significantly the results because of the limited difference in conductivity between sands and shales compared to the salt. Therefore, only a few representative sand units were used in the interbedded sand/shales sequence. This simplified lithology and structure are shown in Figure 4.7. Thermal conductivity values are given in Table 1.

Figure 4.7 

    The boundary conditions are chosen to be representative of a passive margin environment: 1) uniform temperature at the seafloor (18° C), 2) null heat flux across the lateral boundaries (null horizontal gradients) and 3) a uniform vertical heat flow at the bottom boundary (60 mW/m2, from Anderson et al [1991]. The over-relaxation algorithm was implemented on a regular 50 ¥ 50 ¥ 50 elements mesh, of dimensions 300m  300m  120m. The initial temperature distribution (T0) is defined by the seafloor temperature (18° C) and an uniform linear temperature gradient of 30° C/km. The convergence was considered attained when the actual heat influx in any element (the left side of Eq.(17)) was less than 100 W, or about 10 mW/m3.

4.5.2 Results


Figure 4.8 

Figure 4.9

    The results in Figure 4.8 indicate a ~10° C temperature anomaly following the SE-NW pattern of the salt (see Figure 4.2). The influence of salt domes in Block 314, 330-329 and 339 can be identified and correlated with similar features in Figure 4.5. Since we do not include in this simulation the effects of the differences in sedimentation rates between the two sides of the fault, nor the effect of the presence of oil and gas in the sediments of the hanging wall, we could not reproduce the generally lower temperatures on the downthrown side of the fault. Therefore, the features apparent in Figure 4.8 are mostly associated with the salt. The difference between these results and the observed temperatures (Figure 4.5) is shown in Figure 4.9 where a substantial ~10fC residual anomaly remains at all depths in the center of the study area. The superimposition of the fault pattern on this figure (white lines) shows that this residual anomaly follows roughly the SE-NW arc of the "Red" fault, from Block 339 to Blocks 331-314, through Block 330. In the deepest section, the anomaly is mainly located in the southern oil field (Blocks 338-339), while in the shallower surfaces, it is distributed along the arc of the fault zone and centered on the block 330 and 331 oil fields. The general coincidence of these residual anomalies with the regional fault pattern suggests that they could be produced by heat advection through fluid migrations along the active faults.

4.6 Numerical simulation of fluid migration

4.6.1 Indications for fluid circulation within the Red Fault

    The temperature anomalies along the Red Fault are only one manifestation of the fluid migration in the recent history of the Eugene Island area. Hydrocarbon seepages at the seafloor and discrepancies between the respective maturity of the oil and of the reservoirs [Holland et al., 1991], temporal changes in the composition of the oil produced over the 20-years history of the field [Wheelan et al., 1994] and the slow depletion rate of the production during this period [Anderson et al., 1993] were among the factors that indicated the ongoing migration and initiated the Global Basins Research Network interest in EI330. Region growing techniques applied to 3-D seismic amplitudes were used to identify the migration paths characterized by High Amplitude Events [Anderson et al, 1994] left by trapped migrating hydrocarbons. These results were used to determine the trajectory of the Pathfinder and A-10 wells within the fault zone, in which successive Thermal Decay Time (TDT) logs detected migration of fluids that were occurring at the time of the measurements [Losh et al., 1999, Anderson et al., 1994].

    Fluid migrations in EI330 are directly related to the presence of anomalously high pressures that are common to many areas of the Gulf Coast. While various mechanisms can generate geopressures, Hart et al. [1995], and Mello et al. [1994] show that compaction disequilibrium is the principal factor for their development in Eugene Island. Pore pressure increases anomalously when the deposition of low-permeability clay-rich sediments is too fast for the pore fluid to be expelled while the sediments get buried [Hart et al., 1995, Bredehoeft and Hanshaw, 1968]. The location and the depth of the geopressures in the Eugene Island area have been discussed by He [1996], Hart et al., [1995], and Mello et al. [1994]. Roberts et al. [1996] have shown that the periodic release of overpressured fluids could be at the origin of significant temperature anomalies.

4.6.2 3D Numerical simulation of the pressure and temperature

    If the residual anomalies in Figure 4.9 are produced by fluid migration along the fault, the numerical simulation of these anomalies can give an estimate of the duration, fluid volumes and velocities necessary to create them. Roberts et al. [1996] used a two dimensional finite element model to reproduce the observed temperatures. Their results give an indication of the duration and of the strength of the migration, but they are only based on a representative two dimensional cross section of the area. This does not account for the actual three-dimensional nature of the system shown by the distribution of anomalies in Figure 4.9. This figure suggests that the deepest part of the fault is most active in its southern extension (Blocks 338 and 339), and that while migrating upwards, the fluids simultaneously propagate to the North-West in the fault, and spill into the downthrown block, filling the reservoirs in Blocks 330, 331 and 314. That is, if the anomalies are actually generated by fluid migration, this figure shows that the migration is not only ascendant, but also lateral, and intrinsically three dimensional. Therefore, a realistic simulation has to consider all the spatial dimensions of this system, as well as the 4th dimension, time.

4.6.2.1 Numerical formulation

    To simulate time-dependent heat advection requires consideration of pressures in the data cube, in addition to temperature. The additional equation necessary is given by the law for the conservation of pore fluid mass in saturated porous media:

    (20)
with notations similar to (5). The term on the left represents the change of pore fluid mass over time, the first term on the right is the fluid transport through the control surface S, and the last term is an eventual fluid source of volumetric strength G. Once more, the local partial differential equation is the formulation actually used in the numerical discretization:
    (21)
    We use an Integrated Finite Differences Method (IFDM) for this discretization [Narasimhan and Witherspoon, 1976]. Our model was adapted from the code "PT" originally developed and described by Bodvarsson [1982]. Equations (5) and (21) are formulated as follows for a node n connected to an arbitrary number of nodes m :
+ (rfCfVdST)n,m + (HV)n   (22)

   (23)

    (24)

 whereb is the fluid compressibility (b = -1/V.?V/?P), Vdn,m is Darcy's velocity at the interface between nodes n and m, Tm,n is the temperature at this interface, hn,mthe direction cosine of its outward normal (for a vertical surface, hn,m= 0), and Dn,m and Dm,n are respectively the distances from the nodal points in n and m to their shared surface of area Sn,m.

4.6.2.2 Properties at the interface

    rrn,m and mn,m are respectively the fluid density and viscosity at the interface. Both are calculated as a weighted average of the values for each nodes. At each time step, fluid viscosity, density and compressibility are calculated for each node as a function of temperature and pressure after Batzle and Wang [1992].

Tn,m is calculated by an upstream weighting formulation :

Tn,m = uTu+ (1-u)Td    (25)
where u and d designate the upstream and downstream nodes and u is the upstream weighting coefficient, restricted in values to the range 0.5 to 1 for unconditional stability. We use u = 0.7.

    The permeability and the thermal conductivity between nodes are calculated with an harmonic mean, to insure continuity of fluxes at the interface:

   (26)
4.6.2.3 Sources

    In Equations (22) and (23), Hn and Gn are respectively the volumetric heat and mass generation rates of eventual sources in element n.

    Hn was discussed and discarded earlier.

    The possible volumetric mass generation rate Gn could include the generation of hydrocarbon, clay dehydratation, and fluid expulsion by compaction. England et al. [1987] estimate that in a typical oil producing environment, the velocity of fluids expelled by a source rock containing an average potential of 0.02kg of oil per kg of rock is of the order of 10-14 m/s or 3.2 10-7 m/year. This value is several orders of magnitude lower than the velocities estimated by Roberts et al. [1996] in EI330, and proves to be also much smaller than velocities in our own models. Therefore we do not account for this potential source or for similarly slow fluid generation processes (clay dehydratation) and we do not include any mass generation term in our simulations.

4.6.2.4 Lithology description

    The model grid was built from the same surfaces used in the conductive model: the top of the salt, and the Cris S and Sm. Gep 2 transgressive surfaces [Alexander and Flemings, 1995]. For the same reason as in the conductive regime simulation, only four representative sand units were defined in the sand/shale sequence, and the representation of the lithology was similar to the conductive model in Figure 4.7 (see Figure 4.10). Because the calculation of time-dependent pressure and temperature over the entire datacube requires the resolution of large, sparse linear equation systems at every time step, the computational constraints are higher than for the conductive regime simulation. The largest model for the entire data cube had 20 (EW) ´ 20 (NS) ´ 25 (vertical) = 10,000 elements, ranging in size from 1000´ 1000´ 1000 m to 500´ 500´ 100 m. Thermal conductivity values were the same as in the conductive model, and the permeability ranges are given in Table 1. The grid is shown in Figure 4.10.

Figure 4.10 

4.6.2.5 Numerical controls

    For stability purposes, the linear equations are solved implicitly. This means that at every time step, the values used for pressure and temperature in (22) and (23) are the values at the end of the time step. If Tni-1 and Pni-1are the temperature and pressure at the beginning of step (i), the values used for Tni and Pni at step (i) in (22) and (23) are Tni = Tni-1 + DTn and Pni = Pni-1 + DPn . The actual unknown in the equations are DTn and DPn. A sparse matrix solver, MA28 from the Harwell Subroutines Library, is used to solve the linear equations resulting from this formulation.

    While the implicit formulation is unconditionally stable, it is only first order accurate in time and requires time steps to be short enough to reproduce short-lived events. To limit this effect on our results, we limit the length of each time step so that the fluid displacement during one time step doesn't exceed the minimum dimension of any of the elements. From the previous 1-D analysis, velocities of the order of 4´ 10-10 m/s ~ 1.2 cm/yr must exist to generate anomalies competing with the salt. Roberts et al. [1996] simulate velocities of more than 100 m/yr, but create anomalies much stronger and more focused than what we observe in Figure 4.9. This suggests that intermediate values on the order of 1 m/yr can be used as a reasonable order of magnitude. The smallest element dimension in our grid being 100 m this limits the simulation time steps to less than Dtmax ~ 100 years, which we used as a the maximum time increment. Since the maximum velocity is actually not constant, it is calculated after each time step and the time increment is reduced if any "overshoot" was possible over any node. Tests made consequently by allowing longer time steps (up to 5,000 years) provided results similar to the results that we present.

    Initial temperature conditions are the results of the conductive model run on the same grid. Initial pressures are defined as a linear profile from seafloor pressure (1.1´ 106Pa) to a pressure at the bottom boundary equivalent to 90 % of the average lithostatic pressure, similar to the values measured in the Pathfinder well in the overpressured shales [Anderson et al., 1994].

    The boundary conditions are defined by (1) constant and uniform pressure and temperature at the seafloor, (2) constant pressure and steady, uniform vertical heat flow at the bottom, and (3) no heat or fluid flux at the lateral model boundaries, maintained by null horizontal temperature and pressure gradients. Values are given in Table 2.

    The fluid present in the pore space and migrating is assumed to be brine. Salinity is calculated as a linear function of depth derived from the compilation of Gulf Coast data of Batzle and Wang [1992]. The density, viscosity and compressibility of the fluid are calculated from salinity, pressure and temperature after Batzle and Wang [1992].

4.6.2.6 Fault simulation

    The fault zone is defined in the model by importing the main fault surfaces interpreted from 3D seismic: the down to the south "Red Fault", and the counter-regional down-to-the-north "Blue" and "F" faults [Alexander and Handschy, 1998, Rowan et al., 1998]. The fault zone is represented by a subset of variable-permeability "fault nodes", which are the elements of the grid that are crossed by the imported fault surfaces, and the elements directly adjacent in order to assure connectivity between diagonally connected nodes in the rectangular mesh (See Figure 4.10b). These nodes are given the same attributes (conductivity, density, specific heat, porosity, permeability) as the lithological unit (shale/sand) within which they are defined. The fault opening is simulated by increasing their permeability. Roberts et al. [1996] show that the fault must have increased permeability as deep as 2,000 meters below the top of overpressure to reproduce temperature anomalies of amplitudes equivalent to that observed in Eugene Island. Therefore the deepest nodes allowed to open in our models are located at about 5,000 mbsf.

    Wang and Xie [1998] provide a complete discussion on the criteria that can be used to open the fault by compaction-induced hydrofracturing, and on the permeability increase due to the fault opening. The most commonly accepted condition for fault opening in sediments, which we use in our model, is when the pore pressure exceeds the lowest compressive strength of the sediments, roughly 85% of the lithostatic pressure [Wang and Xie, 1998, Roberts et al., 1996, Hart et al., 1995]. When the pore pressure reaches this value in a fault node, its bulk permeability is increased by a factor of 102 to 103 [Wang and Xie, 1998].

    Without simulating the sedimentation and compaction disequilibrium that generate and maintain the overpressure, we could not reproduce the actual cycle of activity of the fault. However, by extending the results of Wang and Xie [1998] to the sedimentation rates of EI330, we can estimate that the shortest time between successive hydrofracturing events should be over 200,000 years, which leaves enough time for the 8 to 10° C anomalies to dissipate almost entirely between two events. Therefore we assume that the observed anomalies are the result of a single expulsion event and we monitor only the temperature produced by the event triggered by the overpressures defined in the initial conditions. The duration of the enhanced permeability within the fault is not imposed, but implicitly controlled by two parameters: 1) the value of the permeability increase in the fault nodes when the fault opens, which determines the time taken for the overpressured fluids to be expelled, and 2) the closing pressure, which is the pressure below which fault nodes are re-assigned their original permeability, simulating the end of the activity period. By introducing a `closing pressure', lower than the opening pressure, we try to reproduce the observation of Wang and Xie [1998] that the permeability within the fault remains high even after the pore pressure has fallen below the opening pressure. They suggest that this could be due to microscopic roughness on the fracture surfaces that would keep some space open along the closed fracture, called `fault smear'. Our closing pressure can be understood as a minimum pore pressure required to support this remaining open space in the fractures. It is expressed as a fraction of the lithostatic pressure, arbitrarily set between 75 and 80%.

4.6.3 Results

    In addition to the permeability in the fault zone and to the value of the closing pressure, another parameter essential to the evolution of the temperature field is the horizontal permeability in the sands, which controls the ease for the ascending fluids to spill into the reservoirs. Otherwise, fluids continue to the seafloor and appear as seeps, which have been observed in abundance in the area [Roberts, 1998]. This variable is defined by the permeability anisotropy in the sands (the ratio of the horizontal to vertical permeability).
 

Figure 4.11
Figure 4.12
Figure 4.13
Figure 4.14

    The results of different simulations for various values of these parameters are shown in Figures 4.11, 4.12, 4.13 and 4.14. In Figure 4.11 the closing pressure is 80% of the lithostatic pressure, the permeability enhancement in the fault is 102, and the permeability anisotropy in the sands is 102 . In Figure 4.12, the parameters are the same except the permeability enhancement in the fault, which is set to 103 in the shales, 102 in the proximal deltaic sands and 101 in the fluvial sands. In Figure 4.13 parameters are the same as in Fig. 4.12, except for a 103 permeability anisotropy in the sands. Finally, the model in Fig. 4.14 is similar to Fig. 13 except for the closing pressure ratio at 75%. In all the figures, (a) shows the evolution through time of the entire temperature field and (b) shows the evolution of the anomaly created by the fluid migration after removal of the salt influence, and therefore should be compared to Figure 4.9.

    In all cases, the migration was mostly confined within the fault at the beginning of the simulation for at least 1,000 years, generating the distinct anomalies in the first steps (left) in figures 4.12b,4.13b and 4.14b. Maximum calculated velocities during this initial phase are up to 1.0 10-6 m/s, or 30m/year. In these first stages, the temperature anomalies are localized only in the fault nodes. As upgoing fluids reach the interbedded sand/shales, they start flowing into the sands and the reservoirs, broadening the temperature anomaly significantly in the cases with high horizontal permeability (Figs 4.13b and 4.14b). After the fault closes, migration is mostly sub-horizontal, filling the reservoirs, and diffusing rapidly the temperature anomaly which becomes a low-amplitude, broader feature after 10,000 years. The exception is in Figure 4.12b where the anomaly remains focused because the sand permeability is lower. Another observation common to the different simulations is that the amplitude of the anomalies is generally higher at shallower depths (1750, 2000 mbsf), because the contrast in temperature between the ascending fluids and the formation temperature is higher in the shallower, colder, formation.

4.6.4 Discussion

    The significant differences between Figures 4.11b, 4.12b, 4.13b and 4.14b show that the three parameters varied between the models all have an influence on the duration, the amplitude and the distribution of the temperature anomalies.

    Figure 4.11 shows the only case where the permeability enhancement at the opening of the fault is 102 in all the formations. In this case, most of the fault nodes remained open to flow during the 10,000 years of the simulation. The anomaly develops progressively to its maximum (6° C) at the end of the simulation, which is lower than after only 1,000 years of any of the other simulations and smaller than the anomalies observed (Figure 4.9). This shows that in this case, the migration is too slow to carry enough heat in the shallower formations to reproduce the anomalies observed in EI330.

    By comparison, Figure 4.12 shows that anomalies of the same amplitude as what is observed in EI330 (i.e. ~10° C) can develop in about 1,000 years of activity in the fault if the permeability in the shales is increased by three orders of magnitude (103) when the fault opens. However, due to the low horizontal permeability in the sands, the ascending fluids are not deviated in the reservoirs and keep migrating upwards to the shallower fluvial sands. After the fault has closed, which is complete after about 5,000 years, the absence of lateral migration maintains the anomalies strong and localized in the fault nodes, because the only dissipation of this heat occurs through diffusion. The amplitude of the anomalies in this figure suggests that the 103 permeability enhancement in the faulted shales could be representative of the processes in EI330, showing that the fault opening has to increase the permeability of the shales by at least three orders of magnitude to allow the overpressured fluids to be released. However the anomalies observed in Figure 4.9 are broader and less focused than in any stage of Figure 4.12b, indicating that some of the heat gets carried into the formations adjacent to the fault.

    In Figures 4.13 and 4.14, where fault permeability is the same as Fig. 4.12, but horizontal permeability in the sands is higher, the anomaly is also mostly limited to the fault nodes until about 2,000 years of simulation. When the fault closes after this time, the anomalies become wider than in Figure 4.12, and more similar to Figure 4.9. In these two cases, the fault closes earlier than in Figure 4.12 because the fluids flowing into the adjacent formation release some of the pressure in the fault nodes. The only difference between Figures 4.13 and 4.14 is in the deeper part of the model (2500 and 2750 m) where the anomaly is slightly higher in Fig. 4.13b (where the closing pressure ratio is 80%) than in Figure 4.14b (where it is 75%). The small difference between the opening and closing pressures in the first case allows the pressure to build up rapidly after the closing of the fault and to reach the opening pressure in a short time, maintaining some semi-continuous activity in the deeper part of the fault for almost 5,000 years before closing definitely. In the case of Figure 4.14, where the closing pressure is 75% of the overburden, the pressure is too far below the opening pressure when it closes after 1,000 years to rebuild to this value.

Figure 4.15 Figure 4.16

    For a finer comparison with the observed data, Figures 4.15 and 4.16 show the simulated temperature distribution and advective anomaly in the case corresponding to the model in Fig. 4.13 after 5,000 years (These figures have the same format as Figures 4.5 and 4.9). Among the different simulations, this case reproduces most closely the field observations. Comparison between Figures 4.15 and 4.5 shows that the main features of the present-day temperature distribution are reproduced by our simulation, including the influence of the salt, the higher temperatures in Block 338 (South, center) in the deepest layers (2,500 and 2,750 mbsf), and in the West part of the field (Blocks 331, and 314) at shallower depths. These two separate domains of activity are emphasized in Figure 4.16 (to compare with Fig. 4.9) where the highest anomaly at depth is in the south of the field, and the anomaly becomes stronger in the western part at 1750 and 2000 mbsf. These two domains correspond respectively to the "Blue" Fault, which seems to be more active at depth, and to the "Red" fault (see Figure 4.3a) which is the dominant migration path to the shallower reservoirs.

4.6.5 Limitations of the results

    Even if the temperatures simulated in Figure 4.15 and 4.16 were reproducing more closely the data in Figure 4.5 and 4.9, it would be unreasonable to assume that these models indicate the exact origin of the observed anomalies. The necessary reservations about the accuracy of these results come mostly from (1) the limited quality of the temperature data and (2) the simplifications of the numerical modeling, and particularly (3) the coarseness of the model grid, especially around the fault.

    Despite the care taken in processing the BHT's, the possible error in the data is on the order of the anomaly amplitudes that we try to explain (4% or about 5° C, Lachenbruch and Brewer, [1956]). The only assurance of the reality of the interpreted temperatures comes from the clarity of some unquestionable features such as the influence of the salt diapirs. Figure 4.9 shows that the anomalies that we attribute to heat advection are indeed located around the faults, but the pattern of these anomalies is far from being as clearly related to the faults as the simulated anomalies in Figures 4.11b, 4.12b, 4.13b and 4.14b or 4.16. Part of this discrepancy must come from the quality of the temperature data. It is unfortunately impossible to do any more with the present data set because it was never collected for such purposes. If our results are valid, the simplicity of our methodology should encourage a systematic record and registration of such data sets in the future, including multiple measurements at the same depth when possible. While it doesn't require any expensive or time consuming procedure, it can provide a comprehensive representation of the hydrologic activity in similar, extensively drilled areas.

    Among the most flagrant simplifications of our model, the uniform permeability attributed to elements not smaller than 500´ 500´ 100m can not account for the heterogeneity in shale distribution and permeability on the scale observed in individual reservoirs [Chapter 3; He, 1996], or for the intricate distribution of individual faults and fractures in the fault zone [Losh et al., 1999, Losh, 1998]. The heterogeneity in shaliness within a reservoir is an important factor in reservoir simulations over short periods (10 to 20 years of production), and it controls the path of the forced migration and the distribution of the various fluid phases [Chapter 3]. However, the influence of this heterogeneity on the temperature distribution is more limited. The diffusion into adjacent formations of the heat carried by centuries-long advection prevents small scale heterogeneity in the temperature distribution. In any case none that could be detected by the low spatial resolution of our temperature field reconstruction. An effective permeability value that takes into account the average sand and shale content and the horizontal anisotropy due to the sediments layering can actually be representative of such large scale permeability [McCarthy, 1991]. More concern regarding the numerical simplifications can come from the choice of a rectangular grid. Narasimhan and Witherspoon [1976] report that for maximum accuracy of the Integrated Finite Difference Method (IFDM) used in the model, interfaces between elements should be perpendicular to the line joining the two center points and intersect that line at a mean position. For all practical purposes, this limits us to rectangular elements, unless we use advanced mesh generators. It forces flow to occur only across vertical and horizontal boundaries and follow step-like paths instead of running along the actual dip of the bedding or of the fault. Practically, this increases the volume of rock interacting with the fluids [Fisher et al., 1994], and dissipates artificially some of the heat carried by the fluids. It eventually leads to a reduction of the amplitude of the thermal anomalies simulated at shallower depth.

    Finally, regarding the coarseness of the grid, Roberts et al. [1996] give a fine scale description of the fault activity by using much narrower elements to define the fault (10 m wide, with fault apertures of the orders of microns), and successfully reproduce in 2D some of the temperature anomalies in Eugene Island. However, similar refinements in the vicinity of the fault planes for a rectangular grid in 3D would multiply the number of nodes and expand to unpractical dimensions the size of the linear equation systems to solve. Moreover, Losh [1998] shows that the fault zone crossed by the Pathfinder well is more than 100 m thick, and the various interpretations of adjacent fault splays [Rowan et al. 1998, Alexander and Handschy, 1998, Holland, 1991] suggest that migration could occur along wide paths, flowing along several fault segments which are only partially constrained [Rowan et al, 1998], or flowing across faults [Alexander and Handshy, 1998]. Therefore, while the results of Roberts et al. [1996] provide the explanation for the expulsion mechanisms along a single fault, the uncertainties in picking individual faults in 3D seismic and the large number of faults in the data-cube could barely justify such high resolution to realistically represent the Eugene Island fault zone in 3D. We assume that the permeability of our fault nodes represents an average bulk permeability, similar to the 1D model of Wang and Xie [1998], which distributes over the entire node the same amount of fluid that is actually confined within the individual fault planes crossing the node.

4.7 Conclusion

    The purpose of these models is to offer a method, relatively simple in its assumptions and realization, for analyzing in 3D the processes occurring in an active basin. Despite all the care that can be taken in processing the BHTs, the way they are collected precludes the possibility of establishing reliably an exact fine-scale temperature distribution. Therefore, even the most legitimate interpretation can be only qualitative. However, the high density of the data in EI330 and in similar densely sampled oil fields can somewhat compensate for the lack in quantitative accuracy of individual measurements. 3D models built on a well-constrained lithology defined by 3D seismic interpretation can provide a complete, semi-quantitative description of the fluid migration activity within the basin. Our results show that the temperatures observed today in the Eugene Island 330 field could be the result of 1,000- to 5,000-years long pulses of fluid expulsion activity within the fault. Our most important conclusions from this modeling is that fluids expelled from the overpressured shales migrate predominantly upwards within the fault during the period of activity, and then flow into the adjacent formation after the fault closes, carrying their heat, and hydrocarbons, into the most productive reservoirs of the field.


Table 1 : Principal physical properties of the sediments in the models


Property
Salt
Shale
Proximal Deltaic Sands
Fluvial  Sands
Thermal Conductivity
6
1.9
1.9 - 2.2
2.0
Vertical Permeability (m2)
10-19
10-17
10-15
10-14
Horizontal Permeability (m2)
10-19
10-15 - 10-14
10-13 -10-12
10-13
Permeability increase
no fracture
103
102-103
101 - 102
% of overburden for opening
no fracture
85%
85%
85%
% of overburden for closing
no fracture
75-85%
75-85%
75-85%
Porosity (%)
10
10
30
40

Table 2: Numerical parameters

Parameter
Value
Seafloor Temperature
18° C
Uniform Vertical Heat Flow at the base
60 mW/m2
Seafloor Pressure
106 Pa
Bottom Pressure
1.4 108 Pa


References

Alexander, L.L., and Flemings, P. B., 1995, Geologic evolution of a Plio-Pleistocene salt withdrawal Mini-Basin: Block 330, Eugene Island, South Addition, Offshore Louisiana, AAPG Bull., 79(12), 1737-1756.

Alexander, L.L., and J.W. Handschy, 1998, Fluid flow in a faulted reservoir system: Fault trap analysis for the Block 330 field in Eugene Island, South Addition, Offshore Louisiana, AAPG Bull., 82(3), 387-411.

Anderson, R.N., He, W., Hobart, M.A., Wilkinson, C.R., and Nelson, H.R. Jr., 1991, Active fluid flow in the Eugene Island Area, offshore Louisiana, Geophysics: The Leading Edge, 10(4),12-17.

Anderson, R.N., 1993, Recovering dynamic Gulf of Mexico reserves and the U.S. energy future, Oil & Gas Journal , April 26, 85-91.

Anderson, R.N., Flemings, P.B., Losh, S., Austin, J., Woodhams, R., and the Global Basins Research Network Team, 1994, Gulf of Mexico Growth fault drilled, seen as oil, gas migration pathway, Oil and Gas Journal, June 6, 97-103.

Armstrong, P.A., D.S. Chapman, R.H. Funnell, R.G. Allis and P.J.J. Kamp, 1996, Thermal Modeling and Hydrocarbon Generation in an Active-Margin Basin: Taranaki Basin, New Zealand, AAPG Bull., 80(8), 1216-1241.

Bachu, S., J.C. Ramon, M.E. Villegas, and J.R. Underschultz, 1995, Geothermal Regime and Thermal History of the Llanos Basin, Colombia, AAPG Bull., 79(1), 116-129.

Batzle, M., and Z. Wang, Seismic properties of pore fluids, Geophysics, 57(11):1396-1408, 1992.

Bodner D.P. and J.M. Sharp, 1988, Temperature Variations in South Texas Subsurface, AAPG Bull., 72(1), 21-32.

Bodvarsson, G. S., Mathematical modeling of the behavior of geothermal systems under exploitation, Ph.D. thesis, 353 pp., Univ. of Calif., Berkeley, 1982.

Bullard, E. C , 1947, The time necessary for a borehole to attain temperature equilibrium , in Mon. Not. R. Astron. Soc. Geophys., Suppl. 5, 127-130.

Burnham, A.K., 1995, Chemical kinetics and shale process design, in C. Snape (ed.), Composition, Geochemistry and Conversion of Oil Shales, Kluwer Academic Publishers, 263-276.

Castoviejo R., R. Gable, R. Cueto, J.C. Foucher, M. Soler, J. Gounot , J.C. Batsale, A. Lopez and M. Joubert , 1996, Ensayo de una metodologia innovadora para la deteccion de masas polimetalicas profundas: Modelo geologico y exploracion geotermica preliminares de la Masa Valverde (Huelva), Boletin Geologico y Minero, 107(5-6), 485-509.

Coelho, D., A. Erendi, and L.M. Cathles, 1996, Temperature, pressure and fluid flow modelling in Block 330, South Eugene Island using 2D and 3D finite element algorithms, in Annual Meeting Abstracts - AAPG and SEPM, v. 5, p. 28.

England, W.A., A.S. Mackenzie, D.M. Mann and T.M. Quigley, 1987, The movement and entrapment of petroleum fluids in the subsurface, J. of the Geol. Soc., London, 144, 327-347.

Fisher, A.T., K. Becker and T.N. Narasimhan, Off-axis hydrothermal circulation: parametric test of a refined model of processes at Deep Sea Drilling Project/Ocean Drilling Program site 504, J. Geophys. Res., 99(B2), 3097-321,1994.

Flemings, P., Zoback, M.D., Bishop, B.A., and Anderson, R.N., 1995, State of stress in the pathfinder well, in Anderson, R.N., Billeaud, L. B., Flemings, P., Whelan, J., and Losh, S., The GBRN/DOE Eugene Island Dynamic Enhanced Recovery project, CD-ROM, Lamont Press, Part IV-3.

Gable, R., 1977, Temperature and heat flow in France, Seminar on Geothermal Energy, Brussels, 6-8 December 1977, in Commission of the European Communities, Directorate General for Research, Science and Education, EURO 5920, Vol. 1: 113-131.

Gable, R., 1986, Acquisition et Interpretation de données Géothermales, Doctor es-Sciences Thesis, Université Pierre et Marie Curie, Documents du BRGM, 104.

Hart, B.S., P.B. Flemings, and A. Deshpande, Porosity and pressure: Role of compaction disequilibrium in the development of geopressures in a Gulf Coast Pleistocene basin, Geology, 23(1), pp 45-48.

Holland, D. S., Leedy, J. B., and Lammlein, D. R., 1990, Eugene Island Block 330 field-U.S.A., offshore Louisian, in Beaumont, E.A., and Foster, N. H., compilers, Structural traps III, tectonic fold and fault traps: AAPG Treatise of Petroleum Geology Atlas of Oil and Gas Fields: 103-143.

Hooper, E.C.D., 1991, Fluid migration along growth faults in compacting sediments, J. of Pet. Geology, 14(2): 161-180.

Horner, D.R., 1951, Pressure build-up in wells, Proc. Third World Petr. Congress, The Hague, section II, 503-521.

Humphris, C. C., Jr., 1978, Salt movements on continental slopes, northern Gulf of Mexico, in A. H. Bouma et al., eds., Framework, facies, and oil trapping characteristics of the upper continental margin: American Association of Petroleum Geologists Studies in Geology 7: 69-85.

Lachenbruch, A.H., and Brewer, M. C., 1959, Dissipation of the temperature effect of drilling a well in Arctic Alaska, U.S.G.S. Bull., 1083-C: 73-106.

Losh, S., 1998, Oil Migration in a Major growth fault: Structural Analysis of the Pathfinder Core, South Eugene Island Block 330, Offshore Louisiana, A.A.P.G. Bull., 82(9), 1694-1710.

Losh, S., L. Eglington, M. Schoell, and J. Wood, 1999, Vertical and Lateral fluid flow related to a large growth fault, South Eugene Island Block 330 field, Offshore Louisiana, A.A.P.G. Bull., 83(2), 244-276.

Mallard, W.G. and P.J. Linstrom, Eds., NIST Chemistry WebBook, NIST Standard Reference Database Number 69, March 1998, National Institute of Standards and Technology, Gaithersburg MD, 20899 (http://webbook.nist.gov).

McCarthy, J.F., 1991, Analytical models of the effective permeability of sand-shale reservoirs, Geophys. J. Int., 105, 513-527.

Mello, U.T., R.N. Anderson and G.D. Carner, 1994, Salt restrains maturation in subsalt plays, Oil and Gas Journal, Jan. 31:101-107.

Narasimhan, T.N. and P.A. Witherspoon, 1976, An integrated finite difference method for analyzing fluid flow in porous media, Water Resour. Res., 12(1), 57-64.

Press, W.H., S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, 1992, Numerical recipes in C: the art of scientific computing - 2nd edition, Cambridge University Press, New York.

Roberts, H.H., 1998, Evidence of episodic delivery of fluids and gases to the seafloor and impacts of delivery rate on surficial geology (Gulf of Mexico), AAPG 1998 annual meeting, expanded abstracts.

Roberts, S.J., J.A. Nunn, L. Cathles, and F.D. Cipriani,1996, Expulsion of abnormally pressured fluids along faults, J. Geophys. res., 101, B12, 28,231-28,252.

Rowan, M.G., B.S. Hart, S. Nelson, P.B. Flemings and B.D. Trudgill, Three-dimensional geometry and evolution of a salt-related growth-fault array: Eugene Island 330 Field, offshore Louisiana, Gulf of mexico, Mar. and. Pet. Geol., 15, 309-328.

Rybach, L., 1981, Geothermal Systems, Conductive Heat Flow, in Geothermal Systems: principles and case histories, edited by L. Rybach and L.J.P. Muffler, John Wiley and Sons, Ltd.

Sharp, J.M., and Domenico, P.A., 1976, Energy transport in thick sequences of compacting sediments, Geol. Soc. of Am. Bull., 87: 390-400.

Wang, C.Y., and X. Xie, 1998, Hydrofracturing and Episodic Fluid Flow in Shale-Rich Basins - A numerical Study, A.A.P.G. Bull., 82(10), 1857-869.

Whelan, J., Eglington, L., Requejo, R., and Kennicutt, M. C., 1995, Pathfinder well organic geochemistry - Indicators of Oil source and Maturity and fluid flow mechanisms, in Anderson, R.N., Billeaud, L. B., Flemings, P., Whelan, J., and Losh, S., The GBRN/DOE Eugene Island Dynamic Enhanced Recovery project, CD-ROM, Lamont Press, Part V.