Mathematical modeling of geomechanical behavior of tarmat during the depletion of giant oil reservoir-aquifer systems

Abstract

In this work, deformation and failure behavior of tarmat layers during depletion of a giant reservoir–aquifer system has been studied. Deformation response of the tarmat to increasing pressure differential caused by continuous depletion of reservoir is examined and a mathematical model is developed for the study of this type of composite systems. The geomechanical failure that takes place when the pressure differential reaches a critical value is also evaluated, along with the characterization of the resulting fracture. Plate theory, maximum shear stress failure criterion, conventional well test model, Perkins–Kern–Nordgren (PKN) and Khristianovic–Geertsma–de Klerk (KGD) models and flow through fractures models are used. The developed sensitivity analysis proposes the proper protocol to be followed in order to undertake production design in such composite systems. The methodology presented in this paper, ultimately, predicts fracture width and fracture permeability that would be developed in a system with a tarmat layer having a certain thickness and a reservoir being produced at a certain production rate and total depletion time.

Introduction

Oil resources are located in various types of reservoir formations, varying with properties, dimensions and architectures. Creating feasible production strategies with a reasonable exploration and development plan is of great importance in the production of these oil sources. The process requires a good understanding of oil, reservoir properties and existing geological architecture of the reservoir of interest. Giant oil fields, being oil sources with high production potentials, contain more than 500 million barrels of recoverable oil as they constitute almost 75% of the recoverable oil resources in the world.

This study focuses on a three-layered composite system, typical of giant oil fields in the Middle East. Upper layer contains mobile reservoir fluids, middle layer is referred as the tarmat, and the bottom layer is a high pressure water aquifer. Tarmat is an extremely viscous hydrocarbon layer, mainly composed of tar or bitumen, which exists between oil and water contact. In many cases, the tarmat acts as a permeability barrier between the reservoir and its aquifer. This composite system is depicted in Fig. 1.

Fig. 1

Schematic representation of the system under consideration

The main objectives of this study are: (1) the characterization of the geomechanical behavior and eventual failure of the tarmat layer as a response to hydrocarbon production and the associated significant increase in pressure differential between the depleting reservoir and its aquifer; and (2) the evaluation of the system behavior after geomechanical failure of the tarmat takes place. The latter part of the analysis involves the characterization of the fracture and its permeability and the resulting communication between the aquifer/reservoir system.

Tarmat deformation analysis

By recognizing the tarmat layer as a flat plate with a thickness significantly smaller than its areal dimensions, plate theory can be used for the analysis of tarmat deformation. Plate theory extends the findings of the theory of beams for these types of structural elements (Boresi and Schmidt 2003; Bickford 1998; Timoshenko and Woinowski-Krieger 1959). The tarmat plate is assumed to be simply supported having rectangular shapes and lateral dimensions that are perpendicular to x and y axes. Deformations can take place in both x and y directions. Tarmat response under uniform and non-uniform loading has been studied, and the resulting forces are a consequence of the developed normal and shear stresses. Figure 2 shows a deformed plate with the resulting forces, the way they act and their relations with each other. These forces include bending and twisting moments (M’s), and shear forces (V), which occur throughout the plane due to the load (q).

Fig. 2

Bending moments, shear forces on deformed plate, force resultants acting on the plate element

For this system, a balance of forces in x and y directions yield the following equilibrium equation relating the resulting moments generated by the applied load (Boresi and Schmidt 2003):

$$ \frac{{\partial^{2} M_{x} }}{{\partial x^{2} }} + \frac{{\partial^{2} M_{y} }}{{\partial y^{2} }} + \frac{{\partial^{2} M_{xy} }}{\partial y\partial x} + \frac{{\partial^{2} M_{yx} }}{\partial x\partial y} = q $$

(1)

where q = load applied to the system. A plate subjected to transverse loading with certain distribution is displaced perpendicular to its middle plane. Strain–displacement relations can be derived following the definitions of normal and shear stresses in terms of plane displacement or deformation (w), as follows:

