Estimation of the Oil-in-Place Resources in the Liquid-Rich Shale Formations Exploiting Geochemical and Petrophysical Data in a 3D High-Resolution Geological Model Domain: Baltic Basin Case Study

The paper discusses the issue of oil-in-place estimation for liquid-saturated shales in Lower Paleozoic (Silurian and Ordovician) organic-rich formations of the Baltic Basin in North Poland. The authors adopted a geochemical method based on Rock Eval results which directly measure hydrocarbon content present in rock samples. Its application on a real data set required the implementation of correction procedures to consider also those oil compounds which were lost before Rock Eval measurements were taken or are not recorded in S1 parameter. It was accomplished through the introduction of two correction coefficients: c1—for evaporation loss and c2—for heavier compounds underestimation. The first one was approximated on the basis of published results and known properties of crude oil, while the second one was addressed with laboratory experimental procedure which combines Rock Eval pyrolysis and rock sample extraction with organic solvents. The calculation formulas were implemented in the 3D geological model of shale formations reproducing their geometry as well as the spatial variability of the petrophysical and geochemical properties. Consequently, the results of oil-in-place estimation were also available as 3D models, ready for visualization and interpretation in terms of delineation of most favorable zones or well placement. The adopted geochemical method and the results of oil-in-place estimation it produced were confronted with standard volumetric method. Although both of them are volumetric methods, the results depend on different sets of rock properties, which is an advantage for result comparison reasons. The study revealed that the geochemical method of oil-in-place estimation in liquid-rich shales after appropriate adjustment, considering shale formation and reservoir fluid dependent conditions, could provide reliable results and be implemented on the early stage of shale exploration process in a condition of production data inaccessibility.


