Thermal modeling and heat flow density interpretation of the onshore Northwest Java Basin, Indonesia

Studies pertaining to heat flow density distribution in Indonesian tertiary sedimentary basins have been conducted intermittently throughout more than two decades. Even so, all of these studies have always concentrated on the basic compilations of thermal gradient, thermal conductivity, and terrestrial heat flow datasets or the study of basin thermal history rather than on the investigations into their present thermal regime. In addition, syntheses on the possible causes to the observed heat flow density distribution in these basins have been largely determined only on the basis of its correlation with geological and geophysical evidences. One basin experiencing this particular problem is the onshore Northwest Java Basin, where knowledge of the present-day subsurface temperature structure is virtually absent and interpretation of heat flow distribution has been mainly qualitative. For these reasons, in this study we have attempted to quantitatively characterize the present subsurface temperature distribution and surface heat flow density of the onshore Northwest Java Basin using a simple numerical thermal modeling approach, based on the most recent heat flow datasets. Taking into account a conductive heat transfer process, the modeling exercise is aimed at deriving a 3-D subsurface temperature distribution. The temperature distribution is translated into the predicted surface heat flow distribution. The modeled temperature field reveals its dependence on the geometrical aspects of the basin’s basement. Comparison between the computed and the observed heat flow densities, taking into account the uncertainties of both, suggests that advective heat transport by groundwater flow may be present, in addition to conduction. The modeling results are thus capable of demonstrating the importance of a quantitative approach in studying the present thermal state of a sedimentary basin.

decade-long time gap since an updated compilation of HFD data was made for sedimentary basins in Indonesia. The most recent compilation was produced in a study by Suryantini (2007). Although the study was carried out for the entire West Java, a more substantial part of it was dedicated to analyzing the HFD distribution in the onshore Northwest Java Basin area (Fig. 1). Nevertheless, the primary goals of Suryantini (2007) were the determination of heat flow values and their use in a thermal history modeling exercise of a single well located at the onshore Northwest Java Basin, rather than the characterization of the entire basin's present-day thermal state. In addition, the interpretation concerning the observed HFD distribution is also limited to spatial correlation on a qualitative basis, i.e., by establishing a relationship between the HFD pattern and geological information and gravity anomaly.
The main components of the thermal regime of a region consist of the present subsurface temperature distribution and the surface HFD. Any standard approaches to assessing these elements in a sedimentary basin would always require that subsurface temperature be estimated first. This is most effectively achieved through quantitative modeling based on the available information of the subsurface thermal property structure and a set of boundary conditions (e.g., Cacace et al. 2013;Noack et al. 2010). A similar study has been performed by Nagao et al. (1995) for the entire Southeast Asia region using a gridded HFD map and applying a one-dimensional thermal modeling approach to each grid. However, local variations in the onshore Northwest Java Basin area could not be resolved due to the applied grid size being larger than the basin itself. On the other hand, Tanaka et al. (1999) constructed a Curie Point Depth map to define the regional thermal structure over East and Southeast Asia, yet another problem was encountered-a remarkably large gap over the western part of Indonesia remains present (cf. Fig. 2 of Tanaka et al. 1999).
Therefore, in this study, we attempt to perform the first numerical modeling study of the 3-D thermal structure of the onshore Northwest Java Basin by taking into account a Fig. 1 Map showing the distribution of Indonesian hydrocarbon-producing tertiary sedimentary basins, showing the location of the offshore and onshore parts of the Northwest Java Basin. Numbers correspond to the basin names. 1 North Sumatra Basin, 2 Central Sumatra Basin, 3 South Sumatra Basin, 4 West Natuna Basin, 5 Sunda Basin, 6 Northwest Java Basin, 7 North Central Java Basin,8 North East Java Basin,9 Barito Basin,10 Kutai Basin,11 Tarakan Basin,12 Banggai Basin,13 Bone Basin,14 Seram Basin,15 Salawati Basin,16 Berau Basin,17 Bintuni Basin,respectively conductive heat transfer. The procedure is based on heat flow density dataset gathered by Suryantini (2007) and a number of publicly available data. The influence of lithology and basement geometry on the thermal structure of the basin is discussed by examining the spatial variation of the modeled temperature distribution. We also compare the observed HFD values with those calculated from the modeled temperature distribution to aid in interpreting the spatial variation of their differences. These differences are then interpreted by attributing them to other subsurface heat transfer-related phenomena, such as advective redistribution of heat by groundwater flow (e.g., Kilty and Chapman 1980) and refraction of conductive heat (e.g., Beardsmore 2005). We demonstrate that the modeling results enable the identification of prominent thermal signatures in relevance to the geological elements of the onshore Northwest Java Basin and provide insight into its present-day thermal regime.

Methods
As mentioned previously, in this study, we aim to generate a conductive subsurface temperature distribution of the onshore Northwest Java Basin, and thus its conductive heat flow density distribution through a numerical modeling procedure. Thermal parameters and geological data are collated based on previous studies concerning their compilation and earlier interpretation. The latter data are utilized to produce a basin fill and basement geometry model to be assigned thermal conductivity values, to construct a subsurface 3-D thermal conductivity model. The numerical modeling is based on a finite difference approximation of the heat conduction equation, to which the thermal conductivity model becomes the primary input. The modeled heat flow density is then compared against the observed one in a similar fashion to that of Cacace et al. (2013), the resulting difference(s) of which is subsequently interpreted in terms of its relationship with the basin's thermal and hydrological regimes. A flowchart is presented in Fig. 2 to provide the overall picture of the research methodology employed in this study. The data, assumptions, procedures, and results relevant to each component of the workflow are elaborated in much more detail in the following sections.

Geological setting of the onshore Northwest Java Basin
The Northwest Java Basin is among several other western Indonesian tertiary sedimentary basins that have been confirmed as hydrocarbon producing (Fig. 1). This has resulted in several studies on its hydrocarbon geology and resource potential, for instance: Kingston (1988), Bishop (2000), Pethe (2013). The onshore Northwest Java Basin is physiographically located on the coastal plain of Jakarta (Suryantini et al. 2006) and encompasses Banten, West Java, and capital city of Jakarta provinces, respectively (Fig. 3). The oldest tectonic event in the basin occurred at the early Oligocene (Fig. 4). It began as a rapid tectonic extension-induced subsidence phase around 36-30 Ma (Koesoemadinata et al. 1994), before ending at the late Oligocene. The cessation of the rifting episode at the Early Miocene was continued by a slower, thermally induced subsidencethe sagging phase (Hall and Morley 2004) about 20-10 Ma-which eventually shifted to compression at the late Miocene-Pliocene, marking the beginning of an uplift phase.
An alternating pattern between several smaller depressions and basement high areas are found to constitute the larger part of the onshore Northwest Java Basin (Bishop 2000;Fig. 5). Systematically, from the west to the east these compartments are Tangerang High, Ciputat Sub-basin, Rengasdengklok High, Pasirputih Sub-basin, Pamanukan High, and Jatibarang Sub-basin, respectively. These compartments are separated from each other by the presence of deep-seated basement faults (Setyowiyoto et al. 2007) of mainly Fig. 3 Digital terrain model of the study area, the central part of which is the onshore Northwest Java Basin area, bounded by the Java Sea and a series of volcanoes to the north and south, respectively. Also shown are several geographic features of this region NE-SW orientations (Hall and Morley 2004;Clure 2005). The northern boundary of the basin is defined by the Java sea as part of the Sunda Shelf, i.e., a shallow marine area situated between the islands of Sumatra, Java, and Borneo (Hall and Morley 2004). A series of quaternary volcanoes and large crescent-shaped thrust faults which run approximately ESE-WNW mark the basin's southern boundary (Fig. 5) that was developed during the Pliocene-Pleistocene transition (Kloosterman 1989). The subsurface lithology of Fig. 4 Generalized stratigraphic column of the onshore Northwest Java Basin area summarized from several references, showing the subsurface lithological composition and important tectonic events that have occurred since the initial formation of the basin. Numbers next to the label of sections of the column refer to the studies from which the information was taken-(1) Bishop (2000), (2) Setyowiyoto et al. (2007), (3) Pethe (2013), (4) Suryantini (2007) the onshore Northwest Java Basin is primarily composed of sandy, shaly, and calcareous rocks, except for deeper sections, which are composed of volcanic rocks underlain by a pre-tertiary metamorphic basement (Fig. 4). Most superficial deposits consist of quaternary alluvium and volcanic products, the proportion of which slightly decreases to the south as the appearance of tertiary and quaternary igneous intrusions with tertiary limestone and sediments becomes more common (Fig. 5).

Heat flow density datasets
The primary data used in this study are based on the most recent compilation of heat flow density, thermal gradient, and thermal conductivity values of the west Java region, each of which was obtained from hydrocarbon wells by Suryantini (2007). Several additional values were also obtained from the collection of open-access world heat flow datasets maintained in the website (http://www.heatflow.und.edu) to make up for a few spatial gaps in the central and western parts of the study area (Fig. 5). The latter data originated from the compilation of thermal gradients by the Indonesian Petroleum Association (IPA) and South East Asia Petroleum Exploration Society (SEAPEX) which were published earlier by Rutherford and Qureshi (1981) and HFD values published by Thamrin (1985).
The principal difference between the two datasets lies in the way thermal conductivity values are reported. In the former dataset, both the averaged and individual formation's thermal conductivity values for every well in the region where available are informed, whereas the latter only addressed the former type of thermal conductivity. Corrected and uncorrected bottom-hole temperature (BHT) and drill-stem test (DST) temperature values were used to calculate thermal gradients by means of least-square regressions (Suryantini 2007), while formation thermal conductivities were either directly measured where core samples are available or, where they are not, were assumed based on local geological similarity between well locations (Suryantini et al. 2006). The individual  Suryantini 2007). The basement depth contour, subsurface deep-seated and near-surface geological structures and subbasin boundaries are overlain on the map. The distribution of heat flow density data points used in this study that were obtained from hydrocarbon wells over the basin is also shown formation thermal conductivities were subsequently converted into vertically averaged conductivities using a harmonic-averaging technique, following that of Nurusman and Subono (1995). Uncertainties surrounding the measured thermal gradients and conductivities are also reported for the Suryantini (2007) dataset, with the error propagation principle (Navidi 2007) being used to compute the final uncertainties of the HFD values. The second dataset does not address measurement uncertainties, so values derived by Thamrin (1986) have to be adopted for each of the measured parameters and combined using the said error propagation theory to yield the HFD uncertainties (Fig. 6a).

General description of heat flow density distribution
The HFD in the onshore Northwest Java Basin features considerably high values (Fig. 6a, b) with an average of 94.05 ± 26.42 mW/m 2 , compared to that of average continental heat flow which is only about 65 mW/m 2 (Stein 1995;Pollack et al. 1993). Although there have been different views pertaining to whether the basins situated behind the Sumatra-Java volcanic arcs are of non back-arc (e.g., Uyeda and Kanamori 1979) or backarc basins (e.g., Koesoemadinata et al. 1994;Suryantini et al. 2006), Matsubayashi and  Nagao (1991) viewed the crust beneath these areas as being a "semi-continent" one, which has HFD similar in characteristics to the basin and Range provinces. The average HFD within this province is actually ~85-90 mW/m 2 (Wisian and Blackwell 2004), slightly lower than the mean HFD of the onshore Northwest Java Basin, where the latter may be affected by a small number of wells with very high HFD (see discussion below and in "Discussion" section).
There are two major ranges into which heat flow densities fall in most parts of the basin, i.e., 75-90 and 90-105 mW/m 2 , respectively ( Fig. 6b). High to extremely high values (105 to > 120 mW/m 2 ) along with the lowest ones (<75 mW/m 2 ) are also distributed in several isolated localities. The higher HFD values are mostly found in the southern part of the basin, the cause of which Suryantini et al. (2006) attributed to the possible location of the boundary between the coastal plain of Jakarta and the Bogor Zone folded units. Although the HFD distribution does not appear to correlate with the superficial lithology of the basin (Fig. 6a), in general, it coincides with the basin geometry and geological structures. The high to anomalously high values, for example, are located in either area of very shallow basement (well SUM244 in the Tangerang High) or close to large faults (e.g., wells JNG-1, JRR-1, JRR-2 southeast of Ciputat Sub-basin, and PJN-P1, PJB-P1 south of Pasirputih Sub-basin), and areas located on both (well RDH-2 in the Rengasdengklok High). On the other hand, the lowest heat flow values (60-75 mW/m 2 ) are located around areas of basement low (e.g., wells SUM245 in the Ciputat Sub-basin, KRW-1 and SKD-1 in the Pasirputih Sub-basin), as well as nearby fault (well TGB-1).
The aforementioned facts suggest that the HFD distribution in the basin might be controlled by both basement topography, i.e., through heat refraction (e.g., Beardsmore 2005;Thakur et al. 2012;Rawling et al. 2013) due to the basement-to-sedimentary fill thermal conductivity contrast and the presence of heat transfer by moving groundwater (Kilty and Chapman 1980) which may recharge or discharge along the basin's bounding faults (Ehlers and Chapman 1999), or both. However, there are also high HFD values that are not directly located in the vicinity of faults nor are they situated over a shallow basement, i.e., wells CLU-5, SNT-1, and PCT-1 (northern Pasirputih Sub-basin, east, and southeast of Pamanukan High, respectively), though Suryantini (2007) interpreted based on gravity studies that several concealed faults possibly exist beneath areas around SNT-1. Moreover, there are not any anomalously high or low heat flow densities in Pamanukan High region where most data are situated near major faults. We will examine this matter of complication in the following sections by the use of numerical modeling of conductive heat transfer in the basin. The modeled temperature field allows for the conductive component of the HFD values to be recovered (Cacace et al. 2013), which in turn enables a subsequent comparison to be made with the observed heat flow values.

Numerical modeling of subsurface temperature distribution
Our numerical modeling procedure is based on several assumptions. First, pertaining to the present thermal state, the basin can be considered to be in thermal steady state since the timing of the last tectonic event. The justification for this assumption is based on the time scale or the thermal time constant (Jaupart and Mareschal 2011;Stuwe 2007): Other authors (e.g., Beardsmore 2004; Stuwe 2007) prefer a different type of proportionality to use as the denominator in Eq. 1a, as that adopted in this study: where l and κ are the characteristic length scale (m) and thermal diffusivity (m 2 s −1 ) of the material (lithology), respectively. The characteristic length scale can be envisaged as the approximate distance that any thermal perturbations may propagate to for a given amount of time. The value of thermal diffusivity is commonly in the order of ~10 −6 m 2 s −1 ; however, it can be made a function of thermal conductivity ( , W m −1 K −1 ), i.e., κ ∼ /(2.3 × 10 6 ) (Beck 1988). The value of τ can thus be used to estimate either the duration of the heat diffusion process or the length of the influenced zone after a thermal event has taken place at a certain time (Stuwe 2007). The latter is of importance in justifying the aforementioned assumption, since we would like to calculate the depth range influenced by a thermal perturbation associated with the last tectonic event, i.e., uplift at the late Miocene-Pliocene. The appropriate τ is thus about 5-6 my BP (Hall and Morley 2004). Taking the average thermal conductivity of the overall Northwest Java Basin to be 1.87 W m −1 K −1 (cf. Table 2 of Matsubayashi and Nagao 1991), κ ∼ 9.7 × 10 −7 m 2 s −1 , the characteristic length scale calculated using Eq. 1b is about 34.11-37.37 km, respectively. This suggests that with respect to the approximate crustal thicknesses of areas within the Sunda Shelf, i.e., ~20 km (Ben-Avraham, 1973), the late Miocene-Pliocene tectonic thermal perturbation must have completely equilibrated through the crust, thus justifying the steady-state assumption.
The second assumption made in this study concerns the mode of heat transfer within the basin. For the purpose of this study, we consider only conduction, since (1) surface manifestations of a hydrothermal system are virtually absent (e.g., warm or hot spring) within the study area, (2) to date, no information on the deep groundwater flow regime in the basin has been obtained, (3) the nature of temperature data from which HFD values were calculated is discrete rather than continuous, thus rendering the detection of groundwater flow and its influence on the subsurface temperature field extremely difficult, and (4) the topography of the area is by itself flat (Fig. 7), thereby reducing the chance of triggering significant and deeply circulating groundwater flow within the basin. In this case, the numerical modeling of temperature distribution may be carried out conveniently by considering only conductive heat transfer (e.g., Beardsmore 2004; Mottaghy et al. 2011).

Governing equation and numerical solution
The equation governing conductive heat transfer through a 3-D, heterogeneous and isotropic, heat-generating solid medium is (Haenel et al. 1988): where ρ, C, and H are each the rock density (kg m −3 ), specific heat capacity (J kg −1 K −1 ), and radiogenic heat production (W m −3 ), respectively. The right-hand term of Eq. 2a can be dropped to yield the equation of steady-state heat conduction for a heterogeneous medium: Since, in a heterogeneous medium, the thermal conductivity is bound to vary spatially (as indicated by the subscripts of each thermal conductivity component), Eq. 2b is non-linear and a closed form analytical solution like that tabulated by Carslaw and Jaeger (1959) for simplified geometries cannot be derived (Nagihara 2003). Instead, the solution is found using a numerical method. In this study, we use a finite difference approximation (Haenel et al. 1988) to solve Eq. 2b for the 3-D conductive temperature distribution, based on a set of material thermal property distribution and boundary conditions (Beardsmore and Cull 2001). For the sake of simplicity, a Gauss-Seidel scheme (Chapra 2011) is adopted to solve for the sets of resulting linear finite difference system of equations (e.g., Kilty and Chapman 1980) of temperature at each grid node, where a fixed number of 1000 iterations is used for each modeling run, ensuring sufficiently small temperature residuals.

Subsurface thermal properties
As the numerical solution to the equation of steady-state heat conduction with heat generation requires the thermal properties distribution of the medium be known, the subsurface 3-D rock thermal conductivity and radiogenic heat production models must be generated. The lithological units involved in the construction of the models are comprised of the five main groups present in the onshore Northwest Java Basin (Fig. 4): Jatibarang, Lower Cibulakan, Upper Cibulakan, Parigi, and Cisubuh, including the basement. Table 1 lists the bulk thermal conductivities assigned to each formation, which are obtained by averaging values of those collated in Suryantini (2007). The thermal conductivities are not corrected for in situ pressure and temperature conditions, because, aside from the lack of data, the effect of sediment compaction and that of in situ temperature may counteract each other, making the overall conductivity more or less constant with depth for each formation (e.g., Pasquale et al. 2012).
To date, both direct measurements and indirect estimations of radiogenic heat production of sedimentary and basement rocks collected from and/or present in Indonesian tertiary basins are evidently absent. Consequently, we must adopt anCross sections extracted from the 3-D conductive temperature established model of the distribution of radiogenic heat generation with depth which has been applied on area(s) with roughly similar geology, and apply it to the entire basin. To this end, we utilize an exponential model of heat generation versus depth: where H 0 and H (z) denote the radiogenic heat generation at the surface and at a particular depth, respectively. Parameter D represents the depth scale for vertical distribution of heat-generating elements. The model was successfully applied to the entire Southeast Asia, covering our region of interest by Nagao et al. (1995), using H 0 = 3.1 μW m −3 and D = 7.8 km. The subsurface geometry of the basin fill, basement, and structures are based on geological cross section of Kingston (1988) as well as a basement depth map utilized by Suryantini (2007). The 3-D thermal conductivity model constructed using values in Table 1 and subsurface geometry is displayed as a pair of cross sections in Fig. 8.

Model domain setup
To apply the finite difference procedure, the 3-D model domain is first discretized into a series of cells and nodes of finite lengths as displayed in Fig. 7. The original region to be modeled is bordered by the solid black line, and we deliberately extend it to the maximum areal extent shown and covered by Fig. 7, by mirroring the edges and sides of the evaluated region to the sides of the extended model domain. This is done to increase the ratio of the total modeled area to the original area of interest (i.e., the outline of the onshore Northwest Java Basin, Fig. 3), so as to suppress the effect of lateral boundary conditions on the modeled temperatures within the latter area (Nagihara 2003;Norden et al. 2008). This practice also has an additional advantage of reducing the complexity of assigning lateral boundary condition values, as the lateral geometry of the total extent of the model domain becomes rectangular. There is insufficient and unreliable information regarding the deep thermal structure of the studied area, even if it is to be derived or inferred from, e.g., other geophysical studies. For instance, as has been addressed in "Background" section, a Curie point depth investigation of SE Asia (Tanaka et al. 1999), which resulted in contours of depth to the 450 °C isotherm, did not cover our region. In addition, a direct thermal modeling study of Nagao et al. (1995), which includes the studied area, utilized a geographic grid size of 5 o × 5 o , larger than the smallest extent of the basin. Owing to these limitations, a constant heat flow value and its uncertainty are  Table 1 employed at the base of the model domain instead of a constant temperature. The basal heat flow and its uncertainty are obtained by subtracting the total radiogenic heat production of the model domain volume from the average surface HFD and its uncertainty (e.g., Limberger et al. 2014). The topography is set to zero meter above sea-level everywhere within the model domain, including the extended area (Fig. 7). Although the initial surface temperature assumed by Suryantini et al. (2006) is 26.67 °C, which conforms to that of the surface air temperature by Hijmans et al. (2005) (http://www.worldclim. org), the average value calculated from the Suryantini (2007) complete HFD dataset is 28 °C. Since some studies have indicated that the ground surface temperature is slightly above that of the surface air (e.g., Smerdon et al. 2006), the latter was imposed as a uniform constant temperature on the top boundary. A complete list of the numerical model domain parameters as well as boundary conditions is given in Table 2. We use a system of finite difference pattern as described by Beardsmore and Cull (2001). In their scheme, temperature is solved at each node, whereas the values of thermal conductivity and radiogenic heat production must be assigned to the cells prior to the temperature calculations. In other words, the discretization approach we adopt in this study makes use of a node-based scheme (Murthy et al. 2006), i.e., the grid system is mesh centered (Faust and Mercer 1980). To explore the uncertainties in the resulting modeled temperature distribution, the numerical modeling is carried out three times; each one is based on the original, maximum, and minimum value of formation thermal conductivity and basal heat flow density.

Synthetic results
To validate the reliability of our finite difference numerical modeling results, we first apply the procedure to a synthetic 3-D data consisting of several layers of different thermal conductivities following the order of each lithology (top to bottom) as presented in Table 1 (Fig. 9a, b). The modeled 3-D synthetic temperature distribution is compared to an analytical solution of 1-D heat conduction through heterogeneous (vertically layered) medium of Hasterok and Chapman (2011), the form of which is written as follows:  where T and q indicate the temperature and heat flow at a point of calculation, while A, , and z denote heat production, thermal conductivity, and thickness of the medium between two successive temperature and heat flow points. Subscripts i and i + 1 denote the position of each parameter relative to the other in terms of depth, such that, for example, T i+1 is the temperature calculated at a point below a depth for which the temperature T i was previously calculated (or given). For the thermal properties, the subscripts translate into the position of each medium relative to one another, i.e., whether it underlies or overlies the other unit. When used in a recursive manner, Eq. 4 becomes an analytical solution to the steady-state heat conduction equation for calculating temperature (and heat flow density) at any particular depths within a 1-D multi-layered medium.
The comparison between 1-D temperature profile extracted from the 3-D finite difference-modeled temperature distribution and vertical temperature calculated using Eq. 4 is displayed in Fig. 9c. It can be observed that the result derived using numerical solution matches that of analytical one. Therefore, our numerical conductive heat transfer modeling procedure is considered reliable and is subsequently employed to the onshore Northwest Java Basin.

Modeled temperature fields
The modeled 3-D subsurface temperature distribution is presented in the form of temperature-depth maps for each 1000 m depth interval (Fig. 10), as well as cross sections (Fig. 11). Here, we only present the original modeled temperature distribution, i.e., that generated using the mean thermal conductivity (Table 1) and a basal heat flow of 78.44 mW/m 2 ( Table 2). As for the results of the modeling using the maximum and minimum formation thermal conductivities and basal heat flow values, we use them only in the validation procedure ("Synthetic results" section) as well as for calculating the maximum and minimum HFD ("Modeled surface heat flow density" section). Figures 10 and 11 demonstrate that the conductive temperature shows an alternating pattern, i.e., higher subsurface temperatures are observed in sub-basinal areas, and conversely, lower subsurface temperatures can be seen to occur in basin high areas. This pattern is more explicitly exhibited in areas crossing the Tangerang High, Ciputat Subbasin, and Rengasdengklok High. The high temperature reaches > 300 °C beneath Ciputat Sub-basin and the low one is 180 °C beneath Tangerang and Rengasdengklok Highs at maximum depth. Specifically, except for the three aforementioned areas, the alternating pattern did not actually make a significant appearance until a depth of 3000-4000 m is reached.

Modeled surface heat flow density
The modeled 3-D temperature distribution allows the calculations of surface HFD values (Fig. 12) at all locations within the study area, with a resolution of the same size as that of lateral grid size used in the numerical modeling, i.e., 2000 × 2000 m 2 . Such calculations require only the near-surface temperature gradient of the model domain (here denoted by dT dz ) and the harmonically averaged thermal conductivity at each lateral grid node position.
While the former can be readily extracted from the 3-D temperature distribution, the latter is calculated using an equation of the form (Beardsmore and Cull 2001): where Kh i,j is the harmonic-averaged column thermal conductivity, d is the cell thickness (in this case the vertical grid spacing), and k is the thermal conductivity of the cell. Equation 5 states that the harmonic-averaged thermal conductivity at every lateral grid position (i, j) is a function of thermal conductivity values at all cells down to the maximum depth at that location as well as their thicknesses. The surface HFD at a node (q i,j ) can then be computed using:

Fig. 10
Temperature-at-depth maps for the entire onshore Northwest Java Basin region extracted from the 3-D conductive temperature distribution that is modeled using the mean formation thermal conductivities and a basal heat flow value of 78.44 mW m −2 . The geological structures are overlain on the maps In addition, the calculation of surface HFD values are also performed using the temperature distributions that have been modeled using the maximum and minimum thermal conductivity and basal heat flow values. These HFD values will be employed in the procedure explained in "Discussion" section.

Discussion
The phenomenon demonstrated by the temperature-at-depth maps ( Fig. 10; "Modeled temperature fields" section) can be explained by examining the characteristic of the subsurface geology, i.e., in terms of the depths where the lateral distribution of lithology starts to change significantly. Except for beneath the Tangerang High-Ciputat Subbasin-Rengasdengklok High regions, the lateral lithology within other areas (eastern Pasirputih Sub-basin, Pamanukan High, and Jatibarang Sub-basin) is more or less similar, i.e., they are comprised of sedimentary units at shallower depths. This occurs until Fig. 11 Cross sections extracted from the 3-D conductive temperature distribution that is modeled using the mean formation thermal conductivities and a basal heat flow value of 78.44 mW m −2 , running through profile lines A-A′ (a) and B-B′ (b) (see Fig. 8a for the location of profile lines). The lithology and their mean thermal conductivities are also superimposed on the temperature cross sections a uniform depth of around 3-4 km is reached beneath these regions. About and below this depth range, the lateral distribution of lithology is alternated between the basement and sediments.
The above suggests that the thermal conductivity contrast between the sedimentary basin filling and the crystalline basement is more significant than that between the sedimentary formations themselves, such that areas with thicker sedimentary piles have their isotherms bent upward, thus implying higher geothermal gradients. The reverse occurs for those of shallower basements (Fig. 11). This means that the conductive temperature distribution is sensitive to the basement's geometry.
There are at least two major features of Fig. 12 that deserve detailed discussions. The first one is that the distribution of modeled surface heat flow density values bears a resemblance with the overall lateral pattern of basement depth. Higher heat flow values are concentrated within areas of basement highs (e.g., Pamanukan and Rengasdengklok Highs) and lower ones are clustered within sub-basinal areas (e.g., Ciputat Sub-basin). This is a direct manifestation of the control of thermal conductivity on the thermal regime, i.e., high heat flow is attributed to high thermal conductivity and, therefore, to locations where thermally conductive rocks (i.e., the basement) are dominant, and vice versa.
A more careful observation leads to the discovery of another HFD pattern, i.e., there are areas where heat flow density values are lower or higher but are not directly positioned above sub-basins or basin highs, but rather, their edges. Two most noticeable examples to this are the low heat flow in the central part of Tangerang High and an exceptionally high heat flow in the southern part of the Rengasdengklok High (between Ciputat and Pasirputih Sub-basins). This demonstrates that the effect of heat refraction (Beardsmore 2004;Thakur et al. 2012;Rawling et al. 2013; "General description of heat flow density distribution" section), i.e., where heat flows more preferably along thermally conductive zones rather than it does along insulating ones, prevails on these regions. This explains why the computed HFD in the southern Rengasdengklok High is Fig. 12 Map showing the predicted surface HFD distribution of the onshore Northwest Java Basin, computed from the three-dimensional conductive temperature distribution that is modeled using the mean formation thermal conductivities and a basal heat flow value of 78.44 mW m −2 . The wells, basement depth contours, and geological structures are superimposed on the map high compared to that of other areas within the basin. There is a short lateral transition between a basin low region (Ciputat Sub-basin), a basin high area (Rengasdengklok), to another basin low (Pasirputih Sub-basin), allowing heat from the low areas to be strongly imparted to the shallow but narrow basement beneath this region. In a reversed manner, the very low modeled HFD in central Tangerang High is probably due to it being situated between two narrow basement highs.
As a final attempt to analyze and interpret the HFD distribution of the onshore Northwest Java Basin, we perform comparisons between the observed HFD with their associated uncertainty (Fig. 6a) and the calculated HFD ( Fig. 12) with their predicted uncertainty. A similar procedure has been attempted by Cacace et al. (2013) for an area in northern Germany. In principle, any deviations of the observed HFD values from the modeled ones that are beyond the range of uncertainty of the latter should be mostly attributed to the redistribution of conductive heat by groundwater flow.
We subsequently define the concept employed for the interpretation, i.e., the HFD of a particular well is considered as being anomalous when: where q is the HFD, and the subscripts max and min correspond to the maximum and minimum observed or modeled HFD. For the observed HFD, the maximum and minimum HFD values are equal to the well HFD values plus and minus their uncertainty (Fig. 6a), respectively. The maximum and minimum predicted HFD can be directly obtained by applying Eq. 6 to the 3-D temperature distributions generated from the modeling runs using maximum and minimum values of thermal conductivity and basal heat flow ("Model domain setup" section).
The result of this procedure is shown in Fig. 13, in the form of point maps of the differences between the maximum and minimum modeled and observed HFD for each well used in the present study. In addition, we also provide a summary univariate statistics (Table 3) and construct histograms based on the resulting HFD differences (Fig. 14). The mean and median values of HFD differences computed using Eq. 7b are relatively similar to each other, in contrast to that of the values calculated using Eq. 7a. This is also emphasized by comparing their standard deviation-to-mean ratios, which indicates that values computed from Eq. 7a possess more variability. The coefficient of skewness is positive for values calculated using Eq. 7b, and lower than zero for those computed using Eq. 7a. These are also reflected in the shape of each histogram, once again owing to the presence of few extreme values in each. This phenomenon, in which few highly negative HFD difference values (lower than the mean minus one standard deviation) appear in the histogram of Fig. 14a, is caused by the presence of exceedingly high observed HFD values at several localities, which implies that the difference between the modeled and observed maximum HFD values there is enormous. Therefore, it follows that these high observed HFD values also cause the histogram of Fig. 14b to display extremely high values, i.e., the minimum modeled HFD values are much smaller than the minimum observed HFD in the said locations. As a result, both histograms are skewed in opposite directions. Taking into account the aforementioned difference in variability between values computed using Eq. 7a and 7b, this means that there are more observed HFD values exceeding their modeled counterparts than there are observed HFD values lower than Fig. 13 Map showing the distribution of the difference between: a the maximum observed and predicted HFD and b the minimum observed and predicted HFD. The groundwater flow regime of a number of areas within the onshore Northwest Java Basin-based hydrogeological studies of Kloosterman (1989) and Lubis et al. (2008). For explanation see text their modeled ones. This is proven by the last row of Table 3, in which there are more numbers of data points whose HFD differences are less than 0 mW m −2 when calculated using Eq. 7a than there are resulting from Eq. 7b. The anomalously high and low observed HFD values (i.e., those yielding negative values when inserted into calculations using Eq. 7a and 7b, respectively) are subsequently interpreted in terms of mode of heat transfer, while the extremely high ones are attributed to a geological cause. These interpretations are addressed in more detail through the following discussions. In relevance to the interpretation concerning the effect of groundwater flow, we provide the distribution of shallow groundwater flow zones across the basin, which were originally constructed by Kloosterman (1989) for the areas adjacent to and surrounding Fig. 14 Histograms and cumulative frequency distributions of a the differences between calculated and measured HFD maxima (values computed using Eq. 7a), and b the differences between calculated and measured HFD minima (values computed using Eq. 7b). See also Table 3 for the univariate statistics Rengasdengklok High, Pasirputih Sub-basin, Pamanukan High, and Jatibarang Subbasin, and by Lubis et al. (2008) for the Ciputat Sub-basin (the capital city of Jakarta region). The differences between the two HFD maxima (Eq. 7a; Fig. 13a) are in agreement with the groundwater regime. The wells RDH-2 (northern Rengasdengklok High), CLU-5 (northern Pasirputih Sub-basin), PRD-1, GTR-1/2/3, PCT-1 (western Jatibarang Sub-basin) all are anomalous according to Eq. 7a, i.e., the observed HFD plus their uncertainties are higher than the maximum predicted HFD. This is due to them being located in or adjacent to regional groundwater discharge areas, such that the observed HFD is higher than what would be expected from a purely conductive or "background" HFD (e.g., Kilty and Chapman 1980;Cacace et al. 2013). Other wells with anomalous HFD, i.e., SUM244 (Tangerang High), JNG-1 (Ciputat Sub-basin), CKR-1 (Rengasdengklok High), JRR-1/2, and PJB-1 PJN-P1 (western and southern Pasirputih Sub-basin) are not covered by the groundwater regime information. For well SUM244, this might be due to heat refraction (due to it being situated adjacent to the southern part of Rengasdengklok High; see Fig. 11) that is higher in magnitude than that can be accommodated by the maximum modeled HFD, given the extremely rugged basement. Well CKR-1 is adjacent to a relatively large basement fault, thus the influence of deep fluid circulation is likely. For the other wells, such as those located in the southeastern edge of the Pasirputih Sub-basin, magmatic activities may fulfill the necessary condition for generating the tremendously high observed HFD, i.e., up to 171.70 and 219 mW m −2 (Fig. 6a), well above a typical conductive HFD value right over active volcanic areas of Southeast Asia, i.e., about 140 mW m −2 (Nagao and Uyeda 1995). These wells are indeed relatively close to the volcanic arc south of a major thrust fault and are not within a groundwater recharge area. Well JNG-1 is neither located in an area where significant heat refraction occurs nor is it within a groundwater discharge area. In addition, it is not located in close proximity to any nearby basement fault that may facilitate deep groundwater circulation and advective heat transport, so it remains uncertain of what effect may be responsible for the anomalously high HFD. A possible explanation would be that the near-surface fault located just to the west of JNG-1 cuts through the confining layer in the area, thereby providing a path for hotter, deeper, groundwater to ascend through to a layer containing shallower groundwater body.
At the other end of the spectrum, we can observe in Fig. 13b that there are also wells whose minimum observed HFD values are actually lower than their modeled minima. Although, these wells are mostly situated in groundwater recharge or mixed rechargedischarge zones, such that their HFD is lowered due to the absorption of heat by fluid flow (Kilty and Chapman 1980), there are also wells which show similarly lower-thanpredicted HFD, but are not directly within recharge areas, or even adjacent to discharge areas. Examples are WNJ-1, BJR-1, CJT-1, and MLD-1 (a group of wells on the northern part of Pamanukan High). Recalling Fig. 12, this could again be due to the effect of heat refraction whose magnitude is even higher than that can be accounted by the modeled HFD, i.e., conductive heat flow beneath the area is diverted to the nearby crest of Pamanukan High. It is also possible that a more dominant advective-convective heat transfer occurs in the form of groundwater recharge in deeper sections of this area, as evidenced by the presence of faults in the proximity of these wells.