$$ \begin{gathered} \sigma_{x} = - \frac{Ez}{{1 - \nu^{2} }}\left( {\frac{{\partial^{2} w}}{{\partial x^{2} }} + \nu \frac{{\partial^{2} w}}{{\partial y^{2} }}} \right),\quad \sigma_{y} = - \frac{Ez}{{1 - \nu^{2} }}\left( {\frac{{\partial^{2} w}}{{\partial y^{2} }} + \nu \frac{{\partial^{2} w}}{{\partial x^{2} }}} \right) \hfill \\ \tau_{xy} = - \frac{Ez}{{1 - \nu^{2} }}\left( {1 - \nu } \right)\left( {\frac{{\partial^{2} w}}{\partial x\partial y}} \right) \hfill \\ \end{gathered} $$

(2)

These expressions, when substituted into definitions of bending and twisting moments, yield:

$$ \begin{gathered} M_{x} = - \int\limits_{{ - {h \mathord{\left/ {\vphantom {h 2}} \right. \kern-\nulldelimiterspace} 2}}}^{{{h \mathord{\left/ {\vphantom {h 2}} \right. \kern-\nulldelimiterspace} 2}}} {z\sigma_{x} {\text{d}}z = D\left( {\frac{{\partial^{2} w}}{{\partial x^{2} }} + \nu \frac{{\partial^{2} w}}{{\partial y^{2} }}} \right)} \hfill \\ M_{y} = - \int\limits_{{ - {h \mathord{\left/ {\vphantom {h 2}} \right. \kern-\nulldelimiterspace} 2}}}^{{{h \mathord{\left/ {\vphantom {h 2}} \right. \kern-\nulldelimiterspace} 2}}} {z\sigma_{y} {\text{d}}z = D\left( {\frac{{\partial^{2} w}}{{\partial y^{2} }} + \nu \frac{{\partial^{2} w}}{{\partial x^{2} }}} \right)} ,\quad D = \frac{{Eh^{3} }}{{12\left( {1 - \nu^{2} } \right)}} \hfill \\ M_{xy} = - M_{yx} - \int\limits_{{ - {h \mathord{\left/ {\vphantom {h 2}} \right. \kern-\nulldelimiterspace} 2}}}^{{{h \mathord{\left/ {\vphantom {h 2}} \right. \kern-\nulldelimiterspace} 2}}} {z\tau_{xy} {\text{d}}z = D\left( {1 - \nu } \right)} \frac{{\partial^{2} w}}{\partial x\partial y} \hfill \\ \end{gathered} $$

(3)

where D is referred as flexural rigidity. Expressions in Eq. 3 can be substituted in Eq. 1 to allow the derivation of the biharmonic equation, shown below as Eq. 4:

$$ D\left( {\frac{{\partial^{4} w}}{{\partial x^{4} }} + 2\frac{{\partial^{4} w}}{{\partial x^{2} \partial y^{2} }} + \frac{{\partial^{4} w}}{{\partial y^{4} }}} \right) = q,\quad D\nabla^{2} \left( {\nabla^{2} w} \right) = q $$

(4)

In this study, all of the edges of the plate are considered to have simply supported boundary conditions. There are two main conditions to be satisfied by simply supported edges. First, the displacement (w) must be equal to zero at the edges and any moment that coincides in direction with the direction of the edge must be equal to zero. Therefore, the following boundary conditions are to be satisfied by the y = constant and x = constant edges:

$$ \begin{gathered} y{\text{ - constant:}}\;w = 0,\quad M_{y} = D\left( {\frac{{\partial^{2} w}}{{\partial y^{2} }} + \nu \frac{{\partial^{2} w}}{{\partial x^{2} }}} \right) = 0 \hfill \\ x{\text{ - constant:}}\;w = 0,\quad M_{x} = D\left( {\frac{{\partial^{2} w}}{{\partial x^{2} }} + \nu \frac{{\partial^{2} w}}{{\partial y^{2} }}} \right) = 0 \hfill \\ \end{gathered} $$

(5)

The biharmonic equation, the fourth degree differential equation in Eq. 4, can thus be solved with these four boundary conditions and a relevant loading expression, to obtain the expression for transverse deformation w(x,y). For the case of a uniform loading, for example, the following solution can be found for plate transverse deformation:

$$ w(x,y) = \frac{{16q_{0} }}{{\pi^{6} D}}\sum\limits_{{m,{\text{odd}}}} {\sum\limits_{{n,{\text{odd}}}} {\frac{{\sin \frac{m\pi x}{a}\sin \frac{n\pi y}{b}}}{{mn\left( {\left( \frac{m}{a} \right)^{2} + \left( \frac{n}{b} \right)^{2} } \right)^{2} }}} } $$

(6)

Once the transverse deformation is found, stresses and bending and twisting moments and twisting moments across the plate can be calculated using Eqs. 2 and 3.

Tarmat failure analysis

In order to evaluate the tarmat behavior at the moment of geomechanical failure, tarmat deformation analysis should be coupled with a failure criterion. There are a number of theories that predict failure as a function of prevailing stresses. For example, the maximum shear stress failure theory, as originally developed by Charles Coulomb and Henry Tresca, indicates that the failure point is reached when the maximum shear stress in the material becomes equal to the value of the shear stress at yielding. This point, which indicates the occurance of failure is referred as yield strength. Yield strenght is a property of the material which needs to be experimentally determined by uniaxial compression or uniaxial tension test. The Mohr–Coulomb failure (internal friction) criterion is another common way of evaluating failure by relating shearing resistance to contact forces, friction and cohesion that is present among the rock grains and it is deemed to be appropriate for the prediction of failure in brittle materials. Other failure criteria include the maximum normal stress failure criterion, Hoek–Brown failure criterion, Von Mises failure criterion, and the octahedral shear stress theory, among others.

If the maximum shear stress failure criterion is utilized, principle stresses must be calculated for the material. They can be estimated by implementing Eq. 7 below:

$$ \sigma_{1,2} = \pm \frac{{\sigma_{x} + \sigma_{y} }}{2} \pm \sqrt {\frac{{\sigma_{x} - \sigma_{y} }}{2} + \tau_{xy}^{2} } $$

(7)

In this equation, maximum values for the stresses must be used. In a simply supported plate, these occur at the center of the structure. At the center of the plate, maximum stresses are given by the expressions:

$$ \left( {\sigma_{x} } \right)_{\max } = \pm \frac{{6M_{x} }}{{h^{2} }},\;\quad \left( {\sigma_{y} } \right)_{\max } = \pm \frac{{6M_{y} }}{{h^{2} }},\quad \;\left( {\tau_{xy} } \right)_{\max } = \pm \frac{{6M_{xy} }}{{h^{2} }} $$

(8)

It can be shown that for the case of interest,

$$ \sigma_{2} \left\langle {\sigma_{1} } \right.\left\langle 0 \right.,\;\sigma_{3} = 0 $$

(9)

Therefore, if $\left| {\sigma_{2} } \right|$ exceeds $\sigma_{\text{YS}}$ of the material, failure occurs (Gere 2001; Bickford 1998).

Fracture width and fracture permeability analysis

For the purpose of fracture width analysis, two different hydraulic fracturing models have been used: the Perkins–Kern–Nordgren (PKN) and the Khristianovich–Zheltov–Geertsma–deKlerk (KGD) models. They relate fracture width to properties of fluid and rock in the fractured system. The working fluid in our model is the aquifer water and the rock is represented by the tarmat. In both models, the most influential fluid properties are viscosity and specific gravity and influential rock properties are Poisson’s ratio, Young’s modulus of elasticty. Additional variables of importance are fluid injection rate and fracture length. Fracture thickness is a property that is not influential in PKN model but is influential in KGD model. Inputs into the analysis protocol used in this study are viscosity and specific gravity of water, Poisson’s ratio and Young’s Modulus of Elasticiy of tarmat, reservoir production rate and thickness of tarmat. Reservoir production rate is assumed close to fluid injection/breakthrough rate since both are expected to create a similar pressure differential in magnitude and direction.

The problem is thus approached as an inverse problem, in which fracture width is estimated as a function of reservoir production rate. In this problem, thickness of tarmat represents the fracture penetration length. In the hydraulic fracturing analog, fracture length is the parameter that helps to express the penetration of the crack created while in the case of our interest, the target crack penetration is the tarmat thickness. Fracture thickness in the KGD model is assumed to correspond to this fracture length. Both the PKN and KGD models are shown to be in agreement when the fracture thickness is assumed to be as long as the line that is drawn at the points on the plate where shear stress is 99% of the maximum shear stress.

The expressions for fracture width calculation for the PKN and KGD models are presented below in Eq. 10, respectively:

$$ \begin{gathered} w = 0.3\left[ {\frac{{q\mu \left( {1 - v} \right)\left( {{{h_{0} } \mathord{\left/ {\vphantom {{h_{0} } 2}} \right. \kern-\nulldelimiterspace} 2}} \right)}}{G}} \right]^{{{1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-\nulldelimiterspace} 4}}} \left( {\frac{\pi }{4}\gamma } \right),\quad G = \frac{E}{{2\left( {1 + v} \right)}} \hfill \\ w = 0.29\left[ {\frac{{q\mu \left( {1 - v} \right)\left( {{{h_{0} } \mathord{\left/ {\vphantom {{h_{0} } 2}} \right. \kern-\nulldelimiterspace} 2}} \right)^{2} }}{{Gh_{f} }}} \right]^{{{1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-\nulldelimiterspace} 4}}} \left( {\frac{\pi }{4}} \right) \hfill \\ \end{gathered} $$

(10)

In these expressions, G is the shear modulus. Figure 3 explains fracture thickness and fracture length orientation for a hydraulic fracturing case and the problem studied here.

Fig. 3

Fracture length, fracture width and fracture thickness in a hydraulic crack

Once fracture width is estimated, the associated fracture permeability can be calculated. For example, fracture permeability can be related to fracture width by equaling the Poiseuille’s Law in parallel plates to Darcy’s law in porous media (Craft and Hawkins 1959) which yields the expression:

$$ q = \frac{{\left( {P_{1} - P_{2} } \right)A^{\prime}w^{2} }}{{12\mu h_{0} }},\quad k = 54 \times 10^{6} w^{2} $$

(11)

A similar expression can be independently derived by considering the case of the flow of hydraulic fracturing fluids through induced fractures (Yew 1997). Yew (1997) assumed that fractures have a narrow opening of constant width all through the fracture thickness. If the flowing fluid is assumed to be an incompressible Newtonian fluid, Yew (1997) shows that fracture width and permeability can be related by Eq. 12:

$$ \frac{{w^{2} }}{12} = k,\quad k = 5.45 \times 10^{7} w^{2} $$

(12)

Reservoir depletion and pressure transient model

At this stage of the analysis, reservoir depletion must be considered in order to recreate the magnitude of the load placed on the tarmat place. In order to relate pressure depletion with time evolution, the standard computational procedure of classical well test model is followed (Earlougher 1977; Lee 1982). Reservoir is single phase and square-shaped. Well is assumed to be located at the center of the drainage area. In this procedure, dimensionless time is converted to actual time and dimensionless pressure is converted into actual pressure. Dimensionless production rate is used as an intermediate step in these calculations. The solutions for dimensionless pressure drop (p_d) versus dimensionless time (t_d) tabulated by Earlougher (1977) for the case pressure variation at the center of the rectangular-shape reservoirs are utilized. Equations 13, 14 and 15 give the expressions for dimensionless production rate, dimensionless pressure and dimensionless time, respectively:

$$ q_{\text{D}} = \frac{{\gamma Bq_{p} \mu }}{{khP_{i} }} $$

(13)

$$ \Updelta P_{\text{D}} = \frac{{P_{i} - P}}{{P_{i} q_{\text{D}} }} $$

(14)

$$ t_{\text{D}} = \frac{\lambda kt}{{\phi \mu cr_{w}^{2} }},\quad t_{\text{DA}} = t_{\text{D}} \frac{{r_{w}^{2} }}{A},\quad t_{\text{DA}} = \frac{\lambda kt}{\phi \mu cA} $$

(15)

Results and discussions

Tarmat failure analysis

Figure 4 shows the expected deformation response as a function of an uniformly applied loading along with the associated failure envelope for two reservoir scenarios with different properties. In these figures, deformation versus loading behavior is investigated until the point where maximum shear stress failure criterion predicts geomechanical failure. The rectangular plate deformation model is used to model the behavior of the tarmat. In these figures, it is readily seen that thicker the tarmat, the larger the pressure drop required to trigger geomechanical failure. At the same time, the thinner the tarmat, the larger the deformation experienced by the tarmat prior to the onset of failure. Inspection of Fig. 4 reveals that as the Young’s modulus of elasticity and yield strength becomes larger, the pressure differential required for the tarmat to fail increases while its deformation is expected to slightly decrease.

Fig. 4

Deformation versus loading with associated failure envelope. aE = 3,000,000 psi, σ_YS = 30,000 psi, a = b = 750 ft, v = 0.30, bE = 5,500,000 psi, σ_YS = 50,000 psi, a = b = 750 ft, v = 0.30

Figure 5 displays the magnitude and nature of the effects caused by the principal parameters of this analysis on failure pressure including lateral dimensions, yield strength and Poisson’s ratio of the tarmat. Each of the figures is drawn considering critical pressures and deformations that occur until critical pressure is reached. It is observed that in smaller drainage areas, more pressure differential is required to fail tarmat, tarmat having higher yield strength, requires more pressure differential until failure and tarmat with smaller Poisson’s ratios, requires more pressure differential until the failure point. All of these parameters indicate a direct proportionality with tarmat thickness and magnitude of pressure required to fail the tarmat.

Fig. 5

a Loading versus thickness for various lateral dimensions (E = 4,000,000 psi, σ_YS = 40,000 psi, v = 0.30). b Loading versus thickness for various yield strengths (E = 4,000,000 psi, v = 0.30, a = b = 750 ft). c Loading versus thickness for various Poisson’s ratios (E = 4,000,000 psi, σ_YS = 40,000 psi, a = b = 750 ft)

Figure 6 displays a comparison between deformation versus loading and associated failure envelope of uniform loading (a) and non-uniform loading (b). For a tarmat with a certain thickness, lateral dimensions, Young’s modulus of elasticity, Poisson’s ratio and yield strength, more pressure differential is required to observe failure in the case of non-uniform loading. A comparison of Fig. 6a and b reveals a slightly more deformation in the case of a uniform loading.

Fig. 6

Comparison of deformation versus loading with associated failure envelope with a uniform, and b non-uniform loading assumption (E = 8,000,000 psi, σ_YS = 75,000 psi, a = b = 750 ft, v = 0.30)

Figure 7 displays a comparison between loading versus thickness for different lateral dimensions of tarmat in the cases of uniform loading (a) and non-uniform loading (b); respectively. A similar comparison can be made for yield strength and Poisson’s Ratio effects. For a tarmat of certain influential properties more pressure differential is required to observe geomechanical failure in the case of non-uniform loading.

Fig. 7

Comparison of loading versus thickness graphs for various lateral dimensions with uniform and non-uniform loading assumption

Reservoir depletion model

This study has been conducted for two different drainage area assumptions. In each case, production rate has been varied between 1,000 and 10,000 STB/day. Reservoir properties assigned in this analysis are given in Table 1 below.

Table 1 Hydrocarbon reservoir properties

Figure 8 provides a comparison of different drainage area assumptions while production rate influences on pressure and time relationship can also be observed. Figure 8a and b represent the analysis with drainage area assumption of 51.65 and 200 acres. As production rate increases, a certain pressure differential is reached in a shorter period of time. For reservoirs with smaller drainage areas, it takes less time to reach a certain pressure differential than it does for those with larger drainage areas.

Fig. 8

Pressure differential versus time differential graphs for various flow rates. a Cross-sectional area = 51.65 acres, b cross-sectional area = 200 acres

Fracture permeability characterization

In PKN and KGD models, which are used to predict fracture width, production rate from the reservoir above that would make a similar effect as injection rate from the aquifer below has been used as an input flow rate. The production rates that have been used in the analysis with conventional well test model is used in this part of the study. The production rates are changed between 1,000 to 10,000 STB/day. The same analysis that is relating production rate and width has been repeated for a possible range of tarmat thicknesses varying between 30 and 100 ft. Figure 9 includes the relationship between fracture width and production rate for PKN and KGD models, for various properties and relationship between fracture permeability and fracture width. It should be noted that Eq. 1 refers to the first method and Eq. 2 refers to the second method.

Fig. 9

Width versus production rate graphs for various tarmat thicknesses for PKN and KGD models. aE = 3,000,000 psi, v = 0.30, h_f = 47.82 ft, bE = 5,500,000 psi, v = 0.30, h_f = 47.82 ft, cE = 3,000,000 psi, v = 0.30, h_f = 47.82 ft, d fracture permeability versus fracture width

In both models, Young’s modulus of elasticity is observed to be inversely proportional with fracture width. A system with known properties encounters a wider fracture width if the reservoir at the top of the system is being produced with a higher production rate. Wider fracture widths would be created in systems with thicker tarmat layers. PKN model predicts larger widths than KGD does. Difference in this prediction is largest in systems with high reservoir production rates and thick tarmat layers.

Coupled analysis protocol

This part of the section outlines the methodology followed in this study coupling each individual step. Suggested protocol is explained through the use of a composite system. Properties of each layer of the composite system are given in Table 2.

Table 2 Properties of composite system

First step involves construction of deformation versus loading with associated failure envelope. This relation is dependent on Young’s modulus of elasticity, yield strength, Poisson’s ratio and lateral dimensions of the tarmat. Figure 10a represents the outcome of the study on this relationship. By entering the chart at the 80 ft tarmat thickness line, magnitudes of deformation and pressure are found to be 15 ft and 432 psi, respectively, at the time of failure. Second step involves the use of conventional well test model to obtain the relationship between pressure differential and time, for various flow rates within the range chosen. Figure 10b represents the results of this analysis (data used in this analysis are from the hydrocarbon reservoir part of Table 2). Two different times chosen are 3 and 7 days. In the analysis, the failure pressure of 432 psi from the first stage is used with the production rates for 3 and 7 days as 5,000 and 2,000 STB/day, respectively. This is followed with the fracture width determination and permeability analysis using the PKN and KGD models. A selected range of production rates and tarmat thicknesses are considered. Relation between fracture width and production rate is the output of Figure 11 is the output (data is presented in tarmat and fluid sections of Table 2). The 80 ft thickness from family of thickness curves in PKN and KGD models is selected. Fracture widths for 2,000 STB/day are predicted to be 0.0140 and 0.0130 in. Fracture widths for 5,000 STB/day are predicted to be 0.0180 and 0.0165 in; respectively. Final step of the analysis is determination of the permeability as related to the fracture widths. Average values of fracture widths of 0.0135 and 0.0173 in yield fracture permeabilities of 10,500 and 16,000 darcy, respectively (Fig. 12).

Fig. 10

a Deformation versus loading with associated failure envelope (E = 8,000,000 psi, σ_YS = 30,000 psi, a = b = 1,500 ft, v = 0.30). b Pressure differential versus time differential graphs for various flow rates (cross-sectional area = 51.65 acres)

Fig. 11

Width versus production rate graphs for various tarmat thicknesses (PKN and KGD model) (E = 8,000,000 psi, v = 0.30, h_f = 47.82 ft)

Fig. 12

Fracture permeability versus fracture width

Summary and conclusions

In this study, we have examined the tarmat deformation and failure behavior of a giant oil reservoir-aquifer system undergoing a depletion process. The fracture that develops after failure is characterized and its permeability is determined. In the analysis, plate theory, maximum shear stress failure criterion, classical well test analysis theory, PKN model, KGD model, and fracture flow models are used. A sensitivity analysis is conducted by varying the reservoir, rock and fluid properties. The proposed methodology, predicts fracture width and fracture permeability that would be created in a system with a tarmat layer of certain thickness as a function of rate of depletion and total depletion time. In the proposed analysis protocol, each analysis step focuses on the behavior of a layer individually. As an alternative approach the composite system in its entirety can be analyzed by coupling of each layer. This will allow computation of the fracture penetration rate through the tarmat layer as a function of time. Also, as a continuation of the work presented in this paper, we would like to take the proposed solution one step further by integrating the formulated geomechanical model with the fluid flow (and/or heat flow) models so that all of the unknowns are solved simultaneously.

Within the bounds of the analysis protocol presented in this paper, following conclusions are drawn:

• 1. As thickness of tarmat increases, total deformation that occurs until the failure point decreases, and pressure differential that is required to fail the tarmat increases.

• 2. As Young’s modulus of elasticity and yield strength of tarmat increase, pressure differential that is required to fail tarmat increases, and magnitude of deformation that occurs until failure of tarmat decreases.

• 3. As lateral dimensions of tarmat increase, pressure differential that is required to fail the tarmat decreases.

• 4. As yield strength of tarmat increases, pressure differential that is required to fail the tarmat increases.

• 5. As Poisson’s ratio of tarmat increases, pressure differential that is required to fail the tarmat decreases.

• 6. A case of non-uniform loading requires more pressure and the system experiences less deformation until failure as compared to a similar case with uniform loading configuration.

• 7. The PKN model predicts larger fracture widths than the KGD model does. The difference in predictions becomes more obvious in the presence of thick tarmat layers and high production rates.

• 8. The PKN and the KGD models predict wider fracture widths for higher reservoir production rates and thicker tarmat layers.