Subsurface temperature from seismic reflection 7 data : application to the post break up sequence 8 offshore Namibia 9

20 Accurate estimations of present-day subsurface temperatures are of critical importance to 21 the energy industry, in particular with regards to geothermal energy and petroleum 22 exploration. This paper uses seismic reflection observations of bottom-simulating reflections 23 and subsurface velocities coupled with an empirical velocity to thermal conductivity 24 transform to estimate subsurface temperature in a process dubbed reflection seismic 25

This manuscript is a preprint and has been submitted to AAPG Bulletin. This manuscript has 1 not undergone peer review. Subsequent versions of this manuscript may have different 2 content. If accepted, the final accepted version of this manuscript will be available via the 3 'Peer-reviewed Publication' DOI link on the right-hand side of this webpage. Please feel free 4 to contact any of the authors directly to comment on the manuscript. We welcome any 5 feedback! 6 Introduction 33 Subsurface temperature is a key parameter in subsurface energy extraction from petroleum-34 and geothermal systems (Harper, 1971; Thompson, 1979;Hunt, 1984;Bonté et al., 2012). 35 Accurate estimations of present-day subsurface temperatures are thus of critical importance 36 to the energy industry. In frontier areas, petroleum source rock maturity is a key uncertainty 37 and without access to bottom hole temperature (BHT) readings from boreholes, source rock 38 characterisation is reliant on estimation and extrapolation. A crucial component of 39 understanding the subsurface temperature field is how heat is transferred (i.e. heat flow and 40 thermal conductivity) (Sclater et al., 1980). Often there is limited understanding of the 41 variation in thermal conductivity both vertically and laterally in the subsurface domain due to 42 the difficulty collecting such data. Prior to entry into frontier basins, it is advantageous to 43 determine the source rock deliverability in the area. This is primarily controlled by two factors 44 namely source rock type and temperature history. Along many passive margins the source 45 rocks are at maximum burial depth and thus maximum temperature at the present day. 46 Present day subsurface temperatures are typically acquired from temperature probe 47 measurements that have been acquired through borehole drilling (Fuchs and Balling, 2016). 48 To estimate heat flow and temperature information in adjacent areas often involves the 49 deployment of multiple seafloor probes prior to the drilling stage (Davis et al., 2003). This 50 then requires the extrapolation of data from nearby boreholes using structural and 51 stratigraphic models (Davies and Davies, 2010). This methodology underutilises the available 52 seismic datasets that are often acquired during the early stages of exploration in frontier 53 basins. This paper presents a reflection seismic thermometry workflow for using seismic 54 reflection data to estimate subsurface temperature before drilling and applies this to a 55 frontier exploration setting. 56 It has been shown that gas hydrate identification on seismic reflection data through detection 57 of bottom simulating reflectors (BSRs) at the base of gas hydrate stability zone (GHSZ), can be 58 used for geothermal gradient estimation (Yamano et al., 1982 (Bagguley and Prosser,  71 1999; Karner and Driscoll, 1999;Schmidt, 2004). The margin offshore Namibia is characterised 72 as having characteristics of both a volcanic passive margin and a non-volcanic margin end 73 members (Light et al., 1993;Gladczenko et al., 1998;Bauer et al., 2000). Asymmetric rifting 74 has resulted in significant variability in the sedimentary and subsidence history between the 75 conjugate margins with this reflected in the nature of hydrocarbons discovered in these areas 76 Basin multiple features such as seaward dipping reflectors (SDRs) and mass transport deposits 82 (MTDs) can be observed (Torsvik et al., 2009). There are however no seamounts within 50 km 83 (~31.1 mi) from the study area, the outer limit for hydrothermal systems to extend from such 84 a system (Sclater et al., 1980;Hasterok et al., 2011). 85 The neighbouring Orange and Walvis Basins both possess working petroleum systems with 86 gas condensate discovered in the Kudu wells and oil in the Wingat well respectively (Fig. 1) 87 (Intawong et al., 2015). The key source interval relevant for the Lüderitz Basin is the Aptian 88 age Kudu shale (Fig. 2d). The source maturity of this interval in the Lüderitz Basin is as yet 89 uncertain. 90  (Frost and Shive, 1986). Thus, 158 the depth corresponding to a lack of magnetism is likely at temperatures in excess of the Curie 159 point of magnetite or a result of compositional changes leading to magnetite poor rocks at 160 depth (Beardsmore and Cull, 2001). However, this method simply defines a solitary 161 subsurface isotherm over a great depth interval from the seabed, thus making any linear 162 geothermal gradient calculated greatly approximated. Furthermore, for high resolution 163 estimation of the Curie isotherm depth, regional scale gravity data would be required. 164 Global thermal data coverage may also suffer from spatial bias, as there tends to be a greater 165 interest for scientists in areas of higher heat flow resulting in a greater concentration of data, 166 with another potential driving factor being interest in areas with geothermal energy 167 application (Davies, 2013). 168   The temperature estimation workflow utilised in this study is outlined in Fig. 4 and described 211 below. The entirety of the workflow has been developed and tested using commercial 212 software developed for use in the petroleum industry (Schlumberger Petrel). 213 The gas hydrate stability field can be utilised to estimate the temperature at the base of the 214 zone of stable gas hydrates, demarcated on seismic by a BSR (Dickens and Quinby-Hunt, 1994; 215 Sloan et al., 1998;Lu and Sultan, 2008). This in turn allows a shallow geothermal gradient 216 across the GHSZ to be estimated and surface heat flow to be estimated (Minshull, 2011;217 Priyanto, 2018). The stability conditions are controlled in part by the geochemical properties 218 of the fluids available to form clathrate hydrates, which in frontier settings with limited 219 ground truthing are generally assumed to be average salinity (33.5 ‰) seawater and pure 220 methane (Sloan et al., 1998). Pore fluid pressure conditions are generally assumed to be 221 hydrostatic, equivalent to 0.0101 MPa m -1 (~0.446 psi ft -1 ). The following relationship (Eq. 1) 222 as defined by (Dickens and Quinby-Hunt, 1994) describes methane hydrate stability: 223 where TBSR is temperature at the base of GHSZ (K); and P is the corresponding pressure (MPa). 225 Assuming hydrostatic pressure at the BSR depth, temperature at the base of the hydrate 226 stability zone has been established: 227 Where TBSR is the temperature at GHSZ (°C); is density (kg m-3 ) (of seawater); g is 229 acceleration due to gravity (m s-2 ) and ZBSR is the depth (m) of the BSR. It must be noted that 230 this is a minimum temperature estimate based on assumed stability field conditions (Dickens, 231 2001 where TSEABED is the modelled hydrothermal gradient temperature (°C) and Z is seabed depth 241 (m). 242 Given both TSEABED, ZSEABED and TBSR, ZBSR at any geographical locality, then the geothermal 243 gradient (dT/dZ) across the GHSZ at that locality is given by the following relationship. 244 Where dT/dZ is geothermal gradient (°C km -1 ); TBSR is temperature at BSR (°C); TSEABED is seabed 246 temperature (°C); ZBSR is depth of BSR (km); ZSEABED is seafloor depth (km). 247 Alongside thermal gradient, two key thermal properties are the heat flow and thermal 248 conductivity. 249 Where Q is heat flow (mWm -2 ); k is thermal conductivity (W m -1 K -1 ) (see Section 3.1) and 251 dT/dZ is geothermal gradient (°C km -1 ). 252 This however does not consider the thermal conductivity structure of the subsurface and how 257 it might be possible to utilise seismic reflection velocity data to do so. 258 Thermal conductivity estimation 259 Thermal conductivity is a measure of how well heat is conducted through a material (Gu et 260 al., 2017). Difficulty associated with measuring thermal conductivity in boreholes arise from 261 poor contact between the measuring tool and the borehole wall (Horai, 1982). Thus, 262 considerable attention has been devoted to determining methods for estimating thermal 263 conductivity through more easily acquired secondary data such as seismic velocity 264 measurements. Experimental studies have shown that primary controls on thermal 265 conductivity include mineral composition, porosity and fractures (Gegenhuber and Schoen, 266 2012). Seismic wave velocity is also largely controlled by the same factors. Early work by 267 (Horai, 1982)  collected under similar parameters. Most studies measured thermal conductivity using the 283 optical scanning method (Popov et al., 1999). There are fewer instances in the source datasets 284 of the use of the divided bar method of measuring thermal conductivity (Hyndman and 285 Jolivet, 1976;Evans, 1977). Only wet samples from these studies were used as our case study 286 is in deep water and thus fully saturated with water, gas and/or gas hydrate. In dry samples, 287 the contribution to thermal conductivity arising from lithological heterogeneities (matrix 288 properties) can be masked by the stronger influence of porosity (Hartmann et al., 2005). In Fractures are known to reduce both P wave velocities and thermal conductivity (Zamora et 298 al., 1993). 299 A regression through the filtered experimental data points taken from the aforementioned 300 studies gives the following empirical relationship for thermal conductivity: 301 Where kV is thermal conductivity from velocity (W m -1 K -1 ) and VP is P wave velocity (m s -1 ). 303 Seismic P wave velocity within the area is converted to thermal conductivity (kV) using the 310 thermal conductivity relationship (Eq. 6), with velocity averaged down to the depth of the 311 BSR, ZBSR. The variation in thermal conductivity with depth can be overlain on a 3D seismic 312 reflection dataset in this manner. Using ZBSR, determined on reflection seismic data, the 313 hydrate stability field can be utilised to compute the temperature at this phase boundary for 314 the base of the GHSZ (using Eq. 2). Temperature at the seabed is known from the 315 hydrothermal gradient (given by Eq. 3). A shallow geothermal gradient may thus be computed 316 between seabed and BSR (Eq. 4). As thermal conductivity has been derived from acoustic 317 velocity data, and with shallow geothermal gradient also available, it becomes possible to 318 reapply Fourier's Law (Eq. 5) to derive heat flow for this area through inverse modelling. 319 Estimating the shallow geotherm and heat flow along the full extent of a BSR helps eliminate 320 the bias in heat flow distribution from direct measurements taken at discrete locations 321 (Shankar and Riedel, 2013). This BSR derived heat flow proxy is used in conjunction with the 322 bulk thermal conductivity volume to generate a volume of average geothermal gradient for 323 the bulk volume (rearranging Eq. 5). 324 Temperature below the seafloor can be summarised as being a function of the depth below 325 the seafloor and the average geothermal gradient. It follows that an estimate of temperature 326 may be arrived at through this simple relationship where the temperature at any given depth 327 point is given by multiplying the average geothermal gradient against the depth to that point: 328 where T is predicted temperature (°C); TSEABED is the temperature at seabed (°C); dT/dZ is the 330 average geothermal gradient (°C km -1 ); and ZSUBSURFACE is the subsurface depth (km). Seabed 331 temperature is added to account for the effect of the hydrothermal gradient on the 332 subsurface temperatures. 333 As the average geothermal gradient is only valid for the subsurface and due to the seismic 334 input volume containing the water column it becomes necessary to negate the latter. Without 335 flattening the volume to the seabed, it is instead possible to use the seabed depth map to 336 derive a depth volume relative to seabed depth. 337

374
The BSR observed in the area has been mapped across the NW and SW quadrants of the 3D 375 reflection seismic coverage (Fig. 6). Though the full extent of the visible BSR was mapped, only 376 the extent corresponding to the highest confidence seismic picks are displayed as the clarity 377 of the BSR degrades towards the edges. This should preclude any resulting anomalous 378 artefacts and edge effects. It is this high confidence extent of the BSR that is referred to in the 379 following sections unless otherwise specified. The BSRs are found to have opposite seismic 380 reflection polarity to the seabed reflection indicating the likelihood of gas hydrate above free 381 gas (Kretschmer et al., 2015). Though there is no record of hydrates from ODP Site 1084, high 382 amplitude reflections are observed to occur in close proximity below the BSR (Fig. 2a), 383 characteristic of the presence of trapped gas. Temperature at BSR depth and the phase 384 relationship used to determine this is shown in Fig. 6. 385

Figure 6 386
Neither the exploration well nor ODP Site 1084 fall within the bounds of the thermal model. 387 As a result, direct calibration is not possible. However well 2513/8-1 contains BHT information 388 that may provide some calibration for the predicted results. Pseudo-wells provide a means of 389 simulating 2513/8-1 at a comparable location along strike (Fig. 1a). P1 is projected into the 390 study area following bathymetric contours as close as possible along strike from 2513/8-1, to 391 maintain structural parity. BHT recordings typically are lower than actual formation 392 temperature due to cooling effect of circulating fluids in a borehole and thus they must be 393 corrected (Deming, 1989). There are insufficient points for a Horner correction (Horner, 1951;394 Bonté et al., 2012) to be applied and hence a rudimentary correction is made for time since 395 circulation (see https://www.zetaware.com/utilities/bht/timesince.html first accessed 396 August 2018). The predicted temperatures are between 17 and 26% higher than the corrected 397 BHT (Fig. 7a). 398

Figure 7 399
On seismic data it was evident that there is a deeply incised canyon like structure trending NE 400 -SW that can be seen in the north-eastern most extent of the seismic volume (Wanke and 401 Toirac-proenza, 2018). This corresponds to the location of P1, which is seen to intersect the 402 channel fill structures of this canyon. It becomes evident then that though P1 was projected 403 into the seismic volume maintaining bathymetric parity, in the subsurface, due to the 404 occurrence of this channel like geometry, it is not possible to maintain stratigraphic parity to 405 2513/8-1. This is surmised to be the primary factor for the misfit with BHT seen. 406 Further pseudo-wells (T1 -3) were modelled to examine the change in thermal profile moving 407 from the proximal section to the distal part of the study area. The results (Fig. 7) display what 408 the thermal profile in these boreholes would be like should a typical geothermal gradient of 409 30 °C km -1 or 40 °C km -1 was applied linearly from seabed. The temperature window 410 considered prospective for reservoirs in present day has been referred to as the Golden Zone 411 However, in the proximal end corresponding to minimal Tertiary cover, there is observed the 421 greatest divergence between isotherms in Mesozoic sediment. Temperature for the Aptian 422 'Kudu shale' source rock in the region has also been mapped (Fig. 8). 423

Figure 8 424
Below both BSRs, but particularly the northern BSR (Fig. 6), the effects of gas blanking were 425 observed in the seismic reflection data. An average interval velocity extraction reveals 426 anomalously low values within this area (Fig. 8) process. For this work, with a lack of well data for ground truthing, the temperature 437 estimation bounds for 95% confidence were used to give an idea of the range within which 438 the estimates can vary. It is important to note the impact of variability in input factors for the 439 component steps. For example, results from the Blake Ridge show that actual temperatures 440 at the BSR depth could be between 0.5 -2.9 °C (32.9 -37.22 °F) lower than the temperature 441 predicted by the hydrate phase relationship for that particular depth and pressure (Wood and 442 Ruppel, 2000). This implies that a significant source of uncertainty in the thermal modelling The quality of the initial velocity model is another source of uncertainty. As thermal 457 conductivity is derived from it using a direct empirical relationship, any anomalies in the 458 existing velocity model or velocity data will be translated into the derived properties. From 459 the low spread of RMS interval velocities for the GHSZ it is apparent the application of a 460 default 1500 m s -1 (~4921 ft s -1 ) velocity above seabed during the velocity model building stage 461 results in a heavily smoothed velocity model. This is expected to be reflected in the nature of 462 the temperature profile generated using velocities as input. 463 In the absence of a reliable heat flow recording for this area, a BSR derived heat flow proxy 464 has been used. This is a shallow heat flow as it uses an average velocity derived thermal 465 conductivity and geothermal gradient valid within the GHSZ (Eq. 5). Unlike in traditional basin 466 modelling the radiogenic heat production of the rock column has not been integrated. 467 Instead, this solitary heat flow proxy has been used to condition the model for an average 468 geothermal gradient. Though hydrothermal fluid circulation in the subsurface can also greatly 469 alter heat flow, both vertically and laterally, the study area is likely to be minimally impacted 470 in this regard. As the study area is sufficiently distant from a neighbouring seamount to negate 471 the convective and advective heat flow impact of hydrothermal fluid circulation, heat 472 transport in this area is predominantly conductive. Therefore, the assumption is of limited 473 lateral heat flow variability, which is backed by the BSR-derived thermal gradients and 474 derivative heat flow estimates. In a separate case study covering the data rich North Sea, it 475 has been shown that the reflection seismic thermometric process can be conducted 476 successfully using laterally varying shallow heat flow as an input (Sarkar, 2020). of this Barremian structure (Fig. 8a), the base of the overlying Kudu shale source rock 492 immediately above would therefore lie in the gas generation window (Bjørlykke et al., 1989). 493 This is consistent with the nearby Kudu fields which produce gas condensate from the same 494 Aptian source interval. The results would suggest then that the Lüderitz Basin has an improved 495 prospectivity outlook with the potential for a working hydrocarbon system with gas charged 496 reservoirs. 497 It was possible in this study to estimate present day temperature at key subsurface target 498 depths in a frontier setting in the absence of any substantive well control. The workflow 499 presented would enable seismic operators to utilise the data libraries of seismic reflection 500 and velocity data available to them to generate present day estimations of subsurface 501 temperature in a non-invasive manner, prior to an expensive drilling campaign. It is hoped 502 that this would help streamline petroleum systems analysis and provide an additional dataset 503 for basin modellers to use with the aim of decreasing the uncertainty with which frontier 504 regions are explored. 505

506
The model proposed in this study is a simple and robust methodology for estimation of 507 present-day subsurface temperature in frontier areas lacking borehole control for 508 temperatures. It makes use of readily available seismic reflection and velocity data in a 509 workflow developed on an industry standard software suite. It highlights how existing 510 workflows for BSR derived heat flow may be combined with existing experimental thermal 511 conductivity and velocity data for various lithologies to develop an empirical transform that 512 may be applied to seismic velocity models. Given thermal conductivity and P wave velocity 513 have sensitivity to similar parameters, this methodology would allow the user to examine the 514 vertical and lateral variability in thermal properties in a frontier basin especially when high-515 quality pre-SDM and FWI velocity models are available. The results of the case study 516 documented here suggest that the main prospect lies just below the golden zone and that 517 sources rocks are in the generative window. 518 References 519 Allen, P., and J. Allen, 2013, Basin analysis: principles and application to petroleum play 520 assessment.  seismic thermometry workflow first presented in (Sarkar, 2020 system, used to compute temperature at the phase boundary (Fig. 6c). A synthetic 846 hydrothermal gradient is shown, computed using the annualised mean temperature data 847 points from the 1 degree resolution dataset of the WOA (Locarnini et al., 2013). The hydrate 848