The Effect of Hydraulic Fracture Geometry on Well Productivity in Shale Oil Plays with High Pore Pressure

We propose three idealized hydraulic fracture geometries (“fracture scenarios”) likely to occur in shale oil reservoirs characterized by high pore pressure and low differential in situ stresses. We integrate these geometries into a commercial reservoir simulator (CMG-IMEX) and examine their effect on reservoir fluids production. Our first, reference fracture scenario includes only vertical, planar hydraulic fractures. The second scenario has stimulated vertical natural fractures oriented perpendicularly to the vertical hydraulic fractures. The third fracture scenario has stimulated horizontal bedding planes intersecting the vertical hydraulic fractures. This last scenario may occur in mudrock plays characterized by high pore pressure and transitional strike-slip to reverse faulting stress regimes. We demonstrate that the vertical and planar fractures are an oversimplification of the hydraulic fracture geometry in anisotropic shale plays. They fail to represent the stimulated volume geometric complexity in the reservoir simulations and may confuse hydrocarbon production forecast. We also show that stimulating mechanically weak bedding planes harms hydrocarbon production, while stimulated natural fractures may enhance initial production. Our findings reveal that stimulated horizontal bedding planes might decrease the cumulative hydrocarbon production by as much as 20%, and the initial hydrocarbon production by about 50% compared with the reference scenario. We present unique reservoir simulations that enable practical assessment of the impact of varied hydraulic fracture configurations on hydrocarbon production and highlight the importance of constraining present-day in situ stress state and pore pressure conditions to obtain a realistic hydrocarbon production forecast.


