Estimation of the Impact of Basement Heterogeneity on Thermal History Reconstruction: The Western Siberian Basin

: Some of the simplifying assumptions frequently used in basin modelling may adversely impact the quality of the constructed models. One such common assumption consists of using a laterally homogeneous crustal basement, despite the fact that lateral variations in its properties may signiﬁcantly affect the thermal evolution of the model. We propose a new method for the express evaluation of the impact of the basement’s heterogeneity on thermal history reconstruction and on the assessment of maturity of the source rock. The proposed method is based on reduced-rank inversion, aimed at a simultaneous reconstruction of the petrophysical properties of the heterogeneous basement and of its geometry. The method uses structural information taken from geological maps of the basement and gravity anomaly data. We applied our method to a data collection from Western Siberia and carried out a two-dimensional reconstruction of the evolution of the basin and of the lithosphere. We performed a sensitivity analysis of the reconstructed basin model to assess the effect of uncertainties in the basement’s density and its thermal conductivity for the model’s predictions. The proposed method can be used as an express evaluation tool to assess the necessity and relevance of laterally heterogeneous parametrisations prior to a costly three-dimensional full-rank basin modelling. The method is generally applicable to extensional basins except for salt tectonic provinces.


Introduction
Temperature is one of the essential parameters for determining the rate and rank of organic matter maturation [1,2]. The reliability of temporal temperature estimations is determined, primarily, by the reliability of the thermal properties used [3], the boundary conditions and the thermal modelling approach [4][5][6]. A reliable thermal model of a basin should consider geological processes of equal importance in both the lithosphere and the sedimentary cover and their effects on each other [1,7].
Basin modellers use a number of simplifying assumptions when reconstructing sedimentary basin evolution through numerical simulation. One common assumption in most basin models is a laterally homogeneous basement. This assumption is usually connected with a scarcity of data or with its poor quality near the bottom of the basin and/or of the crustal basement. This occurs because the resolution of surveys deteriorates with depth due to the natural limitations of geophysical observations, as well as due to economic constraints on costly high-resolution surveys. In some cases this may lead to the introduction of "acoustic" [8] or "economic" [9] basements in the model. Thus, due to a lack of data, basin modellers oftentimes introduce a laterally homogeneous crystalline basement instead of a poorly explored bottom part of the sedimentary lithospheric strata. Such an assumption is actually equivalent to populating the model with incorrect thermophysical properties of the basement, hoping that the effect of this rough approximation on the thermal history could be compensated for by the calibration of the thermal model using wellbore data. However, successful wellbore calibration does not guarantee the exactness of the thermal history reconstruction, especially when it comes to inter-well segments [3].
The heterogeneity of the basement in the basin models is rather scarcely covered in the literature [4,[10][11][12][13][14][15][16]. Basement heterogeneity has only been considered in a few works with respect to its effect on the basin's thermal evolution. Klitzke et al. [10,11], Fjeldskaar et al. [12] and Scheck-Wenderoth et al. [13,14] used seismic imaging to determine the geometry of the basement's heterogeneous blocks, although the time-depth conversion may involve significant uncertainties at such depths. They subsequently used a gravity data analysis to determine the density of rocks. Fattah et al. [15] varied the crust's thickness and composition to calibrate the model with regard to gravity data. However, the impact of the basement's heterogeneity on the heat flow and on the maturity of source rocks was not quantified. In all the above mentioned works [10][11][12][13][14][15], decoupled thermal solutions were used (a comparative analysis of coupled vs. decoupled solutions can be found in [4,5]). To the best of our knowledge, only Hansford [16] and Clark et al. [4] have used a coupled structural and thermal solution to determine the presence of magmatic bodies in the basement and their impact on thermal history. A sensitivity study of the reconstructed thermal properties has not yet been attempted, despite a high variability of the thermal properties and the densities of basement rocks [17,18].
Based on the studied literature, we may conclude that: (1) the heterogeneity of density and of the thermal properties as well as lateral variations of the basement's geometry are usually ignored in the thermal modelling of basins and of petroleum systems, and (2) the impact of basement's heterogeneity on the basin thermal history is scarcely covered in the published studies.
One of the major reasons behind both these observations is the complexity of the inversion procedure aimed at the reconstruction of the basement's properties and of its geometry. The classical inversion procedure is time-consuming and typically requires data from a variety of geophysical studies as well as contributions by experts from a number of fields [19]. At the same time, this classical inversion procedure possesses a certain weakness in terms of basin modelling, whereby it lacks the possibility of a direct reconstruction of the thermal properties, which is crucial for the reconstruction of the thermal history. Thermal properties can only be extracted from subsequent non-unique and highly uncertain interpretations of rock density (see, e.g., in [13,14]). Of course, when wellbore material is available, reliable thermal properties may be obtained through direct measurements on core samples or at least through the interpretation of lithology. A basement is not always drilled, and in the cases in which it is, the wellbore coverage is frequently sparse.
In this study, we propose an express estimation method, which establishes the need for considering the heterogeneous basement for thermal history reconstruction in fullrank three-dimensional (3D) basin models. The express estimation method is based on reduced-rank inversion using the capabilities of two-dimensional (2D) basin modelling simulators. In contrast to the customary inversion, the reduced-rank inversion does not use seismic data as a pre-requisite and allows us to obtain the petrophysical properties of the heterogeneous basement and its geometry. The practical value of the express estimation method is demonstrated in the example of an area in the Western Siberian Basin.

Method and Data
The proposed method aims to estimate the impact of basement heterogeneity on thermal models and source-rock maturity models. The method includes a comparative analysis of a set of basin models with a heterogeneous basement with a reference model built with a laterally homogeneous basement. The construction of each 2D basin model follows the workflow presented in Scheme 1. The workflow, overall, repeats the standard procedure (e.g., [20,21]) and is extended by a newly proposed reduced-rank inversion (Scheme 1, step 3; Appendix A) for basement heterogeneity reconstruction. The applicability of the presented method is limited by extensional types of basins without salt tectonics.
The method assumes that the construction of a set of models includes the following: (1) a model with laterally homogeneous properties of the basement as a reference model (Model 1), (2) a model with directly determined basement heterogeneous properties (Model 2) based on the workflow from Scheme 1, (3) two models with minimum (Model 2: ρ min ) and maximum (Model 2: ρ max ) possible density values to evaluate the uncertainties in the density effect on thermal history reconstruction (following Scheme 1) and (4) a series of models for a sensitivity analysis to assess the impact of uncertainties in thermal conductivity in the heterogeneous basement on the predictions of thermal history and of the petroleum system maturity for all models. Considering the uncertainty evaluations in the heterogeneous basement simultaneously for both density and thermal conductivity allows us to obtain the ultimate possible associated range of the heat flow solutions.
There are four steps that comprise the workflow. In step 1, a classic basin modelling workflow described, e.g., by Hantschel and Kauerauf [21], was used to create a sedimentary cover model. The model was calibrated on measured porosity (ϕ) and pore pressure (P). Finally, a thermally non-calibrated geological model of a basin was obtained. the workflow presented in Scheme 1. The workflow, overall, repeats the standard procedure (e.g., [20,21]) and is extended by a newly proposed reduced-rank inversion (Scheme 1, step 3; Appendix A) for basement heterogeneity reconstruction. The applicability of the presented method is limited by extensional types of basins without salt tectonics. The method assumes that the construction of a set of models includes the following: (1) a model with laterally homogeneous properties of the basement as a reference model (Model 1), (2) a model with directly determined basement heterogeneous properties (Model 2) based on the workflow from Scheme 1, (3) two models with minimum (Model 2: ρmin) and maximum (Model 2: ρmax) possible density values to evaluate the uncertainties in the density effect on thermal history reconstruction (following Scheme 1) and (4) a series of models for a sensitivity analysis to assess the impact of uncertainties in thermal conductivity in the heterogeneous basement on the predictions of thermal history and of the petroleum system maturity for all models. Considering the uncertainty evaluations in the heterogeneous basement simultaneously for both density and thermal conductivity allows us to obtain the ultimate possible associated range of the heat flow solutions.
There are four steps that comprise the workflow. In step 1, a classic basin modelling workflow described, e.g., by Hantschel and Kauerauf [21], was used to create a sedimentary cover model. The model was calibrated on measured porosity (φ) and pore pressure (P). Finally, a thermally non-calibrated geological model of a basin was obtained. Scheme 1. General workflow of the thermal history reconstruction of a riftogenic basin with a heterogeneous basement, where φ is porosity, P is pressure, T is temperature, Ro is vitrinite reflectance, PWD stands for paleo-water depth, MLC stands for Moho-level calibration, GC stands for gravity calibration and TC stands for thermal calibration. The solid arrows show the path for the coupled solution used in the case study, the dashed arrows show the path for a backstripping-based solution, and the bidirectional arrows indicate additional simulations during model calibration.
In step 2, homogeneous basement layers were introduced into the model. Estimations of the present-day Moho depth were used to determine the crust's initial (i.e., before the rifting extension) thicknesses. The thermal calibration of the model, performed in step 4, when avoiding step 3, provided a reference model without basement heterogeneity (Model 1); this model should be used in step 3 as a precursor for constructing models with heterogeneous basements (Model 2, Model 2: ρmin, Model 2: ρmax).
After these steps, the modelling workflow bifurcated according to the selected modelling approach, which may involve either coupled or decoupled thermal and structural Scheme 1. General workflow of the thermal history reconstruction of a riftogenic basin with a heterogeneous basement, where ϕ is porosity, P is pressure, T is temperature, Ro is vitrinite reflectance, PWD stands for paleo-water depth, MLC stands for Moho-level calibration, GC stands for gravity calibration and TC stands for thermal calibration. The solid arrows show the path for the coupled solution used in the case study, the dashed arrows show the path for a backstripping-based solution, and the bidirectional arrows indicate additional simulations during model calibration.
In step 2, homogeneous basement layers were introduced into the model. Estimations of the present-day Moho depth were used to determine the crust's initial (i.e., before the rifting extension) thicknesses. The thermal calibration of the model, performed in step 4, when avoiding step 3, provided a reference model without basement heterogeneity (Model 1); this model should be used in step 3 as a precursor for constructing models with heterogeneous basements (Model 2, Model 2: ρ min , Model 2: ρ max ).
After these steps, the modelling workflow bifurcated according to the selected modelling approach, which may involve either coupled or decoupled thermal and structural solutions. We recommend the coupled thermal and structural approach [7] since it considers sediment stretching, coupled temperature and isostasy modelling with subsidence and does not require information on the paleo-water depth, in contrast to the backstripping-based approach with forward temperature modelling [4,5].
In step 3, reconstruction of the heterogeneous basement following the reduced-rank inversion (see a detailed, step-by-step description in Appendix A) started. The inversion aims to reconstruct the heterogeneous basement properties and geometry (width and thickness of each block). The model with the directly determined properties and geometry of the heterogeneous basement was built (Model 2). Here, the following assumptions were made:

•
The upper part of the basement is heterogeneous. • Under the heterogeneous basement is a crystalline basement.

•
The densities of the lithospheric layers, upper (ρ u_crust ) and lower (ρ l_crust ) crusts and the lithospheric mantle (ρ l_mantle ) are constant and change with temperature.

•
The present-day width of each rectangular block is determined based on the basement maps. • The first approximation of the thickness of each heterogeneous basement block is set.

•
The physical properties of each heterogeneous block are based on the lithological descriptions of maps and their interpretations by standard mixing rules [21].
The physical properties of the rock for lithological mixtures, such as density (ρ), specific heat capacity (C ρ ), radiogenic heat production (A) and thermal expansion (α), were calculated using the arithmetic mean rule, while thermal conductivity (λ) was calculated using the geometric mean rule. Measurements on rock samples from outcrops and data from well logs and core samples may significantly narrow the ranges of inferred petrophysical properties. Whenever possible, the techniques for thermal conductivity determination should be used (see, e.g., [22,23]).
The final geometry and density model of the heterogeneous basement was obtained through the iterative simulation technique, e.g., as in Bott [24], matching the observed gravity anomalies to the modelled ones. The densities of the sedimentary basin and the crystalline basement were fixed. The thickness of a heterogeneous block was inversed, while its density was fixed. The latter should be corrected in the rare cases in which the thickness of a block falls out of the range between the minimum determined block thickness and a thickness greater than the crust thickness beneath the block. The modelled anomaly may be recalculated, e.g., as in Heiland [25], only after correcting all block thicknesses and the model simulation.
At the beginning of thermal calibration (Scheme 1, step 4), the value of radiogenic heat production (A) in the crustal layer was fitted to the wellbore values of temperature (T) and vitrinite reflectance (%Ro) at the regional scale [26,27]. If the thermal calibration does not produce a good fit in some wells, which cannot be explained by uncertainty in thermal conductivity, then the basement model needs to be refined within the possible ranges of density and the geometry of the blocks. Thus, the gravity and thermal calibrations complement each other and constrain the diversity of solutions.
Furthermore, following Scheme 1, two additional models with minimum (Model 2: ρ min ) and maximum (Model 2: ρ max ) heterogeneous basement densities were constructed. In contrast to Model 2, these models had a different last assumption regarding the density of the heterogeneous basement. Here, the minimum and maximum values for the determined rock lithology were based on the database (e.g., [28]).
A sensitivity analysis was performed to estimate the impact of the uncertainties in the thermal conductivity of the heterogeneous basement on thermal history reconstruction with Model 2, Model 2: ρ min and Model 2: ρ max . We considered systematic changes to matrix thermal conductivity (λ) in the ±20% and ±50% range according to [29] and [3], respectively.
When all the models were built and calibrated and once a sensitivity analysis was performed, the significance of the heterogeneous basement's impact on thermal history reconstruction was determined. Following the analysis, whether the heterogeneous basement is needed in constructing a full-rank 3D basin model can be determined, and, if required, additional geophysical studies (e.g., deep seismic sounding, density and thermal property measurements) should be planned.

Central Part of the Western Siberian Basin
The Western Siberian Basin is of great interest for evaluating the presented express estimation methodology's efficiency from two perspectives. Firstly, this basin is one of the largest petroleum provinces in the world [30]. Secondly, most of the basin modellers prefer to assume the pre-Jurassic strata as a crustal laterally homogeneous layer (e.g., [31][32][33][34][35]), although these strata actually have a complex heterogeneous nature [36]. The reason for such an assumption is twofold. Firstly, the seismic data along the seismic study profile do not provide reliable reflections to reconstruct the Palaeozoic sedimentation. Secondly, according to an appraisal of the existing seismic data of entire Siberia by Cherepanova et al. [37], the rifts in our study area cannot be distinguished in the seismic models.
In this work, we considered these Paleozoic strata and syn-rift regional-scale intrusion in rift-graben zones as a heterogeneous basement formed instantaneously before the modelled Mesozoic-Cenozoic sedimentation. In the literature, the heterogeneous basement is referred to as a pre-Jurassic folded basement [38,39].

Geological Framework and Data Set
The geological evolution of Western Siberia has been described in detail by several authors [30,[39][40][41] and is briefly summarised here. This study used regional sub-meridional seismic profile 6 [39], which crosses the centre of the Western Siberian Basin (Figure 1a).

Heterogeneous Basement
The Western Siberian intra-cratonic basin has experienced three phases of evolution with regard to tectonics [30,41]. The first phase is referred to as the Riphean-Vendian and Paleozoic ages. The corresponding basal structural stage was formed during Hercynian compression and was partly eroded by the end of the Permian period. The basin is composed of a peneplenised folded and metamorphosed orogenic system and platformal blocks [39]. In the second phase, the basin extended during the Permo-Triassic mega-rifting. The voluminous intercontinental flood basalts were extruded in the rift system, covering the basin in the northward direction [42,43]. Both of these stages are reflected in the so-called heterogeneous basement structure.
The lithological description of the heterogeneous basement (Table 1) was transferred from the maps [44][45][46] to bulk material properties for each geological unit following the standard procedure [21]. Density (ρ), thermal conductivity (λ), specific heat capacity (C ρ ), radiogenic heat production (A) and the thermal expansion of each block in Model 2 were defined based on the petrophysical database and the application of the mixing rules [21]. The Sekiguchi model [47] is used for the temperature correction of thermal conductivities (λ). Data on the volumetric thermal expansion coefficients (α) are in good agreement with the published values [17,18,48].
Thermal conductivity varies in accordance with lithology (at least composition and porosity), alteration state, in situ pressure, temperature, and saturation (see, e.g., [49][50][51]), and published rock thermal conductivity data often fail to provide actual information (see, e.g., [3]). For example, [49] demonstrates the variability in the bulk thermal conductivity values of up to 60% within the volcanic rock sequences with a thickness of less than 1 km in one well. The bulk thermal conductivity decreases as the chemical index of alteration and porosity value increase [50]. For example, the database value of basalt rocks' matrix thermal conductivity (  [55] and 2.77 W/m/K for the serpentinite sample from the Nizhne-Tabachanskaya area of the Western Siberian Palaeozoic basement studied by Duchkov et al. [54]. In the absence of reliable data, we suggested quantifying the impact of uncertainty in the thermal conductivity of the heterogeneous basement on thermal history reconstruction by performing a sensitivity analysis varying λ of the heterogeneous basement blocks by ±20% and ±50% relative to the values determined from the database.  [39]); (b) present-day simplified stratigraphy of the regional profile 6 linked with the lateral limits of the heterogeneous basement blocks in Model 2 (modified after [6]). The number of blocks and colours correspond to the material properties presented in Table 1. The names of petroleum regions are presented in italics. (c) The geological map of the heterogeneous basement (modified after [44][45][46]); the geological description used in the study blocks is presented in Table 1. The green curve highlights the profile location.  [39]); (b) present-day simplified stratigraphy of the regional profile 6 linked with the lateral limits of the heterogeneous basement blocks in Model 2 (modified after [6]). The number of blocks and colours correspond to the material properties presented in Table 1. The names of petroleum regions are presented in italics. (c) The geological map of the heterogeneous basement (modified after [44][45][46]); the geological description used in the study blocks is presented in Table 1. The green curve highlights the profile location.  [44][45][46]. A Here, ρ is the density; its acceptable range according to [49] is shown in parentheses; α is the volumetric coefficient of thermal expansion; for unit numbers 1-21, A is the radiogenic heat production of the rock matrix; for the crust and lithospheric mantel, A is surface radiogenic heat production (A 0 ) in the equation A(z) = A 0 × exp(z/a r ), where A(z) is bulk radiogenic heat production, z is the negative value of a depth (km) and a r is the e-fold length parameter (km) [27]. C ρ is the rock matrix specific heat capacity, and λ is the rock matrix thermal conductivity at 20 • C. Values are calculated based on the petrophysical database in [21], with the geometric weighted mean rule for λ and the arithmetic weighted mean rule for ρ, C ρ and A. B Questionable values (see the text for details). * Manually truncated edge of density range by [49] to obtain the gravity anomaly calibration for Model 2: ρ min and Model 2: ρ max .
For the same rock type, density values also can exhibit significant variations, see [49] for an example for volcanic rocks and the corresponding ranges for different types of rocks [49]. The variations in density considered in this work were taken from Schön [49] with modifications.

Rift Phases
The Western Siberia Basin is characterised by two rift phases [56]. The first is Permo-Triassic mega-rifting. The second is Neocomian, which is localised in some areas of the basin [56]. According to tectonic analysis in [56], the area of the basin understudy experienced only the Triassic phase. Thus, the rift period was defined from 251 to 228 Ma, with differential thinning of the crust and the mantle [39].

Sedimentary Basin: Stratigraphy, Infill, Erosion
The platformal stage of basin evolution started in the Mesozoic era and was accompanied by uplift and erosion in the Oligocene-Pliocene period. The erosion of sediments with a thickness of 20 to 400 m was set from 23 to 10 Ma [57]. The dataset used in this work was mainly based on the basin model presented in [6]. The section of input stratigraphy is shown in Figure 1b using colour-coded age. Each stratigraphic unit has distinct material properties, as in [6]. The corresponding data are also presented in Table A1.
When modelling the basin evolution, the bulk thermal conductivity of rocks was calculated using the geometric mean rule [58,59]. Athy's exponential low [60] was used for porosity calculations. The Sekiguchi model [47] was applied to consider the temperature dependence of matrix thermal conductivities. The Waples formula was used for a temperature correction of the specific heat capacity [61].

Petroleum Systems
Oil and gas accumulations have been revealed in eight oil and gas plays (see details in [30,40,[62][63][64][65]). The deepest buried play is present at the Paleozoic paleobasin at the top of the heterogeneous basement. The Triassic play is present in the north of the basin and in rift zones. The Togur/Tymen, Vasyugan and Bazhenov-Abalak plays are proven to be in the Lower-Middle Jurassic strata. The Neocomian, Aptian-Albian-Cenomanian and Senonian are plays of the Cretaceous strata. The primary source rock is the Late Jurassic Bazhenov Formation (Fm) (Figure 1b), which plays the most significant role in oil generation [64,65].

Boundary Conditions
At the upper boundary, the temperature decrease was defined in time within the 20-0 • C range, following the method of Isaev [66]. The lower boundary condition for temperature was set at the lithosphere-asthenosphere interface at a constant value of 1300 • C [67].

Gravity Anomaly Data
The reference gravity anomaly values were digitised along the investigated section from a map of the Bouguer anomaly of the Federal Ural district (the map is accessible by a web link in [68]). The map was built with the average topographic density of 2670 kg/m 3 , with the terrain correction radius r = 200 km, the normal gravity value derived from the Helmert formula (1901-1909) and a datum correction of −14 mGal to the International Gravity Standardisation Net [69]. The gravity survey was performed on a regional scale of 1:7,500,000, with a counter step of ε c = 5 mGal. Other details of the gravity anomaly map building are not known. The uncertainty value was not provided; therefore, ε t = ±1.5 mGal that corresponds to a typical requirement for the regional-scale Bouguer anomaly map was assumed [70] (see details in Appendix C).
In this work, the desired accuracy (see details in Appendix A) to fit the model and the measured results of the Bouguer anomaly, of σ = 6.5 mGal, was defined as the sum of the counter step (ε c ) and the uncertainty value of measurements (ε t ) (see details in Appendix C).

Construction and Calibration of Basin Models
All of the built basin models had the same configuration: a sedimentary cover, temperature boundary conditions and rifting setting. The models were calibrated by pore and pressure regimes, Moho depth and a wellbore thermal regime. Referring to Scheme 1, Model 1 resulted from the application of steps 1, 2 and 4. Thus, Model 1 was not fitted with the observed gravity anomaly. Model 2 was built with the heterogeneous basement, following all steps in Scheme 1 and the newly proposed inversion procedure (Scheme A1). For our case study, models labelled Model 2: ρ min and Model 2: ρ max were built with the minimum and maximum density values of the basement's igneous rocks. Based on our modelling results, heterogeneous blocks with a sedimentary genesis had a negligible impact on thermal history reconstruction. Several additional models, which considered the impact of uncertainty in the thermal conductivity (λ) of the heterogeneous basement on thermal history reconstruction, were referred to as Model 2: λ ± 20%; Model 2: λ ± 50%; Model 2: ρ min , λ ± 50% and Model 2: ρ max , λ ± 50%. The matrix thermal conductivity of the crust and the lithospheric mantle remains unchanged.
The initial thicknesses at 251 Ma of the upper crust, the lower crust and the lithospheric mantle were set to H UC0 = 30 km, H LC0 = 10 km and H UM0 = 70 km, as in [6] (Figure 2a).

Calibration of Model 1
After construction, Model 1 was run to check the calibration of the stratigraphy, porosity and pressure regimes, Moho depth and gravity anomaly: • A good match between the modelled stratigraphy and the input stratigraphy was obtained with an approximate convergence misfit of 5%. The result was reached after 15 inversion iterations of the forward modelling [7].

•
The calculated present-day Moho depth is in good agreement with two published interpretations [37,42] within the interval of 190 km to the eastern end of the profile (Figure 3, red line). The absence of the flexural load by the Ural fold belt on the west in the model explains the discrepancy between the calculated and published data within the left interval of 0 to 190 km.

•
The gravity anomaly ( Figure 4, red line) was not fitted since the modelled σ m = 10.5 mGal was greater than the desired accuracy of σ = 6.5 mGal.
The calibration steps outlined above provided only a robust quality check of the structural and sedimentary evolution models; the thermal solution needed additional calibration using the wellbore data. The dataset considered included present-day temperatures in the boreholes at profile locations of 283, 406, 572 and 700 km (Figures 1b and 5a) and vitrinite reflectance data in the profile points at 163, 190, 278, 283, 420 and 558 km (Figures 1b and 5b) from [6]. Thus, only the borehole at 283 km provided both temperature and vitrinite reflectance simultaneously. The thermal solution was only calibrated on a regional scale, i.e., without local tuning. An exponential reduction of the radiogenic heat production of the crust (A) from a surface value (A 0 ) at a rate given by the e-fold length (A r ) with depth (z) was assumed [27] (see the footnote in Table 1). The e-fold length parameter (A r ) selected for Model 1 was equal to 20.4 km. The well thermal data had a reasonably good match with the modelled thermal regime ( Figure 5). Hence, the thermal properties did not require optimisation in this case.
All the modelled temperature values fell within the ±7% interval, relative to the observations (Figure 5a), which is considered acceptable (see details in Appendix C). A simultaneous calibration of temperature and vitrinite reflectance (%Ro) was carried out for three kinetic models, referred to as EasyRo [73], BasinRo [74] and EasyRoDL [75]. Only the first of them provided a successful calibration of the model. Most vitrinite reflectance data fell within ±0.05 %Ro (Figure 5b), which is also acceptable (see details in Appendix C). The only value of vitrinite reflectance measured for the wellbore at 190 km exceeded the model's permissible confidence interval. Two measured %Ro values for this wellbore (0.45% and 0.66%) were significantly different, although both were obtained from almost the same depth (depth difference~100 m). This can be explained by either the random error of measurements or the redeposition of organic matter.  Coloured blocks correspond to the heterogeneous basement with assigned petrophysical properties from Table 1. The layer of the lithospheric mantle is not displayed for the convenience of presenting the picture.
( Figures 1b and 5b) from [6]. Thus, only the borehole at 283 km provided both temperature and vitrinite reflectance simultaneously. The thermal solution was only calibrated on a regional scale, i.e., without local tuning. An exponential reduction of the radiogenic heat production of the crust (A) from a surface value (A0) at a rate given by the e-fold length (Ar) with depth (z) was assumed [27] (see the footnote in Table 1). The e-fold length parameter (Ar) selected for Model 1 was equal to 20.4 km. The well thermal data had a reasonably good match with the modelled thermal regime ( Figure 5). Hence, the thermal properties did not require optimisation in this case.   All the modelled temperature values fell within the ±7% interval, relative to the observations (Figure 5a), which is considered acceptable (see details in Appendix C). A simultaneous calibration of temperature and vitrinite reflectance (%Ro) was carried out for three kinetic models, referred to as EasyRo [73], BasinRo [74] and EasyRoDL [75]. Only the first of them provided a successful calibration of the model. Most vitrinite reflectance data fell within ±0.05 %Ro (Figure 5b), which is also acceptable (see details in Appendix C). The only value of vitrinite reflectance measured for the wellbore at 190 km exceeded the model's permissible confidence interval. Two measured %Ro values for this wellbore (0.45% and 0.66%) were significantly different, although both were obtained from almost the same depth (depth difference ~100 m). This can be explained by either the random error of measurements or the redeposition of organic matter.  All the modelled temperature values fell within the ±7% interval, relative to the observations (Figure 5a), which is considered acceptable (see details in Appendix C). A simultaneous calibration of temperature and vitrinite reflectance (%Ro) was carried out for three kinetic models, referred to as EasyRo [73], BasinRo [74] and EasyRoDL [75]. Only the first of them provided a successful calibration of the model. Most vitrinite reflectance data fell within ±0.05 %Ro (Figure 5b), which is also acceptable (see details in Appendix C). The only value of vitrinite reflectance measured for the wellbore at 190 km exceeded the model's permissible confidence interval. Two measured %Ro values for this wellbore (0.45% and 0.66%) were significantly different, although both were obtained from almost the same depth (depth difference ~100 m). This can be explained by either the random error of measurements or the redeposition of organic matter.
(a) (b) Figure 5. Calibration results of present-day measured and modelled (a) temperatures and (b) vitrinite reflectance by the EasyRo kinetic model [73]. Italic numbers near the data points correspond to the wellbore location (in km) along the profile. Grey zones show the area of ±7% for temperature and ±0.05% (absolute) for vitrinite reflectance.
The final geological section is presented in Model 1 in in Figure 6a.  [73]. Italic numbers near the data points correspond to the wellbore location (in km) along the profile. Grey zones show the area of ±7% for temperature and ±0.05% (absolute) for vitrinite reflectance.  Table 1. The corresponding modelled gravity anomalies are presented in Figure 4. The layer of the lithospheric mantle is not displayed for the convenience of presenting the picture.  (Figure 1b). The full version of the maps are accessible via the web links in [44][45][46]. The material properties from Table 1 were assigned to each block, and the thicknesses of the heterogeneous basement blocks ( Figure  2c) were reconstructed with the help of gravity anomaly fitting (Scheme A1). For Model 2, the density values of units 8 and 18 were manually adjusted from 2870 and 2740 to 2840 and 2840 kg/m 3 , respectively, (Table 1) according to Schön [49] (p. 116).

Construction of
The uncertainty of the determined density of the blocks in the basement may affect the thickness of the reconstructed blocks, impacting the entire thermal regime reconstruction. We considered the uncertainty of the density only for igneous rocks since they have a critical impact on thermal history. We built Model 2: ρmin and Model 2: ρmax using the range of acceptable densities for igneous rocks taken from Schön [49] (p. 116) for units 1, 2, 3, 4, 7, 8, 15, 16, 18 and 19. The basement geometry was obtained following the reducedrank procedure (Scheme A1), varying only the thickness of the heterogeneous blocks when the maximum and minimum acceptable density values were fixed (marked by * in Table 1).

Calibration of Model 2, Model 2: ρmin and Model 2: ρmax
While the calculated gravity anomaly fits the observed patterns for Model 2 with σm = 4.7, for Model 2: ρmin with σm = 5.8 and for Model 2: ρmax with σm = 5.5, slight discrepancies were still visible in some intervals (Figure 4). Hidden bodies beneath the basement surface, the low sensitivity of the large-scale gravity survey to small bodies and/or the effect of density anomalies outside of the studied 2D object could explain such discrepancies.
Once the Bouguer anomaly was fitted, the remaining calibration parameters were rechecked:


The stratigraphy, porosity and pressure remained unchanged.  The Moho depth stayed almost unchanged (Figure 3). A negligible difference with respect to Model 1 was observed only in rift-graben zones.  Table 1. The corresponding modelled gravity anomalies are presented in Figure 4. The layer of the lithospheric mantle is not displayed for the convenience of presenting the picture.

Construction of Models with the Heterogeneous Basement:
Model 2, Model 2: ρ min and Model 2: ρ max Model 1 was used as a precursor for constructing models with a heterogeneous basement. The upper part of the crust was replaced with the heterogeneous basement (Figure 2b), where each block's present-day lateral borders were defined using the pre-Jurassic geological map ( Figure 1c) and linked to the profile (Figure 1b). The full version of the maps are accessible via the web links in [44][45][46]. The material properties from Table 1 were assigned to each block, and the thicknesses of the heterogeneous basement blocks (Figure 2c) were reconstructed with the help of gravity anomaly fitting (Scheme A1). For Model 2, the density values of units 8 and 18 were manually adjusted from 2870 and 2740 to 2840 and 2840 kg/m 3 , respectively, (Table 1) according to Schön [49] (p. 116).
The uncertainty of the determined density of the blocks in the basement may affect the thickness of the reconstructed blocks, impacting the entire thermal regime reconstruction. We considered the uncertainty of the density only for igneous rocks since they have a critical impact on thermal history. We built Model 2: ρ min and Model 2: ρ max using the range of acceptable densities for igneous rocks taken from Schön [49] (p. 116) for units 1, 2, 3, 4, 7, 8, 15, 16, 18 and 19. The basement geometry was obtained following the reduced-rank procedure (Scheme A1), varying only the thickness of the heterogeneous blocks when the maximum and minimum acceptable density values were fixed (marked by * in Table 1).

Calibration of Model 2, Model 2: ρ min and Model 2: ρ max
While the calculated gravity anomaly fits the observed patterns for Model 2 with σ m = 4.7, for Model 2: ρ min with σ m = 5.8 and for Model 2: ρ max with σ m = 5.5, slight discrepancies were still visible in some intervals (Figure 4). Hidden bodies beneath the basement surface, the low sensitivity of the large-scale gravity survey to small bodies and/or the effect of density anomalies outside of the studied 2D object could explain such discrepancies.
Once the Bouguer anomaly was fitted, the remaining calibration parameters were rechecked:

•
The stratigraphy, porosity and pressure remained unchanged.

•
The Moho depth stayed almost unchanged (Figure 3). A negligible difference with respect to Model 1 was observed only in rift-graben zones.

•
The thermal regime was refined by the e-fold length parameter (A r ), which was slightly increased to 21 km (to compensate the decreased radiogenic heat production and thermal conductivity in the heterogeneous basement).
The modelled geological section of Model 2 is presented in Figure 6b, and those of Model 2: ρ min and Model 2: ρ max are presented in Figure 7. The thermal regime was refined by the e-fold length parameter (Ar), which was slightly increased to 21 km (to compensate the decreased radiogenic heat production and thermal conductivity in the heterogeneous basement).
The modelled geological section of Model 2 is presented in Figure 6b, and those of Model 2: ρmin and Model 2: ρmax are presented in Figure 7.  Table 1. The layer of the lithospheric mantle is not displayed for the convenience of presenting the picture.

Results and Discussion
According to the calibration results obtained using the Moho boundary depth ( Figure  3), porosity and pressure, present-day temperatures and vitrinite reflectance ( Figure 5), we concluded that all the models (Model 1, Model 2, Model 2: ρmin and Model 2: ρmax) were acceptable. However, the additional calibration by gravity data demonstrated an oversimplification of Model 1 that affected the modelling of the thermal regime.

Thicknesses of Basement Heterogeneity Blocks
Deep-rooted vertically elongated bodies of magmatic origin mentioned in [42,76,77] are clearly visible at the reconstructed geological section of the basement beneath the rifted zones (coloured violet in Figure 6b). We associated their origin with the latest Palaeozoic to Early Mesozoic times of the global plate boundary and plate kinematics reorganisation [43]. Thus, the magmatic bodies originate from voluminous intracontinental flood basalts.
The basement geometry is in good agreement with studies conducted by Braitenberg and Ebbing [42], who estimated the thickness of intrusive igneous rocks based on gravity and magnetic anomalies of up to 5 km beneath the eastern rift grabens area at around 620-690 km and 720-750 km of the profile intervals. The thickness of the pre-Triassic paleobasin is in good agreement with the interpretation of the seismically derived bottom section of the neighbouring regional seismic profile. The thickness of the paleobasin outside rift-graben zones increases eastward from less than 1 km to ~3.5 km [38]. The modelled heterogeneous blocks' thicknesses matched the minimal thickness of blocks determined by drilling in adjacent territories [45,46]. On the 256-574 km profile interval, the thickness of units 11, 12 and 13 was more than 66, 250 and 1226 m, respectively [45]. On the profile interval of 574-808 km, the thickness of units 13, 15, 17, 18, 20 and 21 was more than 320, 1500, 500, 200, 380 and 100 m, respectively [46].  Table 1. The layer of the lithospheric mantle is not displayed for the convenience of presenting the picture.

Results and Discussion
According to the calibration results obtained using the Moho boundary depth (Figure 3), porosity and pressure, present-day temperatures and vitrinite reflectance ( Figure 5), we concluded that all the models (Model 1, Model 2, Model 2: ρ min and Model 2: ρ max ) were acceptable. However, the additional calibration by gravity data demonstrated an oversimplification of Model 1 that affected the modelling of the thermal regime.

Thicknesses of Basement Heterogeneity Blocks
Deep-rooted vertically elongated bodies of magmatic origin mentioned in [42,76,77] are clearly visible at the reconstructed geological section of the basement beneath the rifted zones (coloured violet in Figure 6b). We associated their origin with the latest Palaeozoic to Early Mesozoic times of the global plate boundary and plate kinematics reorganisation [43]. Thus, the magmatic bodies originate from voluminous intracontinental flood basalts.
The basement geometry is in good agreement with studies conducted by Braitenberg and Ebbing [42], who estimated the thickness of intrusive igneous rocks based on gravity and magnetic anomalies of up to 5 km beneath the eastern rift grabens area at around 620-690 km and 720-750 km of the profile intervals. The thickness of the pre-Triassic paleobasin is in good agreement with the interpretation of the seismically derived bottom section of the neighbouring regional seismic profile. The thickness of the paleobasin outside rift-graben zones increases eastward from less than 1 km to~3.5 km [38]. The modelled heterogeneous blocks' thicknesses matched the minimal thickness of blocks determined by drilling in adjacent territories [45,46]. On the 256-574 km profile interval, the thickness of units 11, 12 and 13 was more than 66, 250 and 1226 m, respectively [45]. On the profile interval of 574-808 km, the thickness of units 13,15,17,18,20 and 21 was more than 320, 1500, 500, 200, 380 and 100 m, respectively [46].

Densities of the Basement
A significant underestimation of density heterogeneity within the basement was observed in Model 1 compared with Model 2 (Figure 8), especially beneath the rift graben zones. The granites of the upper crust of Model 1 were less dense than in Model 2 with regard to basalt rocks formed during mega-rifting. The difference in the basement density distribution between Models 1 and 2 resulted in an insignificant change in the Moho boundary position, up to 450 m. The density distribution in the model by Isaev [78] obtained for a sub-parallel 350 km-long seismic profile located about 100 km north from our model demonstrated a similar lateral increase in density up to 150-300 kg/m 3 beneath the rift structure in its western part (around 80-120 km of the profile interval).

Thermal Regime
Generally, the effective thermal conductivity and radiogenic heat production are lower in the heterogeneous basement domain because the substituted upper crust unit has higher thermal conductivity and radiogenic heat production (Table 1) than the heterogeneous basement's igneous and sedimentary rocks. As a result, in Model 2, the value of radiogenic heat production in the crust increased by the e-fold length parameter (Ar) from 20.4 to 21 km, in order to compensate for the difference in the basal heat flow in the well thermal calibration points. The most significant divergence between the results of Model 2 and Model 1 was observed beneath the rift graben zones in the so-called magmatic roots where thermal wellbore calibration was absent. In this area, in Model 2, the bulk thermal conductivity and radiogenic heat production were 2.0-2.2 W/m/K and 1 μW/m 3 , respectively, while in Model 1, they were 2.2-2.7 W/m/K and 1-2 μW/m 3 , respectively. The difference between the magmatic roots' thermal properties and the crustal layer led to anomalies in temporal variations in the basal heat flow (Figure 9a,b).
In Model 2, the present-day heat flow values in the rift graben zones of the Vartov High (the eastern graben on the profile) are lower by 5 mW/m 2 relative to the surrounding areas (Figure 9b), which is entirely consistent with Kurchikov's [79] geothermal studies. Figure 10a shows the difference in the basal heat flow as a function of time for Models 1 and 2 caused by the basement's thermal properties. The time-averaged difference in heat flow (Figure 10b) was predominantly positive since heterogeneous basement blocks had lower thermal conductivity and radiogenic heat production values than the correspond-

Thermal Regime
Generally, the effective thermal conductivity and radiogenic heat production are lower in the heterogeneous basement domain because the substituted upper crust unit has higher thermal conductivity and radiogenic heat production (Table 1) than the heterogeneous basement's igneous and sedimentary rocks. As a result, in Model 2, the value of radiogenic heat production in the crust increased by the e-fold length parameter (A r ) from 20.4 to 21 km, in order to compensate for the difference in the basal heat flow in the well thermal calibration points. The most significant divergence between the results of Model 2 and Model 1 was observed beneath the rift graben zones in the so-called magmatic roots where thermal wellbore calibration was absent. In this area, in Model 2, the bulk thermal conductivity and radiogenic heat production were 2.0-2.2 W/m/K and 1 µW/m 3 , respectively, while in Model 1, they were 2.2-2.7 W/m/K and 1-2 µW/m 3 , respectively. The difference between the magmatic roots' thermal properties and the crustal layer led to anomalies in temporal variations in the basal heat flow (Figure 9a,b).
In Model 2, the present-day heat flow values in the rift graben zones of the Vartov High (the eastern graben on the profile) are lower by 5 mW/m 2 relative to the surrounding areas (Figure 9b), which is entirely consistent with Kurchikov's [79] geothermal studies. Figure 10a shows the difference in the basal heat flow as a function of time for Models 1 and 2 caused by the basement's thermal properties. The time-averaged difference in heat flow (Figure 10b) was predominantly positive since heterogeneous basement blocks had lower thermal conductivity and radiogenic heat production values than the corresponding values for the crust (Figure 10). The maximum mean value of the difference corresponded to the rift graben zones and reached 10.2 mW/m 2 due to the deepest heterogeneous basement units in this location.   Unit 16 (578-590 km of the profile) had the lowest negative values (down to −1.5 mW/m 2 ) of the difference as it had a higher thermal conductivity than the crust. All the other negative values occurred due to increased radiogenic heat production impacted by the e-fold length parameter (A r ) for the crust radioactivity in Model 2 during its thermal calibration.
Since the primary source rock in the Western Siberian Basin is the Bazhenov Fm, it is necessary to evaluate the impact of basal heat flow changes on its present-day temperatures and vitrinite reflectance distribution (red and blue lines in Figure 11a,b). Discrepancies were observed between Model 1 and Model 2 estimations of present-day temperature and vitrinite reflectance at the Bazhenov Fm surface in areas with insufficient well coverage. Such discrepancies are associated with the rift graben zones (four segments along the profile with a total length of 210 km), where Model 1 provides higher temperatures and a vitrinite reflectance higher than those of Model 2 by up to~9 • C and 0.05 %Ro, respectively. The essential feature of this discrepancy is that Model 1 predicts the Main Oil Window along 75 km of the profile, while Model 2 predicts the Early Oil Window (Figure 11b) according to the definition given by Tissot and Welte [80]. Such a difference between the two predictions may lead to a significant divergence in petroleum system evaluation and may result in errors in the estimated mass generation of hydrocarbons, the critical moment, reservoir volumes (for 3D modelling) and other important parameters. The remaining discrepancies in the estimations of temperature and vitrinite reflectance profiles did not exceed 2 • C and 0.02 %Ro, respectively, which is less than the uncertainty of the calibration values (see Appendix C).

Density Uncertainties
The results obtained for Model 2: ρ max were similar to those of Model 1 in terms of the thermal regime ( Figure 12). The thin heterogeneous basement (with thickness up to 1.8 km) explains the similarity in temperature (Figure 12a) and maturity (Figure 12b) beneath the rift graben zones. The discrepancy between the models was only observed in the profile zone from 490 to 540 km, which is explained by the presence of the paleobasin thick block, with lower thermal conductivity (2.4 vs. 3.0 W/m/K) and radioactive heat production (1 vs. 2 µW/m 3 ) in Model 2: ρ max heterogeneous basement compared with Model 1 crystalline basement.
Model 2: ρ min demonstrates the lowest limit of the thermal regime governed by the thick heterogeneous basement beneath the rift graben zones (up to 27 km) since its blocks have lower values in terms of thermal properties compared to Model 1.
The difference between Model 2: ρ max and Model 2: ρ min in the temperature and maturity of the Bazhenov Fm surface reached 9 • C and 0.06 %Ro, respectively, and the difference in the basal heat flow was 5 mW/m 2 (Figure 12).
The spread between Model 2: ρ min and Model 2: ρ max in terms of temperatures and the vitrinite reflectance on the Bazhenov Fm surface increased with time, while the corresponding residual basal heat flow decreased with time ( Figure 13). At the profile location at 645 km, the spread of temperature increased from 0 to 9 • C ( Figure 13a) and vitrinite reflectance increased from 0 to 0.06 %Ro (Figure 13b). Residual basal heat flow decreased from 8 to 4.5 mW/m 2 (Figure 13c). The latter can be explained by the blanketing effect [81][82][83]. Since Model 2: ρ max has a higher contrast in the total thermal conductivities of the sedimentary cover from the basements, the decrease in the heat flow with the sedimentation rate is more significant.
Due to Model 2: ρ max and Model 2: ρ min , it was possible to determine the range of heterogeneous basement geometry and to quantify its impact on the modelled thermal regime. Such estimates are particularly relevant in the rift graben zones, where the resolution of the available seismic surveys is insufficient for distinguishing certain structural features.

Thermal Conductivity Uncertainties
The results of thermal conductivity uncertainties that impact present-day basal heat flow, surface temperature and vitrinite reflectance of the Bazhenov Fm, are presented in Figures 12 and 13. The difference in thermal conductivity values in Model 2: λ ± 20% relative to Model 2 include minor changes in the basal heat flow of up to ±1.6 mW/m 2 (around 4%), corresponding to changes in temperature and vitrinite reflectance values of ±0.9 °C and ±0.01 %Ro, respectively. The difference in thermal conductivity values of Model 2: λ ± 50% relative to Model 2 resulted in a difference in the basal heat flow of up to ±4.8 mW/m 2 (around 11%), in temperature of up to ±5 °C and in vitrinite reflectance of up to 0.03 %Ro, which can lead to erroneous estimates of the organic matter's maturity rank in some intervals of the profile (e.g., 468-485, 577-587 and 644-654 km). These values may

Thermal Conductivity Uncertainties
The results of thermal conductivity uncertainties that impact present-day basal heat flow, surface temperature and vitrinite reflectance of the Bazhenov Fm, are presented in Figures 12 and 13. The difference in thermal conductivity values in Model 2: λ ± 20% relative to Model 2 include minor changes in the basal heat flow of up to ±1.6 mW/m 2 (around 4%), corresponding to changes in temperature and vitrinite reflectance values of ±0.9 • C and ±0.01 %Ro, respectively. The difference in thermal conductivity values of Model 2: λ ± 50% relative to Model 2 resulted in a difference in the basal heat flow of up to ±4.8 mW/m 2 (around 11%), in temperature of up to ±5 • C and in vitrinite reflectance of up to 0.03 %Ro, which can lead to erroneous estimates of the organic matter's maturity rank in some intervals of the profile (e.g., 468-485, 577-587 and 644-654 km). These values may increase with the thickness of the igneous rock or when using another temperature dependence for thermal conductivity of basement blocks.
One additional study for Model 2 was conducted in order to estimate the effect of possible errors in the matrix thermal conductivity of units 8 and 16 (Table 1) Models providing different estimations of thermal conductivity uncertainties in the heterogeneous basement make it possible to determine the ultimate range of modelled thermal histories. As a result, the maturity rank of the main source rock can be underesti-mated as the Early Oil Window instead of the Main Oil Window along 127 km of the profile, that is, 16% of the total profile length and~75% of the profile length where the Main Oil Window was forecasted in Model 1.

Conclusions
A method was presented for the express estimation of the impact of basement heterogeneity on thermal history and petroleum system modelling.
The method was applied to the case study of the Western Siberian Basin. Reducedrank inversion was added as an extension to the classic basin modelling workflow, and provided a rough model of deep-rooted bodies and the remnants of the pre-Triassic basin, which is otherwise impossible to create as seismic data, as the basement's domain cannot be distinguished in seismic models [37]. The introduction of magmatic bodies in the model over rift graben zones significantly changed the basin's thermal history reconstruction results and forecast kerogen maturity in the rift graben zones. We found that the maturity rank of the Bazhenov Fm can be underestimated as the Early Oil Window instead of the Main Oil Window. The uncertainties in basement characteristics (density, thermal conductivity) result in uncertainties in the temperature and vitrinite reflectance of the Bazhenov Fm in the rift-graben zones, of up to ±5 • C and 0.03 %Ro, relative to the model without the heterogeneous basement.
A full-rank 3D model of the Western Siberian Basin should be constructed accounting for the geometry and properties of the basement in the rift graben zones. Thus, model construction requires additional geophysical studies, such as deep seismic sounding, the building of seismic models, an analysis of gravity and magnetic anomalies, together with the measurement of the density and thermal properties of the rocks. Otherwise, the uncertainty range for the temperature and vitrinite reflectance of the Bazhenov Fm surface may reach 12 • C and 0.1 %Ro.
The presented method can be applied for the extensional basins without salt tectonics. The functionality of the presented method is valid for all similar datasets. The method can be considered most significant for the basin effected by the mega-rifting with deep-rooted mega intrusions.

Appendix A
The reconstruction of the basement's geometry and density presented in Scheme A1 uses an iterative approach similar to Bott's method [24]. Scheme A1 corresponds to step 3 in Scheme 1.
Firstly, when the above-mentioned assumptions about the basement properties and lateral extension were introduced in the model, the first approximation (z) for each (i = 1, . . . n) block's thickness (h i ) was made, e.g., z = 1 km. Further, the model was simulated.
The gravity anomaly was calculated according to the formula of vertically sided two-dimensional blocks presented in the work of Heiland [25] (p. 152).
Further, the difference (∆A i ) between the observed (A obs ) and calculated (A calc ) gravity anomalies was estimated for each i block.
The stopping condition of the inverse problem was used as the standard error, as in [84,85] where σ is the desired accuracy. Some deviations of the modelled gravity anomaly from the measured gravity anomaly may result from compositional differences (e.g., [86]) or 3D density structures that are not considered in the model and may contribute to the uncertainty in the 2D gravity modelling [19]. Otherwise, block thicknesses or, in rare cases, densities need to be refined. The reconstruction of the basement's geometry and density presented in Scheme A1 uses an iterative approach similar to Bott's method [24]. Scheme A1 corresponds to step 3 in Scheme 1.
Firstly, when the above-mentioned assumptions about the basement properties and lateral extension were introduced in the model, the first approximation (z) for each (i = 1, … n) block's thickness (hi) was made, e.g., z = 1 km. Further, the model was simulated.
The gravity anomaly was calculated according to the formula of vertically sided twodimensional blocks presented in the work of Heiland [25] (p. 152).
Further, the difference (ΔAi) between the observed (Aobs) and calculated (Acalc) gravity anomalies was estimated for each i block.
The stopping condition of the inverse problem was used as the standard error, as in [84,85], where σ is the desired accuracy. Some deviations of the modelled gravity anomaly from the measured gravity anomaly may result from compositional differences (e.g., [86]) or 3D density structures that are not considered in the model and may contribute to the uncertainty in the 2D gravity modelling [19]. Otherwise, block thicknesses or, in rare cases, densities need to be refined.
Scheme A1. Inversion procedure of the geometry and density of a heterogeneous basement. Scheme A1. Inversion procedure of the geometry and density of a heterogeneous basement.
Next, the residual error of thickness (∆h i ) was calculated for each i block by the formula from [24].
where G is the gravity constant. Next, the new thickness (h i new ) of the block could be calculated as where ε is the regularising multiplier ranging from 0 to 1 defined by the user. If the new block thickness (h i new ) is greater than the maximum (h i_max ) or less than the minimum (h i_min ) determined thicknesses for the i block, correction of the i block density (ρ i ) is required where, k is a correction factor ranging from 0 to 1 defined by the user and γ is a multiplier equal to +1 if h i new = h i_min and −1 if h i new = h i_max . In addition, h i_max could be limited, for instance, by h i_crust , and h i_min could be limited by drilling data Further, after the correction of values h i and ρ i from 1 to n (as recommended in [24]), the basin model was simulated and the gravity anomaly was calculated. The inversion procedure was rerun while condition (1) was not completed.
There are cases when gravity data fitting for the block is completed along only the part of the block. Therefore, this block needs to be split vertically into several sub-blocks. The ∆h i for each sub-block should be calculated separately. This procedure helps to fit the gravity anomaly better [24] and find the ultimate solutions for further estimations in thermal regime changes. * Here, ϕ 0 is the surface porosity and B is the scale depth factor in the equation ϕ(z) = ϕ 0 × exp(B × z), where z is a negative value of the depth (km); ρ is density; λ is thermal conductivity at 20 • C; C ρ is the specific heat capacity; and A is radiogenic heat production of the rock matrix.

Appendix C
Here, we provide information on the uncertainty related to the data on vitrinite reflectance, formation temperature and gravity data. When calibrating the model to the measured data, the value of uncertainty in the calibration data should be considered. This aspect is often omitted in practice; hence a brief analysis (a comprehensive study is a separate study theme) is included here.
For vitrinite reflectance, we deal with the uncertainty in the experimental determination %Ro, considering the kinetics model used for %Ro calculation. Often, if experimental data are absent, specialists use different standard kinetic models (e.g., [73][74][75]). The choice of model must be supported by the best simultaneous fitting of the modelled vitrinite reflectance and temperature with the massive set of measured data. If a large set is absent, models with different kinetics should be considered. An illustrative example of the difference in applying different kinetic models to evaluate maturity is presented in [74]. A comparison of the BasinRo model and the widely used EasyRo model shows that the EasyRo model can overestimate vitrinite reflectance in the interval of 0.5-1.7 %Ro by up to 0.35 %Ro [74].
Within the same field in identical environments, the Ro values often differ by 10-15%. One of the reasons is the insufficient number of measurements at each micro-inclusion of vitrinite (10-15 instead of 100) [87]. Another reason is the value of uncertainty in determining %Ro, depending on the method and equipment used. For measuring the reflectance of macerals, it is crucial to know what parameter of the vitrinite reflectance value is determined-minimal (R min ), maximal (R max ) or random (R r ) [88][89][90]. However, this is not often mentioned in publications. With the example of Upper Silesian coal, Komorek and Morga [91] demonstrated that the standard deviations in R max and R r are similar and do not exceed 0.05 %Ro for the range of maturity 0.2-0.9 %Ro. For maturity higher than 0.9 %Ro, the standard deviation increases and R max becomes more precise than R r (up to 0.13 %Ro vs. 0.33 %Ro) and, therefore, acts as a more suitable rank indicator. The standard deviation of vitrinite reflectance measured for dispersed macerals in sedimentary rocks also increases with advancing thermal maturity and increasing anisotropy of the organic matter. In [92], vitrinite reflectance for different types of organic matter (by the method in [93]) is demonstrated, where an increase in standard deviation from 0.04 to 0.19 %Ro in the range of 0.31-1.53 %Ro (without mentioning what method of Ro determination has been used) is observed. Houseknecht and Matthews [94], using the example of R r measurements of dispersed matter in carboniferous strata, Ouachita Mountains, established that standard deviation σ = 0.083 + 0.145 × R r (correlation coefficient = 0.808, number of samples = 89). Summarising the review above, for the measured values of vitrinite reflectance 0.45-0.72 %Ro (Figure 5b), the standard deviation in the range 0.01-0.02 %Ro per [91], 0.06-0.09 %Ro per [92] and 0.15-0.19 %Ro per [94] can be estimated.
The uncertainty in the indirect determination of %Ro, e.g., via the recalculation of T max from pyrolysis data, includes the uncertainty of the pyrolysis measurement temperature (about ±2 • C with Hawk and Rock-Eval 6 tools according to [95,96]) and the uncertainty caused by the formula used for recalculation. For example, the widely (see, e.g., [97]) used formula %Ro = 0.018 × T max − 7.16 by Jarvie et al. [98] implies standard deviation σ = 0.23 %Ro at N = 179 (N is the number of samples), i.e., the uncertainty is much higher than the above-described uncertainty resulting from direct measurements.
Uncertainty related to formation temperature includes the instrumental uncertainty of measurement and the uncertainty of the non-equilibrium formation temperature. The first depends on the equipment used, and according to Blackwell and Spafford [99], it is approximately ±0.5 • C for commercial temperature logs and ±0.05 • C for scientific equipment. The second is caused by uncertainty with undisturbed temperature measurements, as the temperature readings in a well during drilling or completion do not reflect the equilibrium formation temperature [100,101], and the measured temperature values can be underestimated by up to~22 • C, with standard deviation σ = 10 • C for depths of 500 to 9000 m [102,103]. Various correction schemes have been presented and discussed in the literature (e.g., [104,105] and references therein); however, the uncertainty is rarely estimated. In [102], the application of corrections improved "a prediction of formation temperatures with an error less than ±10 • C at somewhat deeper depths than the log in several boreholes." The value of the total error in the Bouguer anomaly [106,107] takes into account several components: an error of reference network observation, ε rn ; an error of accounting for terrain correction, ε r ; an error in introducing corrections for free air (depends on the error in the station heights), ε F ; an error in introducing the Bouguer correction (includes the error in determining the topographic relief and average density of rocks), ε B ; and an error in calculating the normal gravitational field (associated with the error in calculating the normal value of the acceleration of gravity and determining the coordinates of observation points, ε γ (see details in [70,106,108,109] and links therein). The total error of the Bouguer anomaly must have standard limitations for different scales of maps [70]; the total error should not exceed ε t = ±0.4 mGal for the 1:5000 scale and ε t = ±1.5 mGal for the 1:500,000 scale, which is consistent with [108].