Introduction
Unconventional shale reservoirs exploration and development resulted in a development of new approaches for interpretation and modeling of reservoir characteristics which takes into consideration specifics of shale formations, its petrophysical and geochemical properties. The issues which also require the adoption of a different approach than conventional reservoirs are resource assessment methods. Literature references regarding the methods for estimation of unconventional resources accumulated in shale formations (gas shales/oil shales) report the existence of several methodological approaches [1][2][3][4][5][6][7][8][9] which can be employed depending on the stage of shale resource exploration or development. Generally, the methods for shale formation resource assessment could be divided into several categories, i.e.: (i) Volumetric methods (based on shale rock petrophysical or geochemical properties and PVT data) (ii) Material balance methods (dynamic petroleum system models) (iii) Production history (based on production data analysis) (iv) Analogy (based on comparative analysis with similar, well-recognized formations) Each of the abovementioned categories requires the availability of a particular data set which in turn depends on the advancement of exploration and development activities in a given shale play. Up to date, except for several most developed shale provinces, mainly in the USA, Canada, and China, the application of methods requiring production data (such as EUR based method) is premature due to the lack of any production history. On the other hand, assessments made on analogies to US gas and oil-producing shale formations turned out to be unrealistic (e.g., estimates prepared for Lower Paleozoic shale formations in Poland), providing overly optimistic outlooks which were not confirmed by the later activity of exploration companies. Therefore, it is a great need for reliable methods to assess oil or gas resources in unconventional shale formations at an early stage of exploration.
In the paper, we presented results of the research project aimed at shale formation characterization in the Baltic Basin in Poland which was accomplished in recent years. One of the challenges considered adaptation and application of geochemical method for oil-in-place estimation of Lower Paleozoic liquid-rich shale formations which have not yet reached the production phase. The specifics of shale formations saturated with liquid hydrocarbons make it possible to undertake the task of oil-in-place assessment taking advantage of direct geochemical measurements of hydrocarbons contained in rock samples collected from liquid-rich shale interval. The laboratory method which allows to assess liquid hydrocarbon content in rock samples and is relatively straightforward as well as inexpensive is Rock Eval (RE) pyrolysis [10,11]. The S1 parameter, resulting from RE measurement, reflects "free" hydrocarbons contained in a rock sample and is the basis for calculation procedure. However, the S1 parameter refers only to the part of hydrocarbons present in reservoir fluid (here: oil with dissolved gas) stored in reservoir conditions. The core sample handling procedure preceding laboratory analysis deprives a sample of the lightest compounds: C1-C4 (gas dissolved in oil) and light liquid hydrocarbons evaporating completely or partially in surface conditions (up to C15). Only hydrocarbons with a carbon number of C15+ up to C25 are reflected in the S1 parameter, which records components evaporating in a pyrolytic oven at the temperatures of up to 300°C. The heavier liquid compounds (C26-C35+) are released at the temperature exceeding 300°C and are recorded as a part of the S2 parameter together with kerogen residues in a temperature range of up to 500°C.
The core of the methodological approach which we used to estimate oil-in-place in a study area consists of an introduction of appropriate correction coefficients expressing the scale of underestimation of the lightest as well as heavier components of crude oil in S1 parameter values as received from Rock Eval measurements. The first one (c1-correction coefficient) takes into account the evaporation process associated with core sample handling and is intended to reflect the loss of lightest hydrocarbons (C15-). The second one facilitates the underestimation of heavier crude oil compounds in S1 values. The value of the correction coefficient for evaporation loss was established on the basis of published data [12], while the correction needed to facilitate the underestimation of heavier components of crude oil was elaborated with laboratory investigation involving Rock Eval pyrolysis and organic solvent extraction method. Each step of correction implementation is followed by an appropriate equation for oil-in-place estimation, including the approximation of free oil, "immobile" or "sorbed" oil, and finally the total oil-in-place.
Apart from methodological studies and laboratory experiments, 3D high-resolution geological models of shale formations being analyzed were developed to enable oil-in-place calculations in the 3D model domain. Hereby, the resolution of the estimation outcome corresponds to dense core sampling and resulting availability of shale formations' properties including Rock Eval data.
The most important issue of each resource or reserve assessment method is its credibility. The geochemical method of oil-in-place estimation in shale formations adjusted and applied in this study was confronted with a standard volumetric method in terms of oil-in-place estimates they produce. Both methods are volumetric; however, they are distinguished by an unlike data requirements. While geochemical methods rely on Rock Eval data, the standard volumetric method is based on porosity and fluid saturation characteristics of a shale formation and reservoir fluid properties. This makes them very suitable for result comparison; however, they could be also used alternatively (according to data availability) or in parallel to get some insight into the uncertainty of resource assessment.
1.1. Geological Background and Data Set Available. The study area is located in northern Poland, and western fragment of Peribaltic Syneclise, which is a part of the Baltic Basin, being an oval recess in the Precambrian formations of the Eastern-European platform [13,14] (Figure 1). Due to the growth of interest in an unconventional shale reservoir in the recent 10 years, several dozen exploration drillings have been carried out in this area with an extensive core collection and borehole logging program. One of the outcome is a relatively wide set of new data providing the basis for more in-depth research and methodology evaluation possibilities.
Shale reservoir exploration in Baltic Basin was focused on Lower Paleozoic (Silurian and Ordovician) complexes which comprise the following formations [15]  of liquid hydrocarbons and is consistent with thermal maturity of major producing shale oil plays in US [19][20][21]. The method proposed in this paper for oil-in-place estimation in liquid-saturated shale formations was applied to the Jantar and Sasino shale formations for which the 3D geological model was developed. For this purpose, an extensive data set was integrated, including laboratory data (petrophysical and geochemical) and well log interpretations from four vertical exploration boreholes (Figures 1 and 2). Laboratory investigation into geochemical and petrophysical properties of Lower Paleozoic formations in the study area included Rock Eval pyrolysis, mercury injection porosimetry, X-ray diffraction analysis (XRD), nuclear magnetic resonance (NMR) measurements, permeability investigation with Pulse Decay Permeameter (micro-and nanoscale), and Dean-Stark measurements. The results of the abovementioned investigations were presented in detail and discussed in a number of publications and conference papers [17,[22][23][24][25][26][27]. In this article, we do not present a detailed geological study but rather focus on the calculation method and its application in the domain of the 3D geological model developed on the basis of the data set presented above. In this method, the 3D high-resolution reservoir model is a tool of result visualization and interpretation of shale oil resource distribution within formations of interest.

3D High-Resolution Geological Models of Shale
Formations Studied. The geological 3D model of the study area was developed on the basis of available data sets consist-ing of laboratory measurements, well logs, and seismic data (structural and seismic attributes) [27,28]. Geostatistical, variogram-based methods of rock property spatial variability modeling was applied [29,30]. 3D distributions of the following petrophysical and geochemical properties of the Jantar and Sasino formations have been modeled: shale rock bulk density, total porosity and effective porosity, organic matter content (TOC), saturation of the pore space of shale rocks with reservoir water (Sw), and the S1 parameter, indicating the free hydrocarbon content of the rock sample undergoing pyrolysis ( Figure 3). The correlation of the S1 parameter with TOC values and TOC with rock density was quantified and used in reservoir modeling procedures. In that way, the 3D grid population with TOC, and finally S1, values were indirectly supported by 3D seismic-derived rock density cube.
Apart from shale formation properties, procedures of oil-in-place calculation require the knowledge of the crude oil properties (Table 1). It applies to standard volumetric methods based on porosity and fluid saturation properties, which were used in this study to quality check the geochemical method.
1.3. Methodology: Adjustment of Geochemical Method Based on Rock Eval Data. The specifics of shale formations saturated with liquid hydrocarbons make it possible to undertake the task of oil-in-place assessment taking advantage of direct geochemical measurements of hydrocarbons contained in rock samples collected from liquid-rich shale interval. An advantage of the oil resource assessment method based on TOC RHOB S1 Figure 3: Reservoir modeling results-3D distributions of properties used for oil resource calculation. For the application of geochemical method, the only properties needed are rock density (RHOB) and S1 parameter; however, modeling task included also other properties: TOC-which was used to support S1 parameter modeling, porosity, and water saturation (Sw), which were used for geochemical method validation with standard volumetric method.
4 Geofluids geochemical data results from the fact that it does not require the availability of porosity and fluid saturation data which are not always available. It applies mainly to fluid (oil) saturation information, usually resulting from well log data interpretation, which could be ambiguous in the case of shale formations. Available laboratory methods to assess the amount of liquid hydrocarbons present in rock samples include [1,6,12,33,34]: (i) the extraction of hydrocarbons contained in the sample (EOM-extractable organic matter) by leaching them using organic solvents and determination of extract amount after solvent evaporation at the temperature of 60°C; the EOM measures C12+ components, including nonvolatile hydrocarbons (ii) pyrolytic analysis (Rock Eval) of shale rock samples-an assessment of liquid hydrocarbons is based on S1 parameter expressing the weight of free hydrocarbons contained in rock sample (mg HC/g rock) which are released in a pyrolytic oven at temperatures of up to 300°C; S1 parameter measures volatile hydrocarbons in the C5-C25 carbon number range The Rock Eval pyrolysis and S1 parameter determination is much more convenient than solvent extraction, since it is a fast and straightforward method which can be introduced early in the exploration process. It can easily be scaled up to reflect oil-in-place potential for a given shale formation thickness. In this study, the authors considered Rock Eval S 1 data as the primary measurement type for oil-in-place estimation methodology being adopted, and solvent extraction procedures were used for the elaboration of one of the correction coefficients.
The calculation procedure of oil resources based on Rock Eval S1 data for shale formation of a specified volume originates from the following formula [1,2,10,11]: where Oil content is expressed in ðm 3 of oil/m 3 of rockÞ; V o is the volume of oil (m 3 ); V r is the volume of rock (m 3 ); S1 is the oil content (mg HC/g of rock) present in rock volume determined with Rock Eval measurement; ρ r is rock density (kg/m 3 ); and ρ o is oil density (kg/m 3 ).
To employ this formula for oil-in-place estimation, it is necessary to address those amounts of liquid hydrocarbons which are not recorded in S1 values. As mentioned above, the underestimation of hydrocarbons present in rock sample with S1 parameter (compared to in situ conditions) results from (i) the evaporation of the lightest compounds of crude oil-published data refers to the loss of up to C15compounds [1,12]: however, partial evaporation of C12-C15 compounds is observed (Figure 4). The loss of lighter hydrocarbons is a consequence of rock samples handling in the RE laboratory procedure, which requires, among others, its grinding (ii) temperature conditions while S1 peak is recorded, which is up to 300°C; hydrocarbons released in a pyrolytic oven in a higher temperatures are not included in S1; however, they represent a fraction of S2 peak of Rock Eval pyrogram [6].
In the methodology we used, those sources of inaccuracies in oil-in-place estimation (if carried out only on the basis of S1 values directly) are facilitated by the introduction of appropriate correction coefficients: c1-correction coefficient for loss of C15-compounds which is dependent on oil API gravity and shale rock maturity level [1,12], c2-correction coefficient for underestimation of heavier compounds which are included in S2 peak (fraction of S2), but through the introduction of c2 coefficient they are expressed as a fraction of S1 values. The scale of correction needed which would compensate for the lightest compound loss due to evaporation process (c1 value) was approximated on the basis of published data [12], which relates liquid hydrocarbon loss due to evaporation as a function of crude oil API gravity and time of its exposure (Figure 5 below).
The estimation of oil-in-place reserves taking into account evaporation loss can be expressed by the introduction of the c1 coefficient, according to the formula: where OIIP Free is oil resources (m 3 ) calculated for given shale rock volume, which approximate "free"/"mobile" oil, S1 is the oil content (mg HC/g of rock) present in rock samples determined with Rock Eval pyrolysis method, V r is the volume of rock (m 3 ), ρ r is the rock density (kg/m 3 ), ρ o is the oil density (kg/m 3 ), and c1 is the correction coefficient for loss of C15-dependent on oil API gravity and maturity level. Oil-in-place resources calculated in that way include hydrocarbons measured in Rock Eval analysis and expressed in S1 values, together with lighter oil compounds lost due to the evaporation process. Generally, this amount would presumably represent "free" or "mobile" oil.
For more accurate assessment of total oil-in-place values, additional correction is needed, as S1 parameter does not take into account hydrocarbons contained in the shale rock which are released in the pyrolytic oven at the temperatures exceeding 300°C, while the heavier fractions 5 Geofluids of crude oil (C25+) are released at higher temperatures (up to c.a. 500°C [6,11]. In this study, the correction for heavier compounds of crude oil was addressed by the development of laboratory procedure which combines Rock Eval pyrolysis results and rock sample extraction with organic solvents. It consists of a double determination of liquid hydrocarbon content of the sample (previously divided into two parts) with Rock Eval pyrolysis, while the first analysis is carried out before hydrocarbon extractions and the second one after extraction with organic solvents. The interpretation of the results includes S1 and S2 peaks of Rock Eval pyrolysis ( Figure 6) and is aimed at determining the portion of S2 parameter which can be attributed to liquid hydrocarbons. The S2 parameter indicates the amount of hydrocarbons released within the temperature range of 300 to 650°C. Usually, it is 27 [12]; in a red rectangle crude oils collected in the study area and its vicinity are presented with the average API = 43.
The assumption was made that the time of core sample handling before Rock Eval analysis exceeds 120 minutes.
6 Geofluids interpreted as a reflection of the residual potential of source rock which is capable of hydrocarbon generation. However, for more mature source rock values of S2 parameter could be partially associated with heavier fractions of already generated hydrocarbons. This part of oil is represented by hydrocarbons in a range of C18-C35, for which boiling point exceeds 300°C. It is assumed that this fraction of oil is essentially trapped in isolated pore space, mainly within a kerogen and could be treated as "sorbed" (adsorbed and absorbed) or "immobile" oil [1,6]. The adsorption is an important storage mechanism in shales. Organic-rich shales contain a certain percentage of microporous kerogen, which exhibits high adsorption potential due to its significant internal surface area [19,35]. Furthermore, in liquid-rich shale plays, the hydrocarbons retained within the kerogen play a key role in estimated ultimate recovery (EUR) [36]. Thus, shale resource assessment methods could not ignore the adsorbed component of hydrocarbon-in-place. In the case of shale gas reservoirs, the quantification of adsorbed component is based on adsorption isotherm development, such as the most known Langmuir isotherm, which represents the adsorption of a monolayer of molecules on the inner surface of the material [37]. The studies focused on the oil adsorption in shales [38] found that there are always multiple adsorbed layers of liquid hydrocarbons. Thus, the assumptions of the Langmuir isotherm are not valid for liquidrich or wet shale formations [35].
As mentioned above, it is believed that sorbed oil is mostly constituted by higher molecular weight asphaltenes and resins [1,6]. Based on that assumption, we have made an attempt to approximate the amount of sorbed component by Rock Eval pyrolysis; in this paper, we denoted this quantity as OIIP Immobile . The amount of additional liquid hydro-carbon compounds (some residues of C18-C25 which did not evaporate totally in a temperature range of up to 300°C and C26-C35) measured using the two-step pyrolysis (before and after solvent extraction) is equal to where S2 0 is the value of S2 parameter before bituminous extraction (mg HC/g rock); S2 E is the value of S2 parameter after bituminous extraction (mg HC/g rock); and ΔS2 is the difference between S2 values before and after extraction. The implementation of the abovementioned correction leads to the formula for heavier liquid hydrocarbons resources calculation: where OIIP Immobile is the oil resources (m 3 ) of heavier fraction of liquid hydrocarbons compounds calculated for given shale rock volume. The calculation formulas presented above produces an assessment of heavier fractions of liquid hydrocarbons in place (nC25-nC35). However, conducting such a procedure is laborious, relatively expensive, as well as time-consuming as it requires the repetition of pyrolysis analysis as well as execution of bitumen extraction with organic solvents for a certain amount of shale rock samples.
For that reason, we have made an attempt to define the function linking the fraction of S2 parameter which reflects heavier liquid hydrocarbons (ΔS2-difference between S2 values before and after extraction) with any parameter received from Rock Eval pyrolysis. This would allow to limit laboratory effort of extraction with an organic solvent and   Figure 7). Relatively significant correlation between values of S10 parameter obtained by Rock Eval pyrolysis (Table 2, column 2) and heavier liquid hydrocarbon content being the difference of S2 peak before and after extraction (Table 2, column 4) was noticed (Figure 7). This creates the basis for determination of correction coefficient value for heavier oil compounds as a fraction of S1 0 , which can be expressed as where c2 is the correction coefficient for heavier oil compounds (part of S2 peak). The idea of c2 coefficient introduction is to replace S2 0 − S2 E /S1 0 expression, which requires the accomplishment of repeated pyrolysis measurements as well as solvent extraction for each sample. If c2 correction coefficient is determined on a limited number of samples, then only S1 values are needed to perform the estimation of oil-in-place resources including heavier oil compounds for the whole volume of the liquid-rich shale formation.
Thus, the calculation formula (4) could be rewritten into the following form, in which only S1 and c2 values are sufficient: where OIIP Immobile is the oil resources represented by heavier fractions (m 3 ) calculated for given shale rock volume; V r is the volume of rock (m 3 ); ρ r is the rock density (g/cm 3 ), S1 is the oil content (mg HC/g of rock) present in rock samples determined with the Rock Eval pyrolysis method; ρ o is the oil density (g/cm 3 ); and c2 is the correction coefficient for C25+ oil compounds (part of S2 parameter). Consequently, total oil-in-place calculation, including lighter (C5-C25) and heavier (C25-C35) compounds of oil can be received as a sum of "free" and "immobile" liquid hydrocarbons: Hereby, the total oil-in-place estimate is expressed with S1 parameter values, together with data on rock and oil density and established values of correction coefficients.

Quality Check with Standard Volumetric Method of
Oil-in-Place Estimation. Taking into consideration the experimental character of some elements of the Table 2: Amount of hydrocarbons measured in a following laboratory procedures: (1) S1 0 -a value of Rock Eval S1 parameter obtained in a measurement carried out for a sample before bituminous extraction (ppm); (2) EOM-amount of hydrocarbons received in bituminous solvent extraction procedure (ppm); (3) S2 0 − S2 E = ΔS2-a difference in S2 values for Rock Eval measurements carried out before and after solvent extraction (ppm); (4) proportional ratio linking heavier liquid hydrocarbons (ΔS2) with S1 0 parameter-correction coefficient (c2) for heavier oil compounds.  8 Geofluids geochemical method adjustments presented in the previous section, an attempt has been made to quality check estimation results it provides with alternative method. The standard volumetric method was used as it reveals some comparability with adopted geochemical method as both of them are volumetric methods. On the other hand, the calculations are carried out using different sets of shale rock parameters, which are also an advantage for result comparison reasons. The classical volumetric approach to oil or gas in place estimation in a conventional reservoir is based on petrophysical properties and PVT data. The free oil and gas volume in a shale formation can be computed in the same way [3]. The PVT data necessary to calculate (oil and dissolved gas volumes) are B o -oil formation volume factor, defined as the ratio of the volume of oil (liquid and gas) at reservoir conditions to the volume occupied by the same unit mass at stock tank (surface) conditions; and R s -the solution gas-oil ratio, defined as the ratio of the volume of gas dissolved in the oil to the volume of oil. For oil-saturated reservoir, the general formula of free oil in place calculation can be simplified as follows [39][40][41] where STOIIP is the volume of free oil (m 3 ); V r is the volume of rock (m 3 ); φ is the effective porosity(fraction); S w is the residual water saturation (fraction); and B o is the oil formation volume factor (dimensionless).
To calculate the volume of gas dissolved in oil following formula could be used assuming the solution gas-oil ratio (R s ) is known as [3,40] where GIIP in oil is the volume of dissolved gas in free oil (m 3 ) and R s is the solution gas-oil ratio (dimensionless). This equation could also be used for the calculation of gas-in-place dissolved in oil for estimations made with the geochemical method presented in the previous section of the paper.

Results and Discussion
The aforementioned methods of oil-in-place resources estimation have been tested on a real data set and highresolution 3D geological model, so that all calculation procedures were conducted for each block (cell) of the 3D model. It makes that the results are available for every grid cell and could be visualized and interpreted in the form of 3D models, average map of the formations of interest, or as cross-sections through the arbitrary part of the model. Figure 8 presents quantities included in calculation formulas, which we used in grid operations. The 3D model of S1 parameter, together with rock density property, as well as constant correction coefficients (c1 and c2 equal accordingly 1.62 and 0.75) produced the 3D distributions of the following resource parameters: OIIP Free , OIIP Immobile , OIIP Total , and GIIP in oil . It is evident that both OIIP Free and OIIP Immobile strongly depend on S1 values, which is a consequence of the procedures used  Figure 8: Exemplary results of oil-in-place estimation using Rock Eval-based geochemical method along one of the boreholes: on the first three tracks input data were posted-accordingly-bulk density, TOC, and S1 parameter. Tracks 3-6 present vertical profiles of OIIP Free , OIIP Immobile , and OIIP Total which is a sum of the previous ones.
9 Geofluids for correction coefficient approximation. On the other hand, in our case, the S1 parameter is strongly dependent on TOC values, which in turn exhibit a correlation with rock density (Figure 2, Figure 8).
In the case of average maps, they reflect the magnitude of oil resources within the full thickness of the given interval of interest over a unit of area (usually equal to the area of a grid cell, but it could be easily scaled to any area unit). Such a form of oil-in-place estimates presentation is convenient for the delineation of most perspective locations of the study area. Figure 9 presents the average map of OIIP Free calculated for Sasino formation-it reflects generally more favorable conditions in the west part of the analyzed area. On the other hand, as it is indicated in Figure 8, a significant variability of OIIP is observed along the vertical profiles of interpreted shale formations. In the case of Jantar formation (Shale I), an interval in the middle of its thickness presents the most favorable resource characteristics, while for Sasino formation (Shale II), its uppermost interval is the most proliferous. For a more detailed interpretation of free oil-in-place distribution within the thickness of shale formations, cross-sections through selected locations were outlined along arbitrary directions (Figures 9 and 10). The spatial form of the resource calculation enables 3D tracking of the extent of the most favorable zones, filtering, and calculation of OIIP for a particular part of the study area (e.g., approximated drainage area of production well), but it could also support horizontal sections designing ( Figure 10).
The results presented in Figures 8, 9, and 10 indicate the existence of considerable lateral and vertical diversity of the oil-in-place quantities within the shale formations being ana-lyzed. The difference between minimum and maximum values observed at the average map presented in Figure 9 demonstrates a six-fold increase. From that point of view, and assuming a correlation between oil-in-place resources and production output, an appropriate designation of borehole location, as well as horizontal section trajectory, could have a decisive impact on the profitability of the liquid-rich shale exploration and production project.
An important issue associated with the proposed geochemical method of oil-in-place estimation is its credibility. In this study, we analyzed the results it provides by comparison with oil-in-place assessment obtained by the implementation of the standard volumetric method based on porosity and fluids saturation property interpretation and PVT data. Using the standard volumetric method, only the hydrocarbon filling pore space (effective porosity) of shale formation could be calculated (in this case oil, and gas dissolved in oil). It is assumed that they represent free oil [3]. In our case study, we approximated the immobile hydrocarbons received by standard volumetric method by scaling STOIIP values by the factor resulting from the relationship between OIIP Free and OIIP Immobile as received by using the geochemical method. Nonetheless, the parameters worth comparing are estimates of free oil-in-place. Figure 11 posts the outcome of oil-in-place estimations performed with both methods (in all variants: free, immobile, and total) which were applied on a 3D grids but also well logs in existing wellbore profiles.
We observed the close correlation between the results of both methods: oil-in-place estimates follow the same vertical trends but also the magnitude of the calculated resources is similar. In the Well-1 extraction of the 3D model ( Figure 11), free oil-in-place calculated with the standard  volumetric method produces slightly larger quantities for Shale II formation, but when the whole study area is considered, the proportion is opposite ( Figure 12).
As the calculation of immobile oil in both methods represents the same fractions of free oil calculated by each method, the relationship between resulting estimates of immobile oil is the same. It refers also to the total oil estimates which is a sum of free and immobile oil.
The convergence of the results obtained with both methods is noteworthy, considering the totally different sets of input data used for the calculation procedures: in the case of geochemical method, the variables used are S1 parameter, rock density and constant values of oil density, as well as the correction coefficients. To perform the estimation with the standard volumetric method, effective porosity and water saturation variables are needed together with the PVT data (oil formation volume factor and the solution gas-oil ratio-if the quantity of gas dissolved in oil is to be estimated). The existing slight discrepancy of the oil-in-place estimates of both methods could result from the overestimation of c1 correction coefficient (defined on the basis of published data) or underestimation of porosity and/or oil saturation. However, experimental results published by Jin et al. [42,43] provide arguments for the geochemical method of oil-in-place estimation in shale oil formations. They studied the correlations between oil recovery performance using CO 2 injection in shale oil reservoirs and various factors characterizing shale formations, such as TOC, porosity, permeability, pore size, and water saturation. Their results showed that the TOC content is the most influential parameter affecting oil recovery. This suggests that the geochemical oil-in-place assessment method based on S1 values, which in our case is strongly dependent on TOC (Figure 8), will provide estimates more convergent with shale oil recovery efficiency. Another advantage of the geochemical method application for liquid-rich shales, and in particular, the implementation of correction for heavier hydrocarbons (c2 coefficient) results from its ability to roughly approximate the quantity of adsorbed oil [1,6], as opposed to the standard volumetric method. Nonetheless, the joint application of both volumetric methods creates an opportunity to consider possible sources of uncertainty, but also provides a preliminary assessment of uncertainty magnitude. Considering the alternative methods of resource assessment available in the early phase of exploration, such as analogy based analysis or petroleum systems modeling approach, volumetric methods seem to be more accurate whenever appropriate wellbore data is available.
In our opinion, the main sources of uncertainty in the geochemical method application are associated with the determination of the correction coefficient for evaporation loss which depends on crude oil API gravity-varying across basins and related to source rock maturity, as well as crude oil and core samples handling procedures. On the other hand, the value of correction coefficient for heavier compounds determined with laboratory procedures may also vary across basins as it depends on crude oil composition (fraction of polar compounds, such as asphaltenes and NSO compounds) [44,45] and sorption capacity of shales dependent on maceral composition of kerogen and its maturity [46,47].
Assuming that correction coefficients for hydrocarbon loss due to the evaporation process, as well as for heavier oil compounds not recorded in S1 parameter, are approximated with sufficient reliability, the Rock Eval-based approach would provide an inexpensive and easy to implement a methodology for oil-in-place estimation in liquidrich shale formations. Its application within the 3D reservoir modeling software allows for spatial analysis of resource distribution within target formation being analyzed. One of the most important advantages of the method results from its applicability on each phase of the exploration and production process, also the preliminary one characterized by data scarcity, particularly lack of production data.

Conclusions
An accurate estimation of oil or/and gas in place (OIIP and GIIP) is one of the priorities of the early stage of shale reservoir exploration and development. In many cases, it is accompanied by the deficiency of data and usually the lack of any production data which limits possible methodical approaches. If core data from exploration vertical wells is available, it enables considering volumetric methods, based on the application of geochemical or petrophysical data.
The geochemical method based on Rock Eval measurement results require core sample availability, while petrophysical properties could be interpreted from wireline borehole geophysics, but on the other hand, it provides a direct estimation of hydrocarbon content in shale rock sample so it is less sensitive to possible interpretation mistakes. This refers mainly to fluid (oil/water) saturation data, usually resulting  Figure 11: Comparison of oil-in-place estimation results obtained with geochemical and standard volumetric methods: the first three tracks refer to the results of the geochemical method, accordingly: OIIP Total , OIIP Immobile , and OIIP Free ; tracks 4-6 present estimates of standard volumetric method, accordingly: STOIIP Free , STOIIP Immobile , and STOIIP Total .  from well log data interpretation, which could be ambiguous in the case of shale formations. One the other hand, the study showed that calculation equations must account for substantial underestimations, as S1 parameter reflects only a part of hydrocarbon compounds originally present in shale rock at reservoir conditions. However, the value of correction coefficients could be determined with laboratory analytical procedures, enabling estimation method adjustment which could also be used in other locations of the basin, assuming similar properties of reservoir fluids as well as the level of source rocks maturity. Appropriately adjusted geochemical method could serve as a reliable, inexpensive, and straightforward approach to assess oil-inplace in liquid-rich shale reservoirs. The application of calculation procedures within the 3D geological model of shale formation allows for the accurate reproduction of oil-inplace 3D distribution which could be used for the interpretation and delineation of formations' sweet spots.
The dissimilarity of data set requested for geochemical and standard volumetric methods causes that they are very suitable for parallel application enabling result comparison, but also preliminary evaluation of uncertainty present in resource assessment. Considering the ambiguity of resource assessments in unconventional shale reservoirs, it seems reasonable to use alternative methods for oil-in-place estimation. This applies especially to new shale plays characterized by a high level of geological uncertainty, which requires knowledge-based decision-making.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.