Introduction
Development of the ultra-low permeability reservoirs, such as oil, condensate and dry gas mudrock ("shale") plays, has provided additional hydrocarbon resources to meet the world's growing demand for energy [1]. The development of these reservoirs relies on a combination of multi-stage completions and hydraulic stimulation along horizontal wells that are typically one to two miles long. Stimulation processes generate dominant propped induced fractures with permeability several orders of magnitude higher than the rock matrix, which enables and accelerates hydrocarbon production [2].
In hydraulic fracturing design, typical hydraulic fracture geometry is assumed to be elliptical, vertical with symmetrical bi-wings, and a propagation direction perpendicular to the least principal stress (S 3 ) assumed to be the minimum principal horizontal stress (S hmin ) [3,4]. However, a growing body of evidence suggests that hydraulic fracture geometries are often far more complex [5][6][7][8]. For example, interactions between induced hydraulic fractures and natural discontinuities, such as natural fractures and bedding planes, alter the dominant fracture geometry [9][10][11][12][13]. Geological and tectonic factors such as pre-existing natural fracture networks and present-day in situ stresses (orientation and magnitudes) determine a hydraulic fracture system's actual geometry. Hydraulic fracture complexity is likely to increase in shale reservoirs with high pore pressure and small differential stress [14]. In these environments, pumping pressures are high, and net pressures are small, making it easier to induce hydraulic fractures and stimulate pre-existing natural discontinuities. This markedly impacts initial production rates and ultimate recovery factors.
Several researchers developed experimental, theoretical, and numerical models to study the factors that influence a hydraulic fracture's geometric complexity in shale plays. Blanton [9,15] investigated the interaction between hydraulic fracture and pre-existing discontinuities (e.g., natural fractures). He identified horizontal differential stress and interaction angle as the most critical factors that control whether induced hydraulic fractures cross or divert from a pre-existing natural fracture. Other researchers validated this finding years later through tri-axial laboratory experiments [16][17][18][19]. These experiments also highlighted that pre-existing discontinuities could arrest hydraulic fracture propagation. Other studies also revealed that operational parameters such as fracturing fluid viscosity, pump rate, and treating pressure are governing parameters for directing hydraulic fracture propagation and final geometry [16,[20][21][22][23].
Numerical simulation is a practical method used to study hydraulic fracture geometry in shale plays. The development of 3D fracture propagation models enabled several researchers to study the interplay between natural fractures, bedding planes, and hydraulic fracture propagation [24][25][26]. These models capture the primary physical mechanisms that govern complexity of hydraulic fracture systems.
Microseismic studies and field pilots have uncovered evidence of hydraulic fracture geometries that appear to deviate appreciably from an idealized geometry and include horizontal fracture components [27][28][29][30]. Such fracture geometries result in production volumes, and rates below customary forecasts [31]. Numerous studies have investigated various hydraulic fracture geometries, but only a few integrate these geometries into reservoir simulations. The integrated hydraulic fracture geometries are often non-planar dominant fractures with varied dipping angle, length, and aperture [32][33][34]. Furthermore, simulation studies include fractures appearing to have the "tree-like" and orthogonal geometries [6,7,[35][36][37]. Such fracture networks are interpreted to result from the interactions among hydraulic fractures and pre-existing natural fractures. The presence and influence of the horizontal, mechanically weak bedding planes have not been discussed yet in these numerical simulation studies. Olorode et al. [32,35] discussed generation of horizontal fractures and assessed their effect on well productivity. They have concluded that reservoir fluid production increases as the surface area of induced fractures increases, but their study does not address the fracture closure mechanism. The studies have demonstrated that complex hydraulic fracture systems are still challenging to model and integrate into reservoir simulation analyses. Consequently, many studies lack an understanding of how these complexities impact production rates and ultimate recovery.
In this study, we conduct a parametric investigation that combines pre-existing natural fracture networks with hydraulic fracture systems. Specifically, we investigate three idealized but viable stimulation geometries (henceforth called "fracture scenarios") that combine hydraulic and pre-existing natural fracture systems as well as mechanically weak horizontal bedding planes reported to exist in shale plays. We examine their effects on wellbore flow performance using a commercial reservoir simulator (CMG-IMEX). We consider the first fracture scenario the reference case since it includes only vertical hydraulically induced fractures. In the second fracture scenario, we add stimulated vertical natural fractures ("secondary fractures") perpendicular to the hydraulic fractures. The third fracture scenario includes stimulated horizontal bedding plane fractures that intersect the vertical hydraulically induced fractures. This last scenario may occur in shale plays characterized by severe overpressure and small differential stresses or by transitional strike-slip to reverse faulting stress regimes in which the least principal stress (S 3 ) is very close to the overburden (S v ). Such conditions are reported to exist, for example, in the Marcellus shale in southern West Virginia [27], and the Tuwaiq Mountain formation in the Jafurah Basin in Saudi Arabia [38,39]. The unique reservoir simulations discussed here enable us to understand the impact of varied hydraulic fracture and stimulated natural fracture geometries and bedding planes on hydrocarbon production. Furthermore, our work emphasizes the influence of constraining present-day in situ stress state (orientation and magnitudes) and pore pressure conditions on hydrocarbon production forecasting.  We apply the dual-permeability (DK) method to model the shale play matrix and the natural fracture system. This method is standard and widely accepted to simulate unconventional shale plays [40]. The natural fracture system only exists within the stimulated reservoir volume and represents critically stressed natural fractures reactivated during hydraulic fracturing due to the increased pore pressure and associated local stress state perturbation [41]. Furthermore, we model hydraulically induced fractures (including horizontal bedding planes) with a locally spaced-logarithmically refined-dual permeability approach (LS-LR-DK) [40]. This approach is based on the DK model, which assumes an orthogonal natural fracture system with uniform fracture spacing. In our simulation model, natural fractures are spaced every 9 m (∼30 ft) in both the I-and J-directions, and their common height is equal to the reservoir thickness (i.e., 41 m or ∼135 ft). Furthermore, the stimulated reservoir volume (SRV) consists of vertical hydraulic fractures with similar properties. These fractures are distributed uniformly along the horizontal wellbore. This model generates a symmetrical reservoir with identical hydraulic fracture stages. Thus, we can simulate a single stage of the hydraulic fracturing process to limit computational time.

Numerical Simulation Model Description
In total, our model has 22 stages (Figure 2), and we multiply correspondingly the volume of reservoir fluids produced from a single stage by the number of stages [42,43]. The contribution of the exterior flow to the stimulated reservoir volume is minimal because we only consider the reactivation of natural fractures near the stimulated hydraulic fractures, and the matrix has ultra-low permeability [41]. Based on the reported massive stimulation jobs performed in mudrock plays worldwide [44][45][46][47][48][49], we estimate a total fracturing fluid volume injected into the rock of 12,700 m 3 (∼80,000 bbls) distributed evenly among the 22 stages. We consider the same fracturing fluid volume for all fracture scenarios and predetermine the stimulated volume's geometry based on each scenario. Then, the total fracturing fluid volume injected into the rock comprises the stimulated fracture volume and the leak-off volume.

Reservoir and Fluid Properties
The modeled formation is undersaturated during the studied production time. Thus, the bubble point pressure is below the reservoir pressure. The shale reservoirs in the U.S. and Saudi Arabia report pore pressure gradients between 0.014 MPa/m and 0.020 MPa/m (∼0.6 psi/ft-0.9 psi/ft) [38,44,46]. Thus, our study considers an overpressured reservoir with a reservoir fluid pressure gradient of 0.0195 MPa/m (∼0.85 psi/ft). The gas-oil and water-oil contacts are assumed to be outside the reservoir interval. The matrix compressibility is 1.0 × 10 −9 Pa −1 (6.9 × 10 −6 psi -1 ) [50][51][52]. Tables 1 and 2 summarize the shale formation and reservoir fluid properties. In Table 2, the natural fracture properties describe the critically stressed fractures. These fractures form a well-connected network that contributes to fluid flow but has low storage capacity.

Relative Permeability Curves
The two-phase flow (oil-water) relative permeability curves in this study are shown in Figure 3a,b for the rock matrix and the fractures, respectively. We assume a mixed to oil-wet target zone, which is characteristic of an organic-rich formation [57]. The organic-rich formations have a total organic carbon (TOC) up to 14 wt. % in the U.S and the Jafurah Basin-Saudi Arabia [45,46]. Experimental studies have shown that kerogen maturation increases meso-and micro-porosity in the rock matrix [58]. These pores increase the residual hydrocarbon saturation and irreducible water saturation [58]. Figure 3a shows that relative permeability values vary in a narrow window of water saturation values (0.4-0.6). It leads to sharp hydrocarbon decline rates and an extended period with low water rates. Figure 3b shows an X-shaped relative permeability curve characteristic of hydraulic fractures. This shape indicates high conductivity with no influence of capillarity.

Pressure-Dependent Permeability
Hydraulically stimulated fractures close when fracturing fluid injection stops and the pressure in the fracture drops below the fracture closure pressure [59]. The fracture closure pressure corresponds to the least principal stress (S 3 or S hmin in normal and strike-slip stress regimes) for the vertically induced hydraulic fractures, the overburden stress for horizontal bedding planes, and the maximum principal horizontal stress (S Hmax ) for the orthogonal natural fracture set ("secondary fractures"). We account for the fracture closure mechanism by implementing compaction curves to the reservoir simulation scenarios. These curves consist of pressure-dependent permeability and porosity multipliers assigned to the hydraulic fractures, bedding planes, and natural fractures grid blocks ( Figure 4). Several studies have examined the relationship between effective stress, permeability, and porosity. The studies described the relationships using empirical models (e.g., power, exponential, and logarithmic models) that consider the confining pressure constant [60,61]. These studies revealed that permeability has a higher stress sensitivity than porosity [62]. We applied the following equations to calculate the permeability and porosity multipliers [55,63]: Here, φ i : Initial porosity, fraction. -k i : Initial permeability, m 2 (mD). -P net : Net pressure, MPa (psi). It is the difference between the initial pressure and the current pressure [55]. -C frac : Coefficient of stress sensitivity, MPa −1 (psi −1 ).
(a) (b) Stimulated fracture permeability variation results in a declining conductivity that depends on the propped or unpropped nature of a respective fracture type. In our scenar-ios, the fracture conductivity reduces more slowly in a propped fracture region than in an unpropped fracture region [64]. Because secondary fractures and bedding planes have small apertures [65], we assume that proppant does not flow into these. Thus, orthogonal secondary natural fractures are self-propped due to naturally occurring roughness and horizontal bedding planes are unpropped. Furthermore, we consider that proppant settling in the vertical hydraulic fracture occurs near the wellbore, because of low viscosity fracturing fluid [66]. In this case, the distal parts of the vertical hydraulic fractures also remain unpropped. Figure 4 shows two compaction curves (red and blue lines). The red line represents the hydraulic fracture's unpropped region (Figure 5b) and the horizontal bedding planes. In contrast, the blue line describes the hydraulic fracture's propped zone (Figure 5b) and the secondary natural fracture set. Proppant support reduces the closure rate appreciably compared with the unpropped zone [55].

Reference Fracture Scenario: Vertical and Planar Hydraulic Fractures
The reference simulation model is an idealized case that discretizes the rock into a matrix and an orthogonal natural fracture network. The model comprises a hydraulic fracture stage of three vertical and planar fractures spaced 30 m (∼100 ft) apart (Figure 5a). The hydraulic fractures have an assumed aperture of 9 mm (∼0.03 ft), a propped fracture half-length of 137 m (∼450 ft), and a permeability of 2000 mD. We assume that hydraulic fracture height is equal to reservoir thickness of 41 m or ∼135 ft (i.e., no out of zone growth). In this fracture scenario, the fracturing fluid volume per fracture is 190 m 3 (∼1200 bbls). The hydraulic fracture dimensions are considered based on the literature reviewed presented in Table 2. We fixed the hydraulic fracture length and the fracturing fluid volume per fracture, assuming a proppant porosity of 0.4 (average sphere packing porosity [67]) we could estimate the hydraulic fracture aperture.
We assume vertical, bi-wing, and symmetrical hydraulic fractures and a secondary orthogonal, vertical natural fracture set (Figure 5b). This fracture scenario represents an ideal configuration that is commonly considered in hydraulic treatment designs [36,68,69]. This fracture geometry is likely to occur in hydrostatically overpressured shale plays, where the tectonic environment is characterized by low or intermediate horizontal compression (i.e., normal (S hmin < S Hmax < S v ), trans-tensional (S hmin < S Hmax ∼S v ), and strike-slip (S hmin < S v <S Hmax ) faulting stress regime) [70]. Under these conditions, hydraulic fractures are vertical and propagate in a direction perpendicular to the least principal stress (i.e., S 3 = S hmin ) [71,72]. Furthermore, they open in the direction of S 3 [73]. In this scenario, net pressures higher than the difference between the principal horizontal stresses and the rock's tensile strength combined with a high pumping rate are required to favor fracturing fluid leak-off into the secondary pre-existing natural fracture set and hydraulically open them [74]. Under these circumstances, the primary hydraulic fracture geometry is preserved [21,75]. A large stress difference between S hmin and the overburden (S v ) impedes fluid leakage from the hydraulic fractures into horizontal bedding planes [76].
The secondary natural fractures with an assumed aperture of 3 mm (∼0.01 ft) are narrower than the primary hydraulic fractures. We also assume that these secondary fractures are self-propped due to naturally occurring roughness and are filled with fracturing fluid at stimulation pressures. In contrast, hydraulic fractures have a propped and unpropped region simulated with two different compaction curves (Figure 5b).

Second Fracture Scenario: Vertical Hydraulic Fractures with Perpendicular Secondary Fractures
The second modeling fracture scenario consists of three vertical and planar hydraulic fractures and eight perpendicular secondary natural fractures ( Figure 6). The hydraulic fracture half-length is 230 m (∼750 ft), which is shorter than the fracture half-length of the reference case. The reduction in length occurs due to the initiation of secondary fractures and their interaction with primary hydraulic fractures. We assume equal primary and secondary fracture aperture (7mm or ∼0.023 ft).
In this fracture scenario, we assume a shale play under a normal faulting stress regime, where the high pore pressure is responsible for the low difference between the principal horizontal stresses. Net pressures larger than the difference between S hmin and S Hmax combined with high organic content (TOC>3 wt% [77]) favor the interaction between vertical hydraulic fractures and orthogonal pre-existing natural fractures ("secondary fractures"). Consequently, stimulation jobs hydraulically open the secondary fractures and create a stimulated interconnected orthogonal fracture network (Figure 7). Similarly to the reference case, this second fracture scenario does not consider horizontal bedding plane stimulation, because of the large difference between S v and S hmin .

Third Fracture Scenario: Vertical Hydraulic Fractures with Horizontal Bedding Plane Fractures
The third simulation configuration consists of three vertical and planar hydraulic fractures, and two stimulated horizontal bedding planes (Figure 8). The primary fracture half-length is 230 m (∼750 ft), and they are spaced uniformly every 30 m (∼100 ft) along the horizontal wellbore. Furthermore, we assume equal apertures for both primary hydraulic fractures and bedding planes (4 mm or ∼0.013 ft). The size of the stimulation job is the same as the other two simulation fracture scenarios presented before.  In this case, stimulation of horizontal bedding planes increases the hydraulic fracture configuration complexity (Figure 9). Horizontal bedding planes, mechanically weak owing to little cohesion, may be stimulated in shale plays in transitional strike-slip to reverse faulting stress regimes, in which S v and S hmin are close in magnitudes and similar to the least principal stress (i.e., S hmin ∼S v ∼S 3 ). Such stress conditions have been reported in some unconventional source rock plays (e.g., the Marcellus Shale in southern West Virginia and the Tuwaiq Mountain formation in the Jafurah Basin, Saudi Arabia [27,38,39]). These environments are characterized by high pore pressure, which highly reduces the effective stress magnitudes, and makes them close to each other [13]. Under these conditions, hydraulic fracturing may stimulate horizontal bedding planes when hydraulic communication occurs with the primary vertical hydraulic fractures. Horizontal bedding planes open to flow at high pumping rates and net pressures [22,76,79]. However, we assume that the aperture is not sufficiently large for proppants to enter, resulting in their closure during fracturing fluid flowback.
This hydraulically stimulated configuration results in simultaneous fracture growth in the vertical and horizontal directions. The fracturing fluid propagates more easily along laminations than across them due to high horizontal permeability in shale plays [11]. Despite this process, hydraulic fractures may continue to cross the bedding planes and promote their vertical propagation when the fracture's pressure is sufficiently high [80].

Results and Discussion
Operators have applied large stimulation jobs in shale plays to achieve production at economical rates. In 2020, the expected mean value (P50) of fracturing fluid injected into U.S. shale plays, was about 60,600 m 3 (∼380,000 bbls) for an average wellbore lengths of 2500 m (∼8200 ft). The mean proppant mass injected was 7500 tons (∼15 MM lbs) [48,49]. Similarly, hydraulic treatments in the Jafurah Basin in Saudi Arabia require large fracturing fluid volumes and proppant amounts. From published data, we deduce that total fracturing fluid volumes and proppant mass per wellbore are about 18,000 m 3 (∼113,000 bbls) and 3600 tons (∼7.2 MM lbs) [44,45,47,81]. The average horizontal wellbore length in the Jafurah basin is approximately 1500 m (∼5000 ft) [45,53]. Some mudrock plays are further characterized by high pore pressure and a small difference between the principal stresses. For example, the Jafurah basin exhibits pore pressure gradients between 0.016 and 0.020 MPa/m (∼0.7-0.88 psi/ft) [38,44], and strike-slip faulting stress regime (S hmin ∼0.96 psi/ft < S v ∼ 1.1 psi/ft < S Hmax ∼1.2 psi/ft) [38,44,82]. These environments combined with large fluid volume injection foster interactions between vertical hydraulic fractures and pre-existing discontinuities (e.g., natural fractures and bedding planes). Based on Haider et al. studies [83], 90% of the fracturing fluid injected into shale systems can be retained in the rock discontinuities and may result in hydrocarbon production below the forecast. Therefore, we have studied the impact of rock discontinuities during hydraulic stimulation on wellbore flow performance in the overpressured shale plays. We have proposed three different hydraulic fracture geometries ("fracture scenarios") and compared two simulation cases with the variable primary fracture conductivity for the scenario with horizontal fractures (fracture scenario #3). Figure 10 compares reservoir fluid production for the three fracture scenarios, assuming equal primary hydraulic fracture permeability of 2000 mD. This case assumes an efficient proppant placement and high conductivity in the vertical hydraulic fractures. The results show that cumulative fluid production slightly differs between the fracture scenarios #2 and #3. Table 3 compares the different parameters of the three proposed fracture scenarios. The reference fracture scenario produces an ultimate oil production of 407 Mbbl and ultimate water production of 564 Mbbl at 15 years. For this scenario, the oil production peaks at 931 bbl/d, and the water production peaks at 18,992 bbl/d. The ultimate recovery factor at 15 years is 8.5%. The table shows that the stimulation of secondary vertical, orthogonal natural fractures increases the initial water production by 29% and the initial oil production by 10%. These fractures are the highly conductive pathways that enhance reservoir production. Despite this, the orthogonal natural fractures close after a few months of production because of ineffective self-propping. Thus, such a configuration reduces cumulative reservoir fluid production after 15 years by an average of 10% compared to the reference fracture scenario. After 15 years, our model predicts an ultimate oil recovery factor of about 7.7%.
On the other hand, stimulated horizontal bedding planes reduce the initial water production by 40% and the initial oil production by 13%. In this case, the unpropped planes experience a rapid closure once pumping stops and fracturing fluid flowback starts. This closure disconnects a huge stimulated reservoir volume, which has a marked impact on well productivity. For this fracture scenario, we forecast a similar recovery factor as scenario #2 because we assume the same primary fracture length and conductivity in both scenarios.  The interaction between hydraulic fractures and bedding planes reduces the hydraulic fracture propagation energy. Thus, it creates narrower and shorter vertical hydraulic fractures [72], and generates poor proppant placement. For this reason, in the second simulation case (Figure 11), we reduce the vertical hydraulic fracture permeability in fracture scenario #3 by 90%. In this simulation case, fracture scenario #3 agrees well with reported wellbore flow performance in shale plays with high pore pressure and small differential stress. Figure 11 compares reservoir fluid production for the three fracture scenarios. Fracture scenarios #1 and #2 are the same as presented in the first simulation case. Despite this, the results show that horizontal fracture generation significantly reduces total fluid production. Ultimate oil recovery for fracture scenario #3 after 15 years is 6.7%. Table 4 compares the different parameters of the three proposed fracture scenarios in the second simulation case. The results show that stimulating horizontal bedding planes reduces the initial water production by about 80% and initial oil production by about 50% compared to the reference fracture scenario. The cumulative oil production after 15 years is reduced by 20% resulting from rapid fracture closure and low hydraulic fracture conductivity. Table 4. Simulation case 2: production comparison after 15 years between fracture scenarios #2 and #3, and the reference fracture scenario #1. The reference scenario comprises 3 vertical hydraulic fractures (HF), fracture scenario #2 has 3 vertical hydraulic fractures and 8 orthogonal natural fractures (NF), and fracture scenario #3 has 3 vertical hydraulic fractures and 2 horizontal bedding planes (BP). In this case, vertical hydraulic fractures' permeability is assumed to be the 2000 mD for scenarios #1 and #2. Scenario #3 has a vertical hydraulic fracture permeability of 200 mD.

Parameter
Reference

Conclusions
• The interactions among vertical hydraulic fractures and the orthogonal pre-existing natural fractures may initiate and hydraulically open secondary natural fractures.
Our results show that these secondary fractures improve reservoir permeability and contribute to the initial hydrocarbon production when they reach conductivity similar to that of primary vertical hydraulic fractures. The presence of stimulated vertical orthogonal natural fractures enhances the initial oil production by about 10%. • We allow the stimulation of secondary orthogonal natural fractures to reduce the vertical hydraulic fracture length and aperture. The shortened primary fracture length combined with the rapid closure of secondary natural fractures, reduce the stimulated reservoir volume and harms long-term production. Thus, our results show that these two factors can reduce the cumulative oil production at 15 years by about 8% compared with the scenario with vertical hydraulic fractures only. • We assume that stimulation of the horizontal bedding planes reduces the primary vertical hydraulic fracture length and aperture. This may produce premature proppant screen-out and inefficient proppant placement into primary hydraulic fractures. Our results show that the lowered vertical hydraulic fracture conductivity combined with rapid horizontal fracture closure can reduce cumulative hydrocarbon production after 15 years by about 20% and the initial hydrocarbon production by about 50% compared with the vertical hydraulic fractures only (" he reference fracture scenario"). • In hydraulic fracturing, it is crucial to constrain the present-day in situ stress state (magnitudes and orientation) and pore pressure conditions to understand the initiation and propagation mechanisms that control the final stimulated reservoir volume geometry. Only then realistic hydrocarbon production forecasts can be obtained. The assumption of vertical, symmetrical, and planar fractures may lead to unreliable hydrocarbon production forecasts. This last finding agrees with the studies performed by previous researchers [6,35,36,69]. The simplification of hydraulic fracture geometry can overestimate hydrocarbon production, primarily if the modeled formation is characterized by high pore pressure and a small difference between the principal stresses. These conditions enhance the fracturing fluid leakage into natural fractures and/or bedding planes. Acknowledgments: The authors thank the Ali I. Al-Naimi Petroleum Engineering Research Center (ANPERC) at King Abdullah University of Science and Technology (KAUST) for supporting this research. We also thank the financial support from KAUST through Professor Tadeusz Patzek baseline. Finally, we thank all the reviewers for their detailed, valuable, and timely suggestions that helped us improve our article's quality.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: