TECTONOPHYSICS
ELSEVIER
Tectonophysics 233 (1994) 177-192
The thermal conductivity and heat flow density of the Jeanne
d’Arc Basin, offshore eastern Canada
Zehui Huang *, Mark A. Williamson zyxwvutsrqponmlkjihgfedcbaZYXW
Atlantic Geoscience Centre, Bedford Institute of Oceanography, P.O. Box 1006, Dartmouth, N.S. B2Y 4A2, Canada
(Received April 8, 1993; revised version accepted December 13, 1993) zyxwvutsrqponmlkjihgfedcbaZYX
Abstract
In an effort to understand the present-day thermal structure of the Jeanne d’Arc Basin, we derive matrix thermal
conductivity values, using a geometric mean model, from thermal conductivity values of thirteen lithologic end-members (of which seven are measurements on materials from the Jeanne d’Arc Basin) and lithological data. Temperature corrections to thermal conductivity are made for both matrix and pore water. In-situ thermal conductivitics
were obtained after correcting for the effect of porosity, which we estimated from both lithological information and
digital well logs. Calculated thermal conductivities of the formations range from 1.46 to 2.65 W mP’K-‘. The
heterogeneity of formation thermal conductivity is caused by at least two factors. The first is lithological variation
due to facies change. The second is due to different degrees of compaction. Under-compacted intervals (overpressure zones) are characteristically of lower thermal conductivity. The surface heat flow density calculated with an
inverse method in seventeen wells varies from 37.8 to 60.7 mW/m2. This study indicates that surface heat flow
density values estimated without porosity corrections may be 15% higher than those with porosity corrections. We
also demonstrate the sensitivity of modelled timing of maturation to thermal conductivity estimations.
1. Intr~uction
A sedimentary basin’s thermal structure, which
includes thermal conductivity of the rock column
and heat flow density, particularly its thermal
evolution through time, is important for understanding its petroleum generation, migration and
accumulation histories.
The Jeanne d’Arc Basin, offshore Newfoundland (Atlantic Canada), is an important hydrocarbon-bearing area and has undergone extensive
* This is Geological Survey of Canada Contribution No. 12092.
exploration (Bell and Campbeh, 1990). Previous
work exploring the basin’s thermal structure can
be found in Correia et al. (1990) who reported
laboratory measurements on the thermal conductivity of seven major lithologies and thermal gradients from corrected bottom-hole temperatures.
Using measured thermal conductivity values as
well as data from literature, and corrected thermal gradients, Correia et al. (1990) estimated
terrestrial heat flow density at 35 well locations in
the basin and proposed a schematic basinal
fluid-flow pattern, which describes a downward
component in the north and an upward component in the south. Their study, however, did not
HO-1951/94/$07.00
0 1994 Etsevier Science B.V. All rights reserved
0040-1951(93)EO278-3
SSDI
178
Z. Huang, MA. Williamson / Tectcmophysics 233 (19941177-192
account for the effect of porosity change with
depth upon thermal conductivity, which is an
important factor controlling thermal conductivity
of sedimental
rock.
To provide a more detailed examination of the
thermal properties of the Jeanne d’Arc Basin, we
established in-situ thermal conductivity profiles at
seventeen well locations (Fig. 1) using lithological
information (percentage of thirteen lithologies
from cutting sample examination) and published
thermal conductivity values. Porosity corrections
were performed using digitized sonic and density
logs. On the basis of porosity-corrected
thermal
conductivity, we used a thermal resistance (Bullard, 1939) and an inversion method (Tarantola
and Valette, 1982; Vasseur et al., 1985; Wang,
48*
490
0
I
0
I
km
I
40
,
n t-23
470
/C-92
n
l c\
Hebron
l-13
46*
Fig. 1. The Jeanne d’Arc Basin and well locations (0 = wells used in this study) (modified from Grant et al., 1986).
Z. Huang, M.A. Williamson / Tectonophysics 233 (1994) 177-l 92
1992) to calculate the heat flow density at well
locations. Derived heat flow densities and thermal conductivity profiles are of use in source rock
maturi~ modelling and can be used to test model
sensitivity to assumed thermal properties. The
in-situ thermal conductivity values were used in a
study (Williamson, 1992) that modelled maturity
of source rocks within the Jeanne d’Arc Basin.
In our study, the stratigraphic framework follows that described in McAlpine (1990). Both
lithological data (results of cutting sample examination), logging data and temperature measurements were retrieved from the Atlantic Geoscience Centre’s database.
2. Geological setting
The Jeanne d’Arc Basin is about 16,000 km2 in
area and is a sub-basin of the East NewfoundIand
Basin. It stretches from southwest to northeast
and is flanked by the Bonavista Platform to the
west and the Outer Ridge Complex to the east
(Fig. 1).
Formed as a result of the Mesozoic extensional rift tectonics that produced the North Atlantic Ocean (Mc~pine,
1990), the oldest sediments of the basin are Triassic red beds and salt.
Jurassic sediments are dominantly shallow-water
carbonates, shales and sandstones. The Lower
Cretaceous strata record several major progradational deltaic pulses. Upper Cretaceous sediments are mainly calcareous shales deposited
during marine transgressions. Tertiary strata are
composed of marine shales with minor sandstones. The lithostratigraphy (Fig. 2) (McAlpine,
1990) shows a thick Mesozoic-Cenozoic sequence
with major hiatuses at Kimmeridgian, late Barremian, late Aptian and late Albian-Cenomanian
times.
The Jeanne d’Arc Basin can be divided into
three structural-stratigraphic
segments by transbasin fauns (Fig. 1). South of the Egret Fault
(Enachescu, 1987), uppermost Jurassic and Lower
Cretaceous rocks are virtually absent because of
erosion and depositional thinning. Between the
Egret Fault and the Hibernia-Ben
Nevis transbasin fault zone this section thickens dramati-
179
cally, North of the trans-basin fault zone, a thickened Middle Cretaceous sequence has been preserved. The Hibernia oil field is located in the
Hibernia-Ben
Nevis trans-basin fault zone. The
Terra Nova oil field lies south of the trans-basin
fault zone.
3. Methods
For each well location, the matrix conductivity
of an interval of rock (at least 3 m thick) was
calculated using the geometric mean model (Eq.
(1); Woodside and Messmer, 1961a,b; Brigaud et
al., 1990):
h, = II@
(1)
h, is the matrix conductivity of aggregates of
several constituents, Ai is the conductivity of the
ith constituent and Pi its fractional volume. The
matrix thermal conductivity of thirteen individual
lithologies (end-members) are from Correia et al.
(1990) and Beach (1985) (Table 1). The values for
limestone, dolomite, anhydrite, shale, sandstone,
siltstone and claystone (the most commonly encountered lithologies in the study area) were
measured on samples from the Jeanne d’Arc
Basin and have been corrected for anisotropy
(Correia et al., 1990). These values therefore represent thermal conductivity perpendicular to the
bedding. The fractional volume of a lithological
component is appro~mated by its percentage estimation through cutting sample examination.
The thermal gradients (G,) for each well location (Fig. 1) are from Correia et al. (1990) which
were calculated from corrected bottom hole temperature (BHT) values (Table 2). These were
used to estimate temperature
(7’) at specific
depths (d) in Eq. (2), for the purpose of correcting temperature
effect on the matrix thermal
conductivity. The temperature at the sea floor
(T,) is taken as 4”C, following Hyndman et al.
(1979). The effect of temperature on matrix thermal conductivity has been accounted for following Chapman et al. (1984) in Eq. (3), in which
bGiE
Ml) I
EPOCH
AGE
FORMAflONS
I
zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJ
IC
2c
/
SC
SANDSTONE
4c
SANOUEREAU
FORMATION
DOMNANTLY SHALE
AND MWSTONE
t
-.
-.
-_
SILtSTONE
SC
so
70
so
DOLCMKE
so
00
SALT
10
VOLCANICS
I I zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCB
h” ‘AYALON
f
20
30
zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
IO
IO
LiMEStONE
BEDS
30
IO
lo
MNGCOMERATE
lo
0
0
0
a
0
_.-.._._.__
-.-._.-.-.-.-.-._
-.-.
-. _ “-7__
9
,
E
[1
e
PRE-MESOZOIC
BASEMENT
D
D
c)
D.W.0
*
s
D
0
D
Z. Huang, M A.
Table 1
Thermal conductivity, sonic travel time and density of thirteen
matrix components
Limestone
Dolomite
Anhydrite
Shale
Sandstone
Siltstone
Claystone
Ma&tone
Chert
Coal
Conglomerate
Salt
Unclassified
181
W illiamson / Tecronophy sics 233 (1994) 177- 192
Table 2
Corrected thermal gradient (“C/km) in wells analyzed in this
study (from Correia et al., 1990)
Conductivi~
(W m- “K-‘1
Sonic travel time
@,‘ft)
Density
(g/cm’)
Well name
Corrected thermal gradient
PC/km)
2.41+ 0.96 (1)
3.04k1.22 (1)
5.23*2.09 (1)
2.30 f 0.92 (1)
4.19+ 1.68 (1)
4.24+ 1.70 (1)
3.11+ 1.24 (1)
3.00+ 1.1 (2)
1.40+0.5 (2)
0.20*0.2 (2)
3.28 + 1.8 (2)
5.70+ 1.0 (2)
2.30 i 2.0 (2)
47.5 (3)
43.5 (31
50.0 (3)
70.0 (7)
53.3 (3)
55.0 (5)
70.0 (6)
48.0 (4)
51.0 (8)
120.0 (11)
54.6 (9)
67.0 (3)
50.7 (10)
2.71 (3)
2.87 (3)
2.96 (3)
2.65 (3)
2.65 (3)
2.65 (5)
2.65 (6)
2.70 (4)
2.65 (11)
1.50 (3)
2.69 (9)
2.17 (3)
2.67 (10)
Ben Nevis 1-45
Beothuk M-05
Flying Foam I-13
Gambo N-70
Hebron I-13
Hibernia C-96
Hibernia K-14
Hibernia K-18
Hibernia O-35
Mercury K-76
Nautilus C-92
Rankin M-36
South Tempest G-88
Terra Nova K-08
West Flying Foam L-23
Whiterose J-49
Whiterose N-22
25.4 rt 4.8
24.9f 3.7
26.7 + 0.0
25.7 it 2.9
24.Ok3.4
25.1 zt 1.3
22.8t 1.1
26.7rt2.1
23.7 i: 4.7
22.6 k 0.3
28.9 + 3.3
21.3+ 1.8
31.6k6.8
25.7+ 1.2 ’
29.8k7.9
31.97t 1.s
31.8+2.0
Notes: (1) from Correia et al. (1990); (2) from Beach (1985);
(3) from Schlumberger (1989); (4) supposed to be near limestone value; (5) from Hagarty et al. (1988); (6) supposed to be
near shale value; (7) from Brigaud et al. (1990); (8) supposed
to be near condensed sandstone value; (9) average value of
type 1 to type 9 (except type 2); (10) average value of type 1 to
type 12 (except type 10); (11) from Doveton (1986).
a From a nearby well Terra Nova K-18.
lithologic end-members. The lithology weighted
Eqs. (5) and (6) (Schlumberger, 1989), respecfor temperature: zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
tively, are used, in which Pi is the fractional
T= T,+G,d
(2)
volume of the ith constituent approximated by its
293
percentage given by cutting sample examination.
Am,T = A(3)
Matrix sonic travel time (to and matrix density
“273+T
(pi) of the thirteen lithologies (end-members) are
Pore-water thermal conductivity (A,,r), which
mainly from Schlumberger (19891, Hagarty et al.
also changes with temperature is estimated from
(19881, Brigaud et al. (1990) and Doveton (1986)
Eqs. (4a) and (4b) (Kappelmeyer and Haenel,
(Table 1):
1974):
A m,T is the matrix thermal conductivity corrected
A w,T= 0.56 + 0.003~“.827
A w,T= 0.442 + 0.0519 In T
0 =GT G 50°C
T > 50°C
(da)
= k AtiPi
(5)
i=l
(4b)
An important step in calculating the in-situ
thermal conductivity of a rock interval is to account for the effect of porosity. We used sonic
and density logs to estimate the porosity-depth
profile. The matrix sonic travel time (At,) and
matrix density (p,) of an interval of rock are
calculated from the lithological data, matrix sonic
travel time (ti>, and matrix density (pi) of the
Fig. 2. Lithostratigraphic
At,
Pm=
kPipi
i=l
(6)
Porosity values are estimated from averaged
readings on the sonic and density logs for each
interval, with Eqs. (7) and (8) (Schlumberger,
1989). In this study we did not use the time-average relationship proposed by Wyllie et al. (1956,
1958) to obtain porosity from sonic log, since, in
scheme used in this study (from McAIpine, 1990).
Z. zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
Huang, M A. W illiamson / Tectonophy sics 233 (1994) 177- 192
182
addition to its theoretical oversimplication, in
practice it is overly conservative in the 5-25%
porosity range and is optimistic for porosities in
excess of 30% (Raymer et al., 1980; RaigaClemenceau et al., 1988). Instead, we employed
an empirical equation based on the comparison
of sonic travel time versus independent porosity
measurements (Eq. 7) (Schlumberger, 1989). In
Eq. (71, 4, stands for porosity from the sonic log,
At,O, represents sonic log response and 0.67 is an
empirical constant. Eq. (8) is a simple mass balance relationship in which 4, stands for porosity
from the density log, plog for density log response, pr for the density of pore water taken as
1.0245 g/cm3 (Brereton, 1992):
+b,, = 0.67(&,,
- A&)/&,,
4hp= t&n -Pi&(&n
-Pf)
(7)
(8)
Two thermal conductivity vahres (A,,, A,) can
be ~al~uIated with the porosity values (&, 4,)
using Eqs. (9) and (101, which followed the same
geometric mean model (Woodside and Messmer,
1961a,b; Brigaud et al., 1990). The average of
these two values is taken as the estimate of the
in-situ thermal conductivity (A) (Eq. 11):
A = A(;
S”
W,TSV
(9)
A = A(~~#~)A~~
(10)
w,T
A”= (A, + Ap),‘2
(11)
where Azj, Ai are the thicknesses and conductivities of lithological layers in a well, respectively.
The therma resistance method has been applied to the Jeanne d’Arc Basin by Correia et al+
(1990) due to its simplici~. With this approach, a
single thermal gradient obtained from a linear
fitting of the subsurface temperature measurements is used. The problem in this approach is
that, in reality, due to radiogenic heat production
and overpressured intervals, the temperature-depth relationship at a well location is not always
linear and generahy, the thermal gradient varies
more or less from interval to interval in each well.
The inversion method is much more sophisticated than the thermal resistance one. It uses a
stochastic inversion algorithm (Tarantola and Valette, 1982; Vasseur et al., 1985; Wang, 1992)
based on the thermal resistance method in one
dimension. This method starts from the relationship between heat flow density (q), temperature
(7’) and thermal conductivity (A):
4 = A(dT/dz)
(14)
For obtaining the subsurface temperature at
depth z CT,), one can integrate Eq. (14) from the
sea floor (z = 0) to depth z if the sea floor
temperature, surface heat flow density and the
thermal conductivities of the sediment layers from
the ground to depth z are known:
3.2. Calculation of heat j?ow density
T(z)
We used two methods to calculate heat flow
density at each well. The first one is the reIatively
simple thermal resistance method. As proposed
by Bullard (1939), where heat production and
transient fluid movement are negligible, heat flow
density (mW/m2) can be calculated with the following equation:
= 7’, +4/e’ dz’/A( z’)
(15)
The inverse problem occurs when we want to
process simultaneousIy all the data, which include
m measurements of T,, 12 estimations on the
thermal conductivities of the lithological layers in
a well, an estimation of sea floor temperature,
and an uncertain estimation of surface heat flow
4 = A,,G zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
T
density, in order to find the optimal values of
(12)
these physical quantities which are subject to
where G, is the average thermal gradient over
errors. Taking all the measurements and estithe driled interval and A, is the average thermates as the a priori information, one can denote
mal conductivi~ of the ~rresponding
interval,
the temperature measurements by the m-dimenwhich can be obtained with Eq. (13):
sion vector d (dT = (ti, t,, . . . t,>) (T indicates
Aave= ~ AZi/C( AZi/A,)
transpose) and the sea floor temperature, ther(13) zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJ
i=l
mal conductivities and surface heat flow density
2. fiuang, M A. W illiamson
/ ~ectff~ophy ~ics
233 (1994) 177- 192
(referred to as parameters thereafter) by the (n +
2)-dimension vector p {pT = (to, A,, A,,...
A,,
4)). The physical law relating d and p as the
forward solution can be written as:
(16)
d = g(p)
where g is a vectorial function defined by Eq.
(15).
Assuming d and p are random vectors and that
the a priori joint probability density function of
the two is Gaussian, the joint density function is
therefore defined by the mean value (do, po), the
covariance matrix of do CC,,) and that of p. (C,,>
(for details see Vasseur et al., 1985). The mean
value do and p. can be defined with the available
a priori information (the temperature measurement set and the parameter set, respectively).
The diagonal elements of the covariance matrices
C,, and C,, are set equal to the a priori errors
which express our confidence in the a priori values. Usually measured values have small variance
but subjectively defined (estimated) value(s),
which in our case is surface heat flow density,
should be assigned with larger variances, The
off-diagonal elements of the covariance matrices
C,, and C,, express the correlation between the
uncertainties in different components of p. or do
and are assumed zero.
Since Eq. (15) is non-linear, the inverse problem is also non-linear. Tarantola and Valette
(1982) present an iteration scheme to search for
the most probable estimate of p and d:
pk+ 1= PO+ C,,G:(Cdd + G&&~)T
x [do - gtp/J + G,(P, - PO>]
dk+l = g(pk)
+ Gk(~k+, - PHI
(171
(18)
where pk is the estimation of p at the kth iteration and G, is a gradient matrix consisting of the
partial derivatives of Eq. (16) with respect to p at
point pk.
The optimal a posteriori information is yielded
when the largest element in the differences between pktl and pk becomes less than a specified
small value.
In the Jeanne d’Arc Basin, the measured subsurface temperature are given a priori error of
f 10°C. Sea floor temperature To is assumed to
183
be 4°C (Evans and Coleman, 1975; Hyndman et
al., 1979), with an a priori error of + 2°C. These a
priori errors for the temperatures may accommodate the variances related to measurement inaccuracy and seasonal changes. The porosity-corrected thermal conductivities are associated with
40% of their values as the relative a priori error,
according to the error estimation by Correia et al.
(1990). A priori surface heat flow density is assumed to be 95 mW/m2, with a large a priori
error of &90 mW/m2, since this is the least
known quantity.
4. Porosity-corrected thermal conductivity
In-situ thermal conductivi~ profiles were derived for seventeen well locations (Table 2) where
both lithological data and digitized sonic and
density logs were available.
In the studied wells, each cutting sample represents a 3-20-m interval. The well log readings
are at a finer interval (1 m). Therefore, the average of the log responses in the interval represented by the cutting sample is used to estimate
the porosity values.
A part of the in-situ thermal conductivity profile established for Beothuk M-05 (from the Ben
Nevis Formation to Hibernia Formation) is illustrated in Fig. 3. The in-situ thermal conductivity
profiles contain errors from various sources. The
first is the published thermal conductivity values
for the litholo~~ end-members (Beach, 1985;
Correia et al., 1990) for calculating matrix thermal conductivities. The second is from the thermal gradients calculated from corrected BHT values (Correia et al., 1990). The third is from our
porosity estimates using sonic and density logs.
We are unable, however, to estimate the error by
comparing systematically the porosity values from
well logs to limited porosity measurements. Generally, porosity values calculated from sonic logs
and those from density logs differ from each
other more or less with different lithologies. Assuming that average porosity value is more representative of the true porosity and the two values
from sonic and density logs represent the maximum and minimum porosities, we observed that
the differences between the average values and
the maximum,Iminimum values are mostly less
than 15% This is in agreement with estimated
error range (l&20%) of in-situ porosity calculation with logs CDella Yedova and Van Herzen,
1987).
In Fig. 3, the range of in-situ thermal conductivity for each interval resulted from the use of
two porosity va[ucs for porosity corrections. The
extent of the range is rerated to the difference of
porosity values from sonic and density togs, which
varies from formation to formation (Fig. 3), The
difference between the two porosity vah~es is
greater in shnly intervals than in other intervals.
We observed that shale/claystone
intervals have
large range of in-situ thermal conductiviw estimates.
At Beothuk M-05, the matrix thermal conductivity for Nautilus Shale and Fortune Bay Shale
varies from 2.3 to 4.1 W m-‘K-l,
and the in-situ
crJnductivity
purosity
WI
20
40
Fig. 3. From feft CO right are shown the Mmlogi4
pmfile, density log, smGc lag, porosity values from kholugical
data and logs,
and thermal conductivity
at Beothuk M-05. The. lithdogical
profile shows the percentage
of carbonate
(A), shaIc,‘clavstone
(B),
sandstone
CC), siltstonc iD1 and coal <El in the Ben Nevis Formatian
CO, Avalon I+rmation
(21, Eastern Shoal Formation
131,
Whiterose Shaie (41, Catalina Formation
(5), “B” marker (6) and Hibernia T;ormation (7). The diagram to the far right shows Iroth
the matrix thermal conductivity corrected for temperature
(triangle) and the in-situ thermal conductivity corrected for temperature
and porosi@ (open square), also with the maximum and minimum values being indicated.
The depth iu kilometer below rotary
table.
Z. Huang, MA
Wiiiiamson / Tectonophysics
233 (1994) 117-l 92
185
matrix conductivities. The difference between the
in-situ thermal conductivity and matrix conductivity in the less compacted sandstone/siltstone/
shale interval of the Catalina Formation is about
1 W m-‘K-l, while that in the more compacted/
cemented limestone interval (lower part of the
Catalina Formation and “B” Marker) is less than
0.5 W m-‘K-i.
For this reason, although the
matrix thermal conductivity of the sandstone/
siltstone/ shale Catalina Formation is higher than
that of the carbonate interval, the in-situ thermal
conductivities of the sandstone/ siltstone/ shale
interval are not much higher than those of the
carbonate interval.
To apply the inversion method for estimating
surface heat flow density, we calculated formation in-situ thermal conductivity at each well. The
range and geometric mean of the formation thermal conductivity (A,) from the seventeen wells
are presented in Fig. 4. Shale/claystone
and zyxwvutsr
thermal conductivity from 1.5 to 2.8 W m- ‘K- ‘.
We also observed lower in-situ thermal conductivity (1.3 W m-“IX-‘) in other wells (e.g., Hibernia K-18). From Fig. 3, it can be seen that the
lithologic variations have a large control on the
variations of both matrix conductivi~ and in-situ
thermal conductivity. For example, the sandstone/ siltstone/ shale Catalina Formation has
higher in-situ thermal conductivity values than
the carbonate “B” Marker. The coal-bearing intervals in the Avalon Formation, Whiterose Shale
and Hibernia Formation usually show low matrix
and in-situ thermal conductivities. In the lower
Hibernia Formation, shale/ claystone alternate
with sandstone/ siltstone, with corresponding alternation of low and high in-situ thermal conductivities.
Porosity also exerts an important control on
the variation of in-situ thermal conductivity. Insitu thermal conducti~ties are always lower than
Thermal conductivity (W/m/K)
1.0
1.5
2.0
2.5
3.0
l....‘....I....‘....I
South Mara Unit
-
Da~~~FM
-
EiderFM
Nautilus Sh.
I
BenNevis
-
Avalon FM
FM
Eastern Shoals FM
Whiiercse SH. (U)
_
Catalina FM
‘B’ Marker
-
Whitercse SH’ (L)
Hiimia FM
Fortune Bay SH.
Jeanne d’Arc FM
-
RankinFM
o
#
sandstone/siltstine
shaleklaystone
I
sandstone
+
carbonate
0
sandstone/siltston&hale
Fig. 4. Geometric mean and the range of formation thermal conductivity for each formation with the dominant lithology specified,
based on the results from seventeen wells.
carbonate formations have greater range than
other types of formation. The geometric mean of
sandstone/ siltstone formations’ therma conductivity progressively increases with depth. However, the geometric mean of the thermal conductivity of shale/claystone
formations first increases with depth and then decreases with depth,
with the Whiterose Shale (Uf as the turning point.
The Fortune Bay Shale has the lowest geometric
mean of formation thermal conductivity. The
sandstone Hibernia Formation, the Jeanne d’Arc
Formation and the Rankin Formation also display a lower geometric mean of formation thermal conductivity than the overlying sandstone
Eider Formation and Ben Nevis Formation.
A
Beothuk Terra Nova
K-01
M-05
1
Hebron
1-13
This vertical heterogeneity of formation thermal conductivity may be due to overpressure of
formations below the “B” Marker. The top of
overpressure is near the top of or in the Fortune
Bay Shale (3-4 km below sea level). The fluid
retained in these overpressured formations explains their reduction in thermal conductivity.
Therefore, the “B” Marker, which has a low
geometric mean of formation thermal conductivity, and the overpressured formations can be considered as the barrier for effective thermal transfer in the Jeanne d’Arc Basin.
In Fig. 5 the spatial variation of formation
thermai conductivity aIong a cross section (the
A-B line shown in Fig. 11 is illustrated. Few
Hibernia
Hibernia
Nautilus
B
Flying Foam zyxwvutsrq
Fig. 5. Cross section showing formation thermal conductivity IW m-t Km’) in Beothuk M-OS, Terra Nova K-08, Hebron I-13,
Hibernia C-96, Hibernia K-18, Nautilus C-92 and Flying Foam I-13 from south to north along line A-B (see Fig. 1). The depth is
kilometer below sea level, and the horizontal space between the wells is in accordance with their distance. Except for FIying Foam
f-13, there are oil occurrences along this cross section (see Fig. 1). The formation boundaries are from McAipine (1990). The wavy
line between formations stands for unconformity and the straight line for conformable contact. The thick horizontal lines indicate
the depth of the top of overpressure zone in Hebron I-13, Hibernia C-96, Hibernia K-18 and Nautilus C-92
Z. Huang, M A.
UW iamson / Tectonophy sics 233 (1994) 177- 192
formations exhibit a N-S trend in formation thermal conductivity. The lateral variation can be
partly explained by facies changes and unconformaties. Both would change the proportion of
lithologies within a formation at different locations. As the shale/claystone
ratio increases, the
formation thermal conductivity decreases. Another reason is overpressure. In Hebron I-13 and
Hibernia K-18, the formations under the top
overpressure have rather low thermal conductivities. Although the Fortune Bay Shale and Jeanne
d’Arc Fo~ation
in Hebron I-13 are buried
deeper than their normally compacted counterparts in Terra Nova K-08 and Beothuk M-05,
they have lower formation thermal conductivities
because of under-compaction.
In Nautilus C-92,
the top overpressure is much higher (in the Nautilus Shale). The sediments in Nautilus C-92 are
deeply buried because of faulting. A comparison
of pressure data at the Nautilus C-92, Hibemia
K-18 and Hebron I-13 locations (Fig. 6) reveals
that Nautilus C-92 has the most overpressured
section, i.e., stratigraphically
and structurally
shallower than the other wells. As the result,
most of the overpressured formations in Nautilus
C-92, although buried deeper, have lower thermal
~nducti~ties
than their counte~arts in Hibemia
K-18, Hebron I-13 and other wells which are not
overpressured.
5. Swface heat flow density
In most of the wells, the uppe~ost
intervals
(mostly of Tertiary age) which are several hundreds to 2000 m thick do not have both detailed
lithological data and logging data. Therefore,
there is no in-situ thermal conductivity estimations on these uppermost intervals. To calculate
the surface heat flow density with the two methods mentioned, the in-situ thermal conductivities
of every lithological interval from the sea floor
down to the bottom of the well are needed. To
solve this problem, we assume that the thermal
conductivity of the uppermost layer is 10% lower
than that of the formation immediately below it,
because of less compaction (higher porosity). This
187
Pressure zyxwvutsrqponmlkjihgfedcbaZYXWVU
(KPa * 1000)
0
.Ir
1000
2000
B
n
3
g
3000
4000
5000
Fig. 6. Comparison of subsurface pressure in the overpressurized wells (Nautilus C-92, Hibernia K-18 and Hebron I-13.
The short horizontal bar indicates the top of Fortune Bay
Shale in Hebron I-13 (a), Hibernia K-18 (b) and Nautilus C-92
(cl.
assumption is compatible with our observation
made in Fig. 4.
The thermal resistance method is quite simple
(see Eq. 12). With the inverse method, the assumed surface temperature, subsurface temperature measurements (BHT in our cases), formation
thermal conductivity estimation and assumed heat
flow density with their assumed standard deviation are taken as the a priori info~ation
and the
calculation is to find the optimal values (i.e., the a
posteriori information, Tables 3 and 4). For example, the inversion reduced the a priori value of
surface heat flow density (95.0 rfr90.0 mW/m*) to
its optimal value of (50.0 + 7.4 mW/m*) (Table
Z. Huang, M .A. W kzamson / Tectonophy sics 233 (19941 177- f 92
188
Table 3
The a priori and a posteriori information of the temperature
(T, in “C) at Hebron I-13
Depth
(m)
a priori
value
a priori
SD.
a posteriori
value
Correction
(“0
779.7
2518.7
2569.7
2700.7
2845.7
3359.7
3723.7
3927.7
4052.7
4072.7
4170.7
4195.7
4323.7
4595.7
48.9
57.7
74.4
68.7
79.4
77.7
99.6
103.0
107.8
105.0
106.0
117.0
107.4
125.5
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
50.6
63.4
73.0
69.3
77.2
80.1
97.7
100.1
108.6
104.2
110.3
111.3
114.8
120.7
1.7
5.7
- 1.4
0.6
-2.2
2.4
- 1.9
- 2.9
0.8
-0.8
4.3
5.7
7.4
-4.8
4). Other parameters (for example, fo~at~on
thermal conductivity) arc also changed with reduction in standard deviation (Table 4).
Our estimates of the surface heat flow density
for seventeen wells in the Jeanne d’Arc Basin
derived from the two methods are summarized in
Table 5. For comparison, the estimations by Correia et al. (19901 are also listed. Our surface heat
flow density values from the thermal resistance
and inverse methods mentioned are quite close.
Usually, the inversion method gives slightly higher
value than the thermal resistance method. As the
inversion method is more rigorous than the thermal resistance method, we will mainly discuss the
results from inversion.
We observed that the surface heat flow density
calculated with the inverse method varies from
37.8 (Rankin M-36) to 60.7 mW/m2 (two Whiterose wells). The majority of the heat flow density
values from this basin is in a much narrower
range from 46 to 59 mW/m*. The average of
surface heat flow density from the seventeen wells
is 51.2 1fl3.1 mW/m2, and that from the fourteen
wells (with the maximum/minimum
values excluded) is 51.8 f 2.8 mW/m2. Therefore,
51
mW/m2 may be taken as the representative basina1 average for surface heat flow density.
In comparison with the results from Correia et
al. (19901, our values from either of the methods,
with few exceptions (e.g., Gambo N-701, are much
lower. According to Allen and Allen (1990), thermally subsiding (post-rift)/ passive margin usually
have heat flow ranges from 35 to 65 mW/m2,
with an average of 50 mW/m2. Our range (37%
60.7 mW/m2) is more compatible with this typi-
Table 4
The a priori and a posteriori information of sea floor temperature (T,, in “C), formation thermal conductivity (A,, in W m-‘K-l)
and surface heat flow density (9, in mW/m*) at Hebron I-13
Depth
cm)
T0
A, uppermost layer
A, South Mara Unit
Ar Dawson Canyon Fm.
A, Nautilus Sh.
A, Ben Nevis Fm.
A, Avalon Fm.
A, Eastern Shoal Fm.
A, Whiterose Sh. Kl)
A, Catalina Fm.
h,“B” Marker
Ar Hibemia Fm.
AF Fortune Bay Sh.
A, Jeanne d’Arc Fm.
A, Rankin Fm.
4
o-1221.7
1221.7-1325.7
1325.7-1582.7
1582.7-1737.1
1737.7-1803.1
1803.7-1867.7
1867.7-1969.7
1969.7-2525.7
2.525.7-2707.7
2707.7-2764.7
2764.7-3367.7
3367.7-3647.7
3647.7-4300.7
4300.7-?
a priori
value
a priori
SD.
a posteriori
value
a posteriori
S.D.
4.0
1.65
1.84
2.10
1.86
2.30
2.12
2.40
2.53
2.65
1.75
2.14
1.62
1.88
2.45
95
2.0
0.66
0.74
0.84
0.74
0.92
0.85
0.96
1.01
1.06
0.70
0.86
0.65
0.72
0.98
95
3.88
1.90
1.91
2.41
1.70
2.32
1.91
2.63
1.99
2.15
1.83
2.12
1.77
1.73
2.34
50
1.99
0.58
0.65
0.68
0.50
0.74
0.57
0.73
0.51
0.64
0.45
0.60
0.50
0.57
0.97
7.4
189
Z. Huang, MA. wiitiamson / Tect~no~~ysic~233 (1994) 177-l 92
cal range than that (48-75 mW/m*) from Correia et al. (1990). If Gambo N-70 is excluded, the
difference in surface heat flow density from the
inverse method and those estimated by Correia et
al. (1990) varies from 2.3 to 20 mW/m2, and that
in standard deviation from 0.3 to 16 mW/m*.
Our values are 4-29%, with an average of 15%,
lower than those estimated by Correia et al.
(1990). The differences are due to our porosity
corrections. Therefore, on average, surface heat
flow densities estimated without correcting for
porosity may be 15% higher than porosity-corrected estimates.
Another observation is that for the seventeen
wells, the difference between the maximum and
the minimum surface heat flow density is 33.1
mW/m2 (thermal resistance method) and 23
mW/m2 (inverse method). The difference of the
maximum and minimum heat flow values calculated by Correia et al. (1990) is 26 mW/m2. We
tend to accept 23 mW/m2 as the representative
basinal difference in surface heat flow density
and it is quite near that which is based on the
estimations by Correia et al. (1990). The cause of
this basinal variation is out of the scope of this
study, as we do not have enough data to probe
Table 5
Surface heat flow density (mW/m’)
further. It is generally believed that a basin’s
thermal structure is influenced by heat convection and that heat flow density signals vary due to
basinal fluid motion. Correia et al. (1990) described an increase in heat flow density with the
depth of the deepest BHT in the northern portion of the Jeanne d’Arc basin and a decrease
with the depth of the deepest BHT in the southern part. They construed from this that the basin
is not in thermal equilibrium and is disturbed by
heat transpo~ation through convection, i.e., currently active fluid flow with a downward component in the north and an upward component in
the south. Our formation thermal conductivity
profiles, corrected for porosity effects, pose a
problem with this interpretation. The presence in
the basin of a strongly overpressured zone, characterized by low thermal conductivities (Fig. 51,
implies by definition separate fluid zones. Furthermore, we noticed that the largest difference
in heat flow density is in W-E direction (Rankin
M-36 verse two Whiterose wells). Therefore,
large-scale N-S basinal fluid movements described by Correia et al. (1990) is unlikely.
Alternative reasons for the variation in surface
heat flow densities might be the combined effects
from two methods used in this study and a comparison with those from Correia et al. (1990)
Well name
Thermal resistive
method
Inverse method
Correia et
al. (1991)
Differences a
Ben Nevis I-45
Beothuk M-05
Flying Foam I- 13
Gambo N-70
Hebron 1-13
Hibernia C-96
Hibernia K-14
Hibernia K-18
Hibernia O-35
Mercury K-76
Nautilus C-92
Rankin M-36
South Tempest G-88
Terra Nova K-08
W. Flying Foam L-23
Whiterose J-49
Whiterose N-22
48.3
53.1
50.4
47.4
47.0
49.4
41.5
51.3
51.2
43.7
50.9
33.4
51.1
50.7
57.6
66.5
60.7
51.7 + 10.0
56.2 * 10.9
52.6 * 16.7
51.3 f 15.0
50.0 + 7.4
53.5 f 9.8
47.4 * 10.5
52.0 + 8.4
50.4 + 10.0
45.6 i 10.9
50.7 +_ 9.6
37.8 + 8.1
58.9 * 10.5
52.7 rt 9.1
57.9 a 18.5
60.7 * 13.3
60.7 f 11.9
57 * 22
65+ 17
64& 17
48+ ?
58+ 16
64k 17
50 + 16
57 + 20
71 +23
58 + 15
61+ 21
49 f 12
70 f 27
55* ?
70 + 28
75 + 24
72 + 23
5.3 (12)
8.8 (6.1)
11.4 (0.3)
- 3.3 (?)
8 (8.6)
10.5 (7.2)
2.6 (5.5)
5 (11.6)
20.6 (13)
12.4 (4.1)
10.3 (11.4)
11.2 (3.9)
11.1 (16.5)
2.3 (7)
12.1 (9.5)
14.3 (10.7)
11.3 (11.1)
a Difference between the surface heat flow density values from the inverse method and those estimated by Correia et al. (1990).
The first value is difference in surface heat flow density and the second in the brackets is that in S.D.
Temperature
(degC)
zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
0
20
40
60
80
100 120 140 zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
u
IEBRON
j
r
0.5
i
I
*
I
c
I
/
t
J 140 120 100 80
161
5
’
I
60
’
I
40
’
I
r
20
Time(Ma)
Maturity(%Ro)
b
Fig. 7. (a) Plots of model-predicted maturity/temperature/depth
profiles for the Hebron I-13 well using two sets of thermal
conductivity assumptions, The crosses are measured values for comparison. (b) ~odei-predicted
maturity/ time curves for the
Egret source rock using the two conductivity sets. Curve I is based on thermal conductivities from literature. Curve 2 is based on
thermal conductivities from this study.
of some local factors (e.g., variations in basement
heat flow flux, small-scale thermal convection and
radiogenic heat production) and errors in data.
The effect of each factor may not be very substantial. For example, we have neglected radiogenie heat production in calculation, as there is
no published data on heat production in the
Jeanne d’Arc Basin. According to Keen and Lewis
(1982), the measured heat production in the adjacent continental margin off Nova Scotia ranges
from 0.5 to 1.6 pW,/m3. If these values are applicable in the Jeanne d’Arc Basin, we estimated
that in a drilled interval of 5000 m, containing
50% shale/dlaystone of high radiogenic heat production value (1.6 pW/m3>, the increase of heat
flow density due to heat production in shale/
claystone would be 4 mW/m*. That value alone
can not explain why the heat flow at Rankin
M-36 is about 23 mW/m* lower than at the
Whiterose wells even if we assume that the former had no shale package.
6. Sensitivity of basin modelling to thermal conductivity estimation
Part of the impetus behind our examination of
this basin’s thermal structure is the need to constrain the thermal assumptions necessary to derive information on the timing of hydrocarbon
generation in the region,
This involves translating sediment compaction
histories into time temperature info~ation
by
viewing the rock column as a series of stacked
beds, each with its own thickness and thermal
conductivity (which varies through time with com-
ual lithological layers ( 2 3 m) in the Jeanne d’Arc
Basin, varies from 1.3 to 2.8 W m- ‘K- ‘. At each
well location, in-situ thermal conductivity is controlled by lithological composition. Porosity is also
an important control on the in-situ thermal conductivity. In intervals with higher porosity value,
in-situ thermal conductivity corrected for porosity
can be only about 50% (or lower) of the matrix
thermal conductivity.
(2) The formation thermal conductivity in the
Jeanne d’Arc Basin is heterogeneous.
Shale/
claystone formations are more heterogeneous
than sandstone and sandstone/ siltstone formations. The range of formation thermal conductivity is 1.46-2.65 W m-‘K-l. From the South Mara
Unit down to the Catalina Formation, the formation thermal conductivity generally increases with
depth. However, the formations below do not
follow this trend.
(3) There are at least two reasons for the
heterogeneity of the thermal properties in this
basin. The first is lateral variation in lithological
composition caused by facies change and the second is the development of overpressure zones in
some areas of the basin. The heterogeneity in
sandstone/ siltstone formation may be mainly explained by the first reason. The heterogeneity of
shalc/claystone
fo~ation
is because of both faties change and overpressure zone and the latter
may have played a more important role.
7. zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
Summary
(4) The surface heat flow density in the studied
seventeen wells varies from 37.8 to 60.7 mW/m2,
Understanding
a basin’s thermal structure,
according to the results from the inverse method.
particularly the factors that control it, provide
The basinal average of the surface heat flow
important constraints on models that attempt to
density is 51 mW/m”.
predict the basin’s thermal evolution through ge(5) This study indicates that on average surface
ological times. This in turn influences the quality
heat flow density estimated without porosity corand accuracy of a range of models that require
rection may be 15% higher than the one with
this information as an aid to understanding for
porosity correction.
example the basin’s rifting history or the genera(6) In basin modelhng, different thermal contion, expulsion and migration of hydrocarbons
ductivity estimates are not particularly significant
within the basin. Our study of the basin’s thermal
for present-day predictions but clearly have substructure is based on porosity-corrected
thermal
stantial impact on modelled timing of maturation,
conductivity and used the inversion method to
calculate the surface heat fIow density for seventeen wells. The major results and implications are
as the follows:
This study was financially supported by a
(1) The in-situ thermal conductivity of individCanadian Government Laboratory Visiting Fel-
paction). With a constant heat flow through time
(see Williamson, 1992), surface temperature
(which can vary through time) and assuming
steady-state thermal equilibrium, it is possible to
calculate temperature-depth
profiles for each
time increment from the simple formulation expressed earlier. This is then used as a basis for
computing levels of organic maturity versus depth
at present-day or through time for a specific
datum. The level of sensitivity of these maturity
predictions to assumed thermal conductivities is
demonstrated in Fig. 7. Curve 1 represents maturity predictions based on matrix conductivity values cited in the literature (e.g. Clark, 1966; Gretener, 1981) and based on investigations in other
areas. Curve 2 represents maturity predictions
using conductivities derived in this study, i.e.,
correcting for compaction and based on measurements performed on sediments from the Jeanne
d’Arc Basin, The differences in these two assumption sets are not particularly significant for
present-day predictions but clearly have implications for modelled timing of maturation. The
modelling indicates that the source rock reached
a certain maturity level about 50 Ma earlier with
the thermal conductivities from this study than
with the matrix thermal conductivi~ values from
literature.
192
Z. zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
Huang, M A. W illiamson / Tectonophy sics 233 (1994) 177- l 92
lowship to Z, Huang at the Atlantic Geoscience
Centre of the Bedford Institute Oceanography,
Dartmouth, Nova Scotia. Our study has benefited
from fruitful discussion with Serge Levesque. Kelin Wang kindly provided the program for calculating heat flow density with the inversion method.
We thank Serge Levesque, Sonya Dehler and
John Shimeld for their critical reading of our
manuscript. Drs. Cermark and Mareschal are
thanked for their reviewing. This is Geological
Survey of Canada contribution No. 12092.
References
Allen, P.A. and Allen, J.R., 1990. Basin Analysis; Principles
and Applications. Blackwell, Oxford, 451 pp.
Beach, R.D.W., 1985. Sedimentary heat flow in Alberta. M.Sc.
Thesis, Univ. of Alberta.
Bell, J.S. and Campbell, G.R., 1990. Petroleum resources. In:
M.J. Keen and G.L. Williams (Editors), Geology of the
Continental Margin of Eastern Canada. pp. 679-720.
Brereton, N.R., 1992. Physical property relationship from Site
765 and 766. Proc. ODP, Sci. Results, 123: 453-468.
Brigaud, F., Chapman, D.S. and Le Douaran, S., 1990. Estimating thermal conductivity in sedimentary basin using
lithologic data and geophysical well logs. Bull. Am. Assoc.
Pet. Geol., 74: 1459-1477.
Bullard, E.C., 1939. Heat flow in South Africa. Proc. R. Sot.
London, Ser. A, 173: 474-502.
Chapman, D.S., Keho, T., Bauer, M. and Picard, M.D., 1984.
Heat flow in the Uinta basin determined from bottom hole
temperature (BHT) data. Geophysics, 49: 453-466.
Clark, S.P. (Editor, 1966. Handbook of Physical Constants.
Geol. Sot. Am. Mem. 97, 587 pp.
Correia, A., Jones, F.W. and Fricker, A. 1990. Terrestrial
heat-flow density estimates for the Jeanne d’Arc Basin,
offshore eastern Canada. Geophysics, 55: 1625-1633.
Della Vedova, B. and Von Herzen, R.P., 1987. Geothermal
heat flux at the COST B-2 and B-3 Wells, U.S. Atlantic
Continental Margin. Woods Hole Oceanogr. Inst., Tech.
Rep. WHOI-87-27.
Doveton, J.H., 1986. Log Analysis of Subsurface Geology.
Wiley, New York, NY, 273 pp.
Enachescu, M.E., 1987. Tectonic and structural framework of
the Northeast Newfoundland continental margin. In: C.
Beaumont and A.J. Tankard (Editors), Sedimentary Basins
and Basin-forming Mechanisms. Can. Sot. Pet. Geol.
Mem., 12: 117-146.
Evans, T.R. and Coleman, N.C., 1975. North Sea geothermal
gradients. Nature, 247: 28-30.
Grant, A.C., McAlpine, K.D. and Wade, J.A., 1986. The
continental margin of Eastern Canada: geological frame-
work and petroleum potential. In: M.T. Halbouty (Editor),
Future Petroleum Provinces of the World. Am. Assoc. Pet.
Geol., Mem., 40: 177-205.
Gretener, P.E., 1981. Geothermics: Using temperature in
Hydrocarbon Exploration. Am. Assoc. Pet. Geol., Educ.
Course Note Ser. 17, 170 pp.
Hagarty, K.A., Weissel, J.K. and Mutter, J.C., 1988. Subsidence history of Australia’s southern margin: constraints
on basin models. Bull. Am. Assoc. Pet. Geol., 72: 615-633.
Hyndman, R.D., Jessop, A.M., Judge, A.S. and Rankin, D.S.,
1979. Heatflow in the Maritime provinces of Canada. Can.
J. Earth Sci., 16: 1154-1165.
Kappelmeyer, 0. and Haenel, R., 1974. Geothermics with
special reference to application. Geoexploration Monogr.
Ser. 1, No. 4, Borntraeger, Berlin-Stuttgart, 238 pp.
Keen, C.E. and Lewis, T., 1982. Measured radiogenic heat
production in sediments from continental margin of Eastern North America: implications for petroleum generation. Bull. Am. Assoc. Pet. Geol., 66: 1402-1407.
McAlpine, K.D., 1990. Mesozoic stratigraphy, sedimentary
evolution, and petroleum potential of the Jeanne d’Arc
Basin, Grand Banks of Newfoundland. Geol. Surv. Can.
Pap. 89-17.
Raiga-Clemenceau, J., Martin, J.P. and Nicoletis, S., 1988.
The concept of acoustic formation factor for more accurate porosity determination from sonic transit time data.
Log Anal., 29: 54-60.
Raymer, L.L., Hunt, E.R. and Gardner., J.S., 1980. An improved sonic transit time-to-porosity transform. SPWLA
Trans., 21st Annu. Log. Symp. Pap., pp. 8-11.
Schlumberger, 1989. Log Interpretation
Principles/Applications. Schlumberger, 198 pp.
Tarantola, A. and Valette, B., 1982. Generalized nonlinear
inverse problems solved using the least squares criterion.
Rev. Geophys. Space Phys., 20: 219-232.
Vasseur, G., Lucazeau, F. and Bayer, R., 1985. The problem
of heat flow density determination from inaccurate data.
Tectonophysics, 121: 25-34.
Wang, K., 1992. Estimation of ground surface temperatures
from borehole temperature data. J. Geophys. Res., 97:
2095-2106.
Williamson, M.A., 1992. Subsidence, compaction, thermal and
maturation history of the Egret Member source rock of
the Jeanne d’Arc Basin. Bull. Can. Pet. Geol., 40: 136-150.
Woodside, W. and Messmer, J.H., 1961a. Thermal conductivity of porous media, I. unconsolidated sands. J. Appl.
Phys., 32: 1688-1699.
Woodside, W. and Messmer, J.H., 1961b. Thermal conductivity of porous media, II. consolidated rocks. J. Appl. Phys.,
32: 1699-1706.
Wyllie, M.R.J., Gregory, A.R. and Gardner, L.W., 1956. Elastic wave velocities in heterogeneous and porous media.
Geophysics, 21: 41-70.
Wyllie, M.R.J., Gregory, A.R. and Gardner, G.H.F., 1958. An
experimental investigation of factors affecting elastic wave
velocities in porous media. Geophysics, 23: 459-493.