Impact of Flowback Dynamics on Fracture Conductivity

Journal Pre-proof Impact of flowback dynamics on fracture conductivity A.A. Osiptsov, I.A. Garagash, S.A. Boronin, K.I. Tolmacheva, K.E. Lezhnev, G.V. Paderin PII:

S0920-4105(19)31241-0

DOI:

https://doi.org/10.1016/j.petrol.2019.106822

Reference:

PETROL 106822

To appear in:

Journal of Petroleum Science and Engineering

Received Date: 28 June 2019 Revised Date:

28 September 2019

Accepted Date: 16 December 2019

Please cite this article as: Osiptsov, A.A., Garagash, I.A., Boronin, S.A., Tolmacheva, K.I., Lezhnev, K.E., Paderin, G.V., Impact of flowback dynamics on fracture conductivity, Journal of Petroleum Science and Engineering (2020), doi: https://doi.org/10.1016/j.petrol.2019.106822. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

Impact of Flowback Dynamics on Fracture Conductivity A.A. Osiptsova , I.A. Garagasha,c , S.A. Boronina , K.I. Tolmachevaa , K.E. Lezhnevb , G.V. Paderinb a Multiphase

Systems Lab, Skolkovo Institute of Science and Technology (Skoltech), Russia b LLC Gazpromneft Science and Technology Center, Russia c The Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences, Moscow, Russia

Abstract We aim at the development of a general modeling workflow for design and optimization of the flowback operation on a multistage fractured well. In our methodology, we employ a coupled modeling approach, where geomechanics phenomena are considered together in the framework of a global model for fluids displacement in a fracture. Among the solid mechanics effects considered are: proppant embedment with plastic deformations and creep of the rock near the fracture surface, proppant crushing, proppant flowback, proppant pack compaction under the action of the closure stress with elastic compaction of proppant grains and frictional rearrangement of the grains, tensile rock failure at high filtration rates from the reservoir to the fracture (when drawdown pressure drop is excessively high). These geomechanics effects are implemented into the global fluid flow model in a propped fracture, which takes into account: displacement of fracturing fluid by oil, two-phase effects, suspension filtration with particle trapping and mobilization, resulting in permeability damage, colmatation and skin buildup inside the fracture in the near-wellbore zone (where particles are mainly formation fines and crushed proppant). In our numerical experiments and parametric study, we made the following observations, structured by the discipline. Geomechanics: proppant pack compaction and proppant embedment due to rock plasticity and creep may reduce the fracture width by tens percent from the original (ideal) estimates, thus having a profound impact on fracture conductivity. Tensile rock failure in the near-wellbore zone may result in fracture being disconnected hydraulically from the wellbore. Fluid mechanics: colmatation of the proppant pack by fines (formation fines or crushed proppant) may result in further decrease in the conductivity. We developed a modeling workflow for the transient flowback from multistage-fractured wells. Based on the results of parametric studies, we provide evidence that a ”smooth” scenario of piece-wise constant choke opening is preferable in terms of the preserved nondimensional fracture conductivity and increased cumulative well production, as compared to an ”aggressive” scenario of rapidly opening the choke. Implications for field operations are discussed along with modeling-based scenarios to be evaluated in a field-testing campaign. The novelty of our approach stems from the ambition to construct a general framework for a coupled geomechanics-fluid mechanics model of flowback, accurately equipped with proper submodels for each important physics phenomenon peculiar to the flowback operation. Keywords: fracturing, flowback, tensile rock failure, formation damage, colmatation, fracture conductivity, proppant embedment, proppant pack compaction, proppant crushing, sanding, formation fines

1. Introduction Hydraulic fracturing is a technology for stimulation of oil and gas production based on injection of particle-fluid suspension into a well for creating fractures in the hydrocarbonbearing rocks to provide highly conductive channels (see Fig. 1). The process of fracture creation and propagation and proppant placement is carefully designed with the help of simulators of hydraulic fracturing based on advanced continuum mechanics models [1]. The subsequent long term production is also thoughtfully planned with the help of 3D reservoir simulators. At the same time there is relatively short operation (on the time scale of long-term production), called flowback, which importance has long been underestimated. By flowback we understand the process of cleanup of the “well-multiple frac-

Email address: [email protected] (A.A. Osiptsov) Preprint submitted to Journal of Petroleum Science and Engineering

tures” system to set up the well for steady-state production after cleanup is complete. 1.1. A problem in field In the field, fracture flowback and well startup operation has long been regarded as a relatively easy job of cleaning the system to meet the following criteria: (1) the well should be cleaned up from the proppant that did not enter the fracture; (2) the fracture should be cleaned up from the fracturing fluid (e.g., broken gel); (3) the well can be connected to the petroleum pipeline system only when the fracturing fluid is cleaned up from the system, which requires additional expenses. Hence, the flowback is simply managed by optimizing the time taken to cleanup the multistage-fractured well system from the fracturing fluid. This means the faster the better. Consequently, flowing back at excessively high rates (large choke sizes or even no choke at all) may lead to excessive drawdown in the December 19, 2019

Figure 1: Scheme of multistage fracturing completion of a near-horizontal well, with placement and flowback phenomena shown on two separate transverse fractures.

near-wellbore and, as a result, to undesired geomechanics phenomena such as: (1) proppant embedment with plastic deformation of the fracture face and creep; (2) proppant pack compaction with crushing of proppant grains and their rearrangement; (3) tensile rock failure due to excessive drawdown in the near-wellbore reservoir area (near-fracture zone); (4) proppant flowback from the fracture into the wellbore with formation of empty cavities, which then close up, pinching out and disconnecting fractures from the wellbore. These geomechanics phenomena are usually neglected, as short-termed effects are favored. However, in long-term perspective this preference might cause the loss in the cumulative volume of oil produced and the overall income deficiency. For these reasons, the full spectrum of the phenomena should be considered during the planning of the well startup. A summary of the geomechanics, fluid mechanics and geochemics phenomena critical for fracture conductivity is given in Tab. 1 [2, 3, 4, 5, 6]. In Fig. 1, we show a detailed scheme of the multi-fractured near-horizontal well. We specifically demonstrate various fluid mechanics and geomechanics phenomena peculiar separately to proppant placement/fracture growth and flowback in the scheme shown in Fig. 1.

1.2. Geomechanics complications of flowback Some of these risks have been recognized earlier. In particular, even in late 1980s [7] it was proposed that severe proppant crushing and proppant flowback damage fracture conductivity. Rapid increase in the choke size may expose the proppant pack to excessively increasing closure stress, which results in crushing. Thus, dynamics of the flowback in time is important. Slow controlled flowback operations with gradual choke size increase (spread over days or weeks with small increment) were recommended for routine use in the field, yet apparently it has not become a widely accepted industry standard. Theoretical studies of well startup are scattering, though some of the work is narrowly focused on the problem of flowback specifically. The work [8, 9] based on a number of field case studies for early time production management for shale gas wells proposed that excessive flow rates during well start up have adverse impact on production. In the last decade, operators and service companies begin to realize the importance of this short early production stage for the cumulative production from the well. In the open scientific literature related to hydraulic fracturing, the vast amount of papers are devoted to the mechanics of fracture creation and propagation (the recent review [10] presents a detailed summary), with much less attention paid to the fluid mechanics part (see the review of the fluid mechanics aspects in [11]), and the fracture flowback and well startup 2

1 2 3 4 5 6 1 2 3 4 5 6 1 2

Geomechanics phenomena Proppant embedment with rock plasticity and creep Elastic compaction of proppant pack Proppant crushing and frictional rearrangement Tensile rock failure at fracture face Proppant flowback from fracture into well Fracture closure on a proppant-free cavity Fluid mechanics phenomena Suspension filtration and colmatation in proppant pack Two-phase frac fluid-oil displacement Fingering in two-dimensional displacement Yield stress of cross-linked frac fluid Transport of fines and development of filter cake in formation and proppant pack Transient pore pressure diffusion in surrounding rock Geochemical phenomena Proppant dissolution Scale formation at proppant grains

Table 1: Key geomechanics, fluid mechanics and geochemical phenomena which have impact on fracture conductivity.

being least studied (see the summary in [11]). The research into a family of geomechanics phenomena critically affecting the fracture conductivity during flowback and production is presented in [12, 13]. In particular, these studies demonstrated an elastic model for proppant embedment based on Hertzian contact problem, a model for elastic compaction of a proppant pack, and a model for a conductivity of a fracture partially filled with proppant (heterogeneous proppant placement resulting from alternate-slug pumping during proppant placement into a fracture). These papers are probably the most relevant to our present work. The deformation of proppant pillars (resulting from alternate-slug fracturing) was also considered in a non-linear formulation in [14], complemented by lab testing of proppant pillars compaction under the closure stress. Conductivity of the fracture propped with rod-shaped particles with effects of embedment is studied in [15]. A lab study of proppant embedment into shale rock is reported in [16].

which reduces the proppant permeability [18]. The effect of fines migration on the reduction of proppant pack permeability was studied experimentally in laboratory conditions. In particular, it was shown that even small volume fraction of crushed proppant grains can reduce the the proppant pack conductivity by above 50% [19]. Nevertheless, there is no general mathematical model capable of predicting the effect of different types of fines mentioned above on the proppant permeability dynamics during fracture cleanup and consequent long-term production. Hydraulic fracturing fluid typically contains polymer additives. Due to fluid leak-off into the rock formation, polymer molecules accumulate at the fracture face and form a polymer cake [22]. After the end of the fracturing treatment, some crosslinked gel may not be fully broken back into the linear gel for efficient cleanup (e.g., some encapsulated breakers did not release, or concentration was incorrect, etc.). Both unbroken gel and polymer cake show yield-stress behavior so that they are displaced from the proppant pack only when a local pressure gradient exceeds a certain threshold value. Development of unyielded zones due to insufficient pressure gradient in a hydraulic fracture reduces significantly the effectiveness of cleanup and the overall productive area of a hydraulic fracture [23, 24, 25]. A technical discussion of the flowback from a multistage fractured well in shales was given in [26]. There are also papers that looked specifically at the transient phenomena during cleanup of a near-horizontal well, without expanding the study into the near-wellbore area of the fracture and focusing on the multiphase flow transients in the pipe [27].

1.3. Fluid mechanics factors of flowback As it is shown in Table 1, there are several fluid mechanics phenomena accompanying fracture flowback, which can potentially reduce the effectiveness of fracture cleanup and long-term fracture conductivity. Different sources of fines are identified during fracture flowback, namely, (i) fines produced from fracture face due to proppant embedment (spalling) or rock failure [17, 16, 6]; (ii) fines mobilized in pores of the formation and transported into the hydraulic fracture (e.g., minerals and swelled or deflocculated clay particles) [18, 16]; (iii) crushed proppant grains due to applied rock confinement stress [19]; (iv) formation of precipitates inside the proppant pack due to incompatible physical and chemical conditions at high pressure and temperature (deposition of insoluble salt or formation of scale due to remineralization of dissolved proppant grains) [5, 20, 21]. Transport of fines through the proppant pack typically leads to colmatation of pores and formation of internal filter cake,

1.4. Our approach In our methodology, we employ a coupled modeling approach, where geomechanics phenomena are considered together in the framework of a global model for fluids displacement in a fracture. Among the solid mechanics effects considered are: proppant embedment with plastic deformations of the rock near the fracture surface, proppant crushing, proppant 3

flowback, proppant pack compaction under the action of the closure stress, tensile rock failure at high filtration rates from the reservoir to the fracture (when drawdown pressure drop is excessively high). These geomechanics effects are implemented into the global fluid flow model in a fracture, which takes into account: displacement of fracturing fluid by oil, two-phase effects, suspension filtration with particle trapping and mobilization, resulting in permeability damage, colmatation and skin buildup inside the fracture in the near-wellbore zone (where particles are mainly formation fines and crushed proppant). Compared to earlier work (e.g., [12, 13]) we take into account a number of significant complexities peculiar to the geomechanics of flowback, namely, (i) in the problem of proppant embedment, we include the plasticity and creep of the rock (where plasticity adds irreversible deformations and creep adds time-dependent deformations at a given loading, both not accounted for earlier in purely elastic models), (ii) in the proppant pack compaction problem, we take into account non-elastic deformation of a packed proppant with account of rearrangement of the grains and their crushing. The ultimate goal is to solve an optimization problem to obtain the most effective regime of the well flowback and startup operation. Aggressive regime is now perceived in the field (e.g., in some oilfields of Western Siberia in Russia) as a better one from the fluid mechanics point of view (more efficient cleanup from the fracturing fluid which was not fully degraded with a breaker), whereas it is dangerous for geomechanics (may cause irreversible geomechanics deformations in and around the fracture). The coupled model presented in the paper is capable of balancing the fluid mechanics and geomechanics phenomena and solving the direct problem of well startup after fracturing. Some preliminary aspects of this research have been presented at a conference and published in a brief communication in proceedings [28]. The applied aspects of flowback technology in the field and the plans for the field testing campaign based on the present modeling workflow are discussed in a conference presentation [29].

2.1. Proppant embedment Below we present a family of models of proppant embedment with increasing complexity, beginning with illustrative numerical simulations and then presenting analytical models for embedment into a rock with effects of plasticity and creep. The proppant is assumed spherical. 2.1.1. Embedment problem with effects of rock plasticity A fracture closed on proppant is shown in Fig. 2, a. Once the well is opened for flowback and cleanup, the pressure is decreased inside the fracture and fracture closes on proppant so that some proppant grains touching the fracture walls may embed into the fracture faces (Fig. 2, b,c), which will result in the reduction of the effective magnitude of the propped fracture width, and, hence, in the decrease of the fracture conductivity.

Figure 2: Proppant placement in a hydraulic fracture (a); fracture closure (b) and proppant embedment into fracture faces (c) (illustrations taken from [30]).

During embedment of a spherical grain into a rock there forms a zone of limiting equilibrium state at some depth into the rock from the fracture face (Fig. 3, a). Then the limiting equilibrium zone expands and covers up almost entire volume underneath the sphere Fig. 3, b).

2. Geomechanics models for fracture flowback

In this section we present a family of geomechanics models for the phenomena that eventually lead to degradation of hydraulic fracture conductivity: (i) proppant embedment into fracture walls with effects of plasticity of the rock; (ii) tensile rock failure at the fracture face, when the drawdown is over-critical; (iii) proppant pack compaction; and (iv) proppant flowback from the fracture into the well. These geomechanics submodels will be then implemented into the global modeling framework for multiphase flow in a closed propped fracture for simulation of the process of fracture flowback and well startup. Note that in what follows we use SI units for all physical quantities throughout the manuscript (unless it is stated otherwise) with the exception of permeability, which is measured in Darcy units (1 D = 0.986 · 10−12 m) convenient to filtration theory.

Figure 3: Formation of the fractured zone at the initial (a) and final (b) stages.

Assume that the proppant grain is an elastic sphere, which enters into the elastoplastic half-space (Fig. 4, a). An appropriate model for this process is a non-associated plastic flow law with the Mohr-Coulomb yield criterion. In order to make sure that the model accurately captures the behavior of the rock it is required to conduct a series of three-axial tests to determine the dependence of the model parameters (Young’s modulus, Poisson’s ratio, cohesion, friction and dilatancy) on the parameters of state (pressure and cumulative plastic deformation). As a result of such pre-processing it is possible to take 4

Friction angle φ 0 4.587 7.584 8.582 11.356 13.297 14.655 15.556 16.125 16.457 16.519 16.445 16.278 16.019

Plastic deformation εipl 0 0.000028 0.000061 0.000082 0.000342 0.000764 0.001298 0.001931 0.002669 0.003695 0.004886 0.005858 0.006876 0.007942

non-associated plastic medium with the Mohr-Coulomb yield criterion as shown in Fig. 4, b. Input parameters of the simulations are as follows: Young modulus E = 1.1 · 1010 Pa, bulk modulus K = 3 · 1010 Pa and cohesion c = 1.5 · 107 Pa. The grain is assumed elastic with the shear modulus E = 8 · 1010 Pa and the bulk modulus K = 1.3 · 1011 Pa.

Figure 5: Numerical scheme for the embedment of an elastic sphere into an elasto-plastic half-space (a) using the strain-softening model.

In Fig. 5 we show the distribution of accumulated plastic deformations at different stages of grain embedment. The most intensive plastic flow is shown by red color. Compared to Figs. 5, a and b, one may conclude that the volume of medium being in plastic flow increases with embedment. These calculations confirm the formation of areas of fractured medium at the initial and final stages of embedment, shown in Fig. 3. However, the use of numerical simulation to determine the depth of grain embedment and the resulting decrease in the fracture width is not convenient as it requires heavy computations. Hence, in what follows we develop an analytical model for grain embedment into the fracture face.

Table 2: The dependence of the friction angle φ on the cumulative plastic strain pl intensity εi .

into account the properties of a specific oil field [31] and use the Mohr-Coulomb model with strengthening, the so-called strainsoftening model [32].

2.1.3. Analytical modeling of proppant embedment into a fracture wall With increase in the embedment depth of a proppant particle, the local deformation zone at the contact zone goes through three stages: elastic, elastoplastic, and plastic [34, 16]. Experimental studies lead to a conclusion that proppant embedment is accompanies with large plastic deformations and this process cannot be described within the framework of the elasticity model [35]. In order to overcome this difficulty one could use the empirical law [30], which related the force and the embedment depth of an indentor, proposed earlier in [36]. From here it follows that theoretical and experimental investigation of the process of proppant embedment in the case of significant decrease in the fracture width can be reduced to consideration of the Brinell problem. Thus embedment of a spherical grain is similar to the problem of the Brinell probe in material characterization. There is vast amount of literature on this characterization technique, e.g. see [37], where the author manages to integrate numerically the equations of axisymmetric problem in the conditions of full plasticity. The paper [38] presents the solution based on plane theory of slip lines in the form of an algebraic relation between the limiting force and the radius of contact for embedment of a hard sphere into the ideally plastic half space. The authors used the distribution of slip lines shown in Fig. 6 taken from [39] for the plane problem on embedment of a cylinder. The slip lines,

Figure 4: Numerical scheme of an embedment of an elastic sphere into an elasto-plastic half-space (a) with the use of the strain-softening model (b).

A typical value of the internal friction angle used in simulations may result in significant errors. We will demonstrate that on an example of processing the results of a lab experiment on three-axial compression of a rock core sample. In Fig. 4, b, we show the curve of core deformation in a real experiment in the press (grey line) and the results of numerical modeling of three-axial core testing. Dashed line corresponds to the numerical result for the core deformation with a constant internal friction angle of φ = 27◦ (according to the lab testing standards). The blue curve corresponds to the numerical result for the deformation curve with account for the dependence of the internal friction angle on the accumulated plastic deformation (Table 2). 2.1.2. Numerical modeling of proppant embedment into a fracture wall In this section we provide illustrative results of direct numerical simulations of the proppant embedment using the programming code FLAC3D, developed for solving 3D geomechanics problems [33]. For calculations, we use the model of 5

on which the shear stresses reach their maximum value of τ = k according to the Mises criterion, are governed by the following equations: σ + 2k∆φ slip = C1 along the α − line

(1)

σ − 2k∆φ slip = C2 along the β − line

(2)

Finally, the limiting traction can be represented as ! " 4π 1 π F = σ s √ r2 − + 3 6 3    a 2 !3/2 1 −1 a − 1 + cos 1− − 3 r r !#  a 3 a a 1 − + 4 −3 4r 36 r r

(6)

In Fig. 7 we show the dependence of the non-dimensional traction on the radius of the contact surface a, created as a result of solving the problem with finite elements for σ s = 108 Pa and E = 3 · 1011 Pa, according to formula (6) and take from the studies [40], [41]. It is seen from the plots that the numerical result is in a very good agreement with the analytical solution (6) and is close to the results by other authors.

Figure 6: Embedment of a sphere into an ideally plastic substrate.

√ Here σ is the stress, k = σ s / 3, σ s – yield strength, C1 C2 constants, ∆φ slip – angle between the x-axis and the tangential curve to the slip line. The angle ∆φ slip is equal to θ shown in Fig. 6. Figure 7: Comparison of the dependence of the non-dimensional force F/πa2 σ s on the non-dimensional contact surface, obtained as a result of the numerical solution (shown with crosses), according to the formula (6) and following [40], [41].

According to the widely accepted assumption that the plane model is applicable to the axisymmetric stress-strain state, the pressure exerted on the contact spherical surface is given by [39]: p = 2k + 2kθ (3)

Solution (6) cannot be used directly for the problem of proppant embedment as the Mises criterion τ = k does not take into account the dependence of hardness on the internal friction, which is peculiar to reservoir rocks. Besides, expression (6) does not include the traction on a free surface which needs to be taken into account when considering the pressurized hydraulic fracture. Consider an embedment of a hard cylinder into an ideally plastic saturated half-space with the Mohr-Coulomb limiting condition (Fig. 8):

Relation (3) allows one to calculate the vertical limiting force: F = 2π

Za xp sin(θ)dx

(4)

0

Taking into account the relation (3) and the expression x = r cos(θ) the integral (4) can be rewritten as:

F = 2πr

Zθ1

σ + σ 2  σ − σ 2 x y x y + τ2xy = + ce f cot φ sin2 φ 2 2

(2k + 2kθ) sin θ cos θdθ = 2

2 π/2

" 4kπr

2

sin θ 1 + (12θ sin3 θ + 9 cos θ − cos 3θ 3 36 3

where ce f = c− pr tan φ – effective cohesion, c – cohesion in the conditions of drainage at p = 0, φ – angle of internal friction, pr – pore pressure. As it is shown in Fig. 9, the cylinder with a radius R is introduced into the half-space with a certain penetration depth and

#θ1 ,

(7)

(5)

π 2

where θ1 is the angle at x = a. 6

expression (10) can be rewritten in the form Pu = −2πr

Zθ0 " 2 π/2

Eventually, after integrating we obtain the expression for the limiting force shown in Eq. (14). The plots of the limiting force as a function of the depth of embedment of the sphere are shown in Fig. 10. According to Fig. 10, the depth of embedment depends significantly on the internal friction angle. The analysis carried out so far allows us to determine the depth of proppant embedment into fracture faces and the associated fracture closure (the decrease in the fracture width). We note that there is no analytical solution available to this problem in the literature. Typically one would use empirical relations based on the Hertz solution for an elastic sphere and a half space, which is then tuned to match lab data. As an example, we will consider the recent study on fracture closure [43] with a monolayer of proppant (Fig. 11). The Hertz solution for a half-space was used first !2/3  !2/3 !2/3   1  1 − ν2 1 − ν2s 1 − ν2 3  (15) + − eh = 2 pD1  8 D1 E Es E

Figure 8: Mohr circle for the fracture face.

a uniform pressure distribution over the surface. The slip lines (where the shear stress reaches its maximum value) are governed by the equations [42]: σ0 e−2ω tan φ = C1 along the α − line

(8)

σ0 e−2ω tan φ = C2 along the β − line

(9)

where σ0 = ce f cot φ + σ, c and φ – the cohesion coefficient and the internal friction coefficient, respectively, and σ = (σ x + σy )/2 – mean stress, C1 and C2 are constants, ω is the angle calculated clockwise from the x-axis to the main principal stress plane (Fig. 8).

where E and ν are Young modulus and Poisson’s ratio of the proppant particle, E s and ν s are Young modulus and Poisson’s ratio for the fracture wall (formation). Other notation in Fig. 11. Calculations according to formula (15) and lab data on embedment are shown in Fig. 12, a. It can be seen that the numerical results (blue curve) differ significantly from the lab data (blue dots). Hence it was proposed to modify equation (15) and introduce extra terms b0 and b1 ,   !   1 − ν2 1 − ν2s 2/3 + − eh = b0 + b1 1.04D1 (K 2 p)2/3  E Es  (16) !2/3    1 − ν2 p  + D2  − E E

Figure 9: The field of slip lines in the rock, according to the Mohr-Coulomb criterion.

The limiting traction (maximum for this type of rock) can be expressed in the following form: Pu = 2π

which finally help to ensure there is a match with experiment (Fig. 12, b). In conclusion we will note that, based on solving the boundary-value plasticity problem, we obtained the final expression (13), which related the fracture opening with the solid mechanics properties of the rock, the pressure inside the fracture, proppant grain size and proppant concentration. The geomechanics properties of the rock are obtained based on the three-axial core testing. In contrast to the frequently employed semi-empirical correlations based on the modified elastic Hertz solution, we propose the solution (14), which is adequate to the realistic behavior of the embedment process, as it is confirmed by the numerical simulation of the penetration of a grain into the fracture face (Fig. 4). We will now consider the embedment of a layer of proppant

Za xσE,y dx

(10)

0

where σE,y is the vertical component of a contact stress σOE . The vertical component of the stress on the contact surface takes the form [42]: σE,y = τ sin φ + σ cos φ = [σOE (1 + sin φ) − ce f cot φ] cos φ.

(11)

Taking into account that σOE =

p s + ce f cot φ exp(2θ tan φ) 1 − sin φ

p f + ce f cot φ (1 + cos θ) exp(2θ tan φ)− 1 − sin φ (13) i −ce f cot φ sin2 θ cos θdθ

(12) 7

h   pi + ce f cot φ e2 tan φθ0 2 tan φ cos3 θ0 + 3 cos2 θ0 sin θ0 − 2 (1 − sin φ)(4 tan φ + 9) # # 2 4 tan φ + 3 2 tan φθ0 4 tan2 φ + 3 tan φπ 1 3 (2 ) − c cot φ sin φ(1 − sin θ ) e tan φ cos θ + sin θ + e − ef 0 0 0 3 4 tan2 φ + 1 4 tan2 φ + 1 "

Pu = 2πr2

(14)

Eq. (14) to the form ∆p s =

A C

 [ph + (c cot φ − pr )] − B /

L2 A + 2 C 2πr

! (18)

h   A = e2θ0 tan φ 2 tan φ cos3 θ0 + 3 cos2 θ0 sin θ0 + 4 tan2 φ + 3 tan φπ e − 4 tan2 φ + 1 # 4 tan2 φ + 3 2θ0 tan φ (2 tan φ cos θ0 + sin θ0 ) − e 4 tan2 φ + 1 1 B = ce f cot φ sin φ(1 − sin3 θ0 ) 3 C = (1 − sin φ)(4 tan2 φ + 9). +

Figure 10: The dependence of the normalized force Pu /cπa2 on the depth of embedment e0 /r for φ = 0 (curve 1) and φ = 15◦ (curve 2).

The expression (18) provides the relation between the angle of embedment θ0 and the applied pressure drop ∆p s , while the expression for the embedment depth in terms of the angle θ0 is as follows: r e0 θ0 = arcsin 1 − r We will now define typical values of the geomechanics properties of the rock. Assume the values typical of an oilfield in Western Siberia, Russia at the depth of the order of 3 km. For clays: ν = 0.25, c = 14.4 · 106 Pa, φ = 30◦ , ρ = 2000 kg/m3 , pr = 3.3 · 107 Pa. For sandstones: ν = 0.25, c = 36.4 · 106 Pa, φ = 50◦ , ρ = 2500 kg/m3 , pr = 3.3 · 107 Pa. The initial stress-strain state at the depth of about 3 km is evaluated using the Dinnik formula [44]:

Figure 11: Model of a fracture with a monolayer of proppant.

σ033 = σ022 = into the fracture wall (see Fig.13) with account for Eq. (14). When the pressure inside the fracture is decreased, the proppant particles are exposed to the load equal to the pressure drop ∆p s . We will assume that at the initial pressure ph the proppant fills the entire open fracture but experiences no loading from the fracture faces (the pressure inside the fracture is in equilibrium with the closure stresses). When the pressure inside the fracture is decreased by the magnitude of ∆p s , the fluid pressure inside the fracture is equal to p f = ph − ∆p s (Fig.13, b). In this case the loading from the fracture faces is transmitted to proppant grains and each particle is under the action of the embedding force equal to Pu = ∆p s L2 (17)

ν σ0 . 1 − ν 11

For clays: σ011 = −6 · 107 Pa, σ022 = σ033 = −2 · 107 Pa. For sandstones: σ011 = −7.5 · 107 Pa, σ022 = σ033 = −2.5 · 107 Pa. For calculating the hydraulic fracturing pressure we will use the formula [45]: ph = −

ν (σ0 + pr ) + pr . 1 − ν 11

(19)

This formula yields ph = 4.2 · 107 Pa for clays and ph = 4.7 · 107 Pa for sandstones. The plots of the proppant embedment depth e0 /r with the decrease in the fluid pressure inside the fracture p f are shown in Fig. 14. It is clear from the plots that the dependent of embedment depth on pressure is nonlinear. Proppant is embedded deeper into a soft clay than into a hard sandstone, at the same pressure drop. In a loosely packed layer

For calculation of the dimensionless embedment depth e0 /r under the pressure drop ∆p s , we will transform the expression 8

Figure 12: Dependence of the depth of proppant embedment on the pressure drop: numerical results and lab data (blue dots). Blue curve according to the classic Hertz model (a) and red curve according to the Hertz model with empirical extra terms (b).

Figure 13: Regular proppant distribution in a fracture (a) and the numerical scheme (b). Figure 14: Depth of proppant embedment into the fracture face e0 /r at the total vertical depth of 3 km with the decrease in the fluid pressure p f for clays (orange lines) and sandstones (gray lines) for the close packing (L = 2r, solid lines) and for the loose packing (L = 4r, dashed lines).

of proppant (L > 2r), the depth of embedment is larger than in the close packing. We note that with the decrease in cohesion, the threshold pressure, at which proppant noticeably penetrates into the rock, decreases. Thus, the analytical solution to the boundary-value problem of plasticity theory obtained above allows us to get the algebraic relation for the fracture opening (which is the initial opening minus the twice the depth of proppant embedment) with the geomechanics properties of the rock, pore pressure, fluid pressure in the fracture, proppant size and its concentration.

elastic properties of the rock [46]. An example of the calculated creeping curve obtained using the FLAC3D programming code [33, 47] is shown in Fig. 15 (after [48]). However, the process of particle embedment is not elastic. It is particularly clear from the microphotograph of the fracture face near the proppant grains penetrating into the rock surface (Fig. 16, a, after [17]). An explicit evidence for plastic deformations is irreversible residual deformations and rock spalling (Fig. 16, b, [49]), the latter being the source of formation fines (small particles), which are then transferred during flowback towards the well and may result in proppant pack colmatation, fracture conductivity damage and skin. Based on this evidence it is not unreasonable to assume that the transient evolution of proppant embedment in time will occur due to the gradual microcracking of the rock around proppant grains and, hence, the decrease in the rock cohesion. Experiments with rock samples demonstrate that rock failure under the conditions of controlled deformation occurs gradually and is related to the accumulation of microcracks [50]. In order to take into account this phenomenon in constitutive relationships, the material damage function ωdmg ≥ 0 was intro-

2.1.4. Proppant embedment with effects of rock creeping Embedment of proppant into the fracture face may be qualitatively split into two stages. First, particles penetrate into the rock at a certain depth under the action of the loading induced by the pressure drop. This embedment occurs instantaneously with the decrease in pressure. Then, at the second stage during the long-term operation (production through the fracture) the geomechanics properties of the rock change with time, even under the action of the same loading (pressure drop). This change in the mechanical rock properties induces gradual penetration of the particles into the rock as a result of rock creeping. Usually one assumes that the initial embedment is due to elastic and plastic deformations of the rock, while the further embedment due to creep is conditioned by the visco9

Figure 17: Sketch of the lab setup for creep tests.

∆ε0 = ε∗0 − εini 0 0.2 0.3 0.4 0.5

Figure 15: Numerical results of finite-element modeling of the fracture closure on proppant with effects of rock creeping in the Haynseville shale, after [48].

ccr , Pa 1.063 · 107 9.902 · 106 9.242 · 106 8.594 · 106

p. Pa 2.702 · 107 3.365 · 107 3.775 · 107 3.98 · 107

Table 3: Embedment increments as function of the loading and cohesion.

and the rate of embedment is proportional to the initial embedment eini 0 (Fig. 17), a, eini k= 0 (22) Tr

Figure 16: Microphotograph of the fracture wall after proppant embedment with effects of spalling [17] (a); the sketch of rock spalling [49] (b).

where T r = 0.25 years is a characteristic time scale.

duced [51, 52] (this concept was later summarized in the publication in English [53], [54]). This material damage function ωdmg is equal to zero at the initial time instant and is equal to unity at the time instant of failure. If ωdmg does not depend on time, then we have a plasticity problem. If the dependence of failure on time is specified, then the relations obtained describe the material creep problem. The kinetic equation of failure accumulation with time can be represented in the form [53]: Qωdmg = F(ei ),

Figure 18: Creep curves (a) and plots of the variation of cohesion with time (b). Curves 1-4 correspond to initial embedment of the grain εini 0 = 0.1, 0.2, 0.3, 0.4.

(20)

Variation of cohesion with time may be calculated according to the formula, following from (14):

where Q is a certain differential operator in time, and ei is the strain intensity. According to Eq. (20), we will assume that creep is due to the decrease in cohesion between the material particles in time c(t) with gradual accumulation of microcracks under the load and the increase in the strain intensity ei (t). As a characteristic measure of the accumulated strain we may adopt the nondimensional embedment ε0 (t) = e0 (t)/r and assume that c(t) = c(ε).

c=p

C , 2 cot φ(A − BC)

p=

Pu 2πr2

(23)

Corresponding curves for the material with parameters c = 1.44 · 107 Pa and φ = 30◦ are shown in Fig. 18, b. The time of embedment of proppant in case of constant initial loading Pu is fully defined by experimental creep curves, shown in Fig. 18, a. However, if the loading is changing at some time instant, then the embedment should be recalculated with the new revised value of ccr . Consider a synthetic case. Assume that embedment reaches the maximum value of ε∗0 = 0.6 at the time instant t = 0.225. Table 3 contains the list of embedment increments ∆ε0 = ε∗0 − εini 0 with the increase in the load p and cohesion ccr , which are in agreement with the magnitude of maximum embedment of ε∗0 = 0.6

(21)

To evaluate the time tmax required for the proppant grain to penetrate into the rock at the distance equal to its radius (so that ε0 = 1), one needs to plot a series of creep curves ε0 (t) for various initial instantaneous values of εini and respective magnitudes of the pressure drop ∆pini s . Consider the embedment of a rigid sphere into the rock sample under the constant loading, starting from the initial embedment at a specified depth (Fig. 17). Assume that the creep curves are linear (as shown in Fig. 15) 10

2.2. Tensile rock failure ν21 ν31 1 (σ11 + pr ) − (σ22 + pr ) − (σ33 + pr ) , E1 E2 E3 1 ε23 = σ23 , 2G23 ν12 1 ν32 (σ22 + pr ) + (σ22 + pr ) − (σ33 + pr ) , ε22 = − E1 E2 E3 (25) 1 ε23 = σ13 , 2G13 ν23 1 ν13 (σ11 + pr ) − (σ22 + pr ) + (σ33 + pr ) , ε33 = − E1 E2 E3 1 ε12 = σ12 , 2G12 ε11 =

For solving various applied problems of oilfield development, there is a need to determine the initial stress-strain state of the reservoir, and also its change in the course of pressure increase (e.g., during fracturing) and decrease (e.g., depletion during production). Estimates for initial state of the ”dry” reservoir rock are known for decades, beginning with the Dinnik formula [44] for an isotropic medium, which was later generalized for anisotropic distribution of the geomechanics properties [55], [56]. In the reservoir saturated with the fluid, the full stress depends on the variation of pore pressure. There are direct methods to determine the stress at various stages of field development, based on fracturing technology. For example, for an unconsolidated sand in the Gulf of Mexico located at a depth of 3700 m an empirical correlation was obtained, which relates the current value of the horizontal stress σ11 with its initial value σ011 and the pore pressure decrease in the reservoir ∆pr in the process of development [57]: σ11 = σ011 + ∆pr /2

E1 ν21 = E2 ν12 , E2 ν32 = E3 ν23 , E3 ν13 = E1 ν31 . Since ε11 = ε22 = 0, ε12 = ε13 = ε23 = 0 (Fig. 19, b), from (25) the following equations can be obtained: σ11 = λ1 (σ33 + pr ) − pr , σ22 = λ2 (σ33 + pr ) − pr

(26)

E2 ν32 + ν12 ν31 E1 ν31 + ν21 ν32 , λ2 = . E3 1 − ν12 ν21 E3 1 − ν12 ν21

(27)

λ1 =

(24)

For a transversely-isotopic medium E1 = E2 = E, E3 = E 0 , ν21 = ν12 = ν0 , ν31 = ν32 = ν0 and, hence, horizontal stresses are equal σ11 = σ22 . In this case

We consider a hydraulic fracture in a productive later with initial stress state σ011 , σ022 , σ033 (σ033 < σ011 , σ022 ) and pore pressure pr . During the pumping process and subsequent field development, the stress state and pore pressure vary up to the current values σ11 , σ22 , σ33 , pr , respectively. Pumping of suspension into a fracture is carried out at the pressure ph > |σ033 | (Fig. 20). After the pumping process is complete, the pressure is decreased by the magnitude ∆pr – this load is exerted on proppant which prevents the fracture from complete closure. We distinguish two different pressure drops: ∆p s – at the fracture wall, while ∆pr is the variation of the pore pressure.

E ν0 . (28) E0 1 − ν In the case of an isotropic medium E 0 = E, ν0 = ν and from (28) one could find a widely known lateral pressure coefficient [44]: λ1 = λ2 = λ =

ν . (29) 1−ν Note, that the decrease in the pore pressure leads to the increase of the principal shear stresses and, consequently, to the decrease in the rock strength. λ=

We will assume that the stress state at the fracture walls at a certain distance from the fracture tip depends only on the coordinate x3 and is statically determinate.

σ11 − σ33 1 = (λ1 − 1) (σ33 + pr ) , 2 2 σ22 − σ33 1 T2 = = (λ2 − 1) (σ33 + pr ) 2 2

T1 =

(30)

The decrease of pressure by the magnitude ∆p leads to the increase in shear stress and the decrease in the rock strength. We will now compare the decrease in the stress σ11 with empirical relation (24). In the case of an isotropic medium, the initial value of the horizontal stress σ011 takes the form according to Eq. (26):

Figure 19: Numerical scheme for the fracture face (a) and an element of the porous medium saturated with fluid (b).

σ011 =

Consider an orthotropic medium with principal axes along the coordinate axes of the half-space. In this case equations in the expanded form can written as [58]:

 ν  0 σ33 + p − pr 1−ν

(31)

Hence, the current value of the stress σ11 for p = p0 − ∆p can be rewritten as 11

1 − 2ν ∆pr (32) 1−ν Comparing the relations (32) and (24), we note that these expressions almost coincide, if the Poisson’s ratio is in the range 0.3 ≤ ν ≤ 0.35, typical of reservoir rocks. σ11 = σ011 +

The hydraulic fracture opening pressure can be determined from the formula (19). Substituting Eq. (35) into the condition (33), we find the equation for ∆pcr s in terms of the initial rock stress components σ011 and σ022 : AP2 + BP + C = 0,

2.2.1. Tensile rock failure in 1D formulation

(38)

where 1 A = [(1 − λ)2 − α(1 + 2λ)2 ], 3 0 0 B = [(1 − λ)(σ11 + σ22 ) − 2cα(1 + 2λ)+ 2 + α(σ011 + σ022 )(1 + 2λ)] 3 1 2 0 0 C = −[3c − 2cα(σ11 + σ22 ) + α(σ011 + σ022 )2 ]+ 3 ((σ011 )2 − σ011 σ022 + (σ022 )2 )

P = −∆p s + (ph − pr ),

Figure 20: A hydraulic fracture.

Proppant distribution in a fracture after closure may be heterogeneous, including proppant pillars and proppant-free cavities (e.g, an alternate-slug fracturing technique [59]). Along the unsupported cavity, one can neglect the variation of the stress state Fig. 20, a. Consider the stress state of the fracture wall after the pressure is dropped down to the value of p f = ph − ∆p s Fig. 20, b. As a result of the tensile load (“stretching” the rock) ∆p s , the shear stress is increased and the rock is close to the failure threshold. As long as the strains are elastic, the fracture face remains undamaged. The tensile failure occurs when the strains become plastic. We will use the Drucker-Prager yield criterion for analysis T + ασe f f = c,

σe f f =

1 ef f σ δi j 3 ij

Solving Eq. (38), we find √ B + B2 − 4AC cr + (ph − pr ) ∆p s = 2A

(39)

We will now make illustrative simulations using representative values of the geomechanics parameters of the rock, typical of an oilfield in Western Siberia, Russia at the depth of 3 km (the same parameters as those used in Sec. 2.1).

(33)

where 1 p (σ11 − σ22 )2 + (σ22 − σ33 )2 + (σ33 − σ11 )2 T= √ 6

(34)

Here T is the shear stress intensity, c – cohesion, α = sin ϕ – internal friction coefficient, ϕ – internal friction angle, σe f f is the mean stress. The stresses at the fracture wall are determined by the formulas:

Figure 21: Dependence of the nondimensional pressure drop pcr s /ph on the 0(e f f ) initial stress |σ22 | for clays (a) and sandstones (b) of the Priobskoe field. Yellow color shows the zone of fracture face failure.

σe33f f = −(ph − ∆p s ) + pr , f f) σe11f f = σ0(e − λ1 (ph − ∆p s − pr ), 11

(35)

f f) σe22f f = σ0(e − λ1 (ph − ∆p s − pr ), 22 f f) f f) where σ0(e , σ0(e are initial stresses due to the rock weight 11 22 and tectonic stresses, and λ = ν/(1 − ν). The initial stress state is determined using the Dinnik formula [44]: f f) σ0(e 22

ν = σ0(e f f ) , 1 − ν 11

Figure 22: Dependence of the nondimensional pressure drop pcr s /ph on the 0(e f f ) initial stress |σ22 | for clays (a) and sandstones (b) of the Priobskoe field (the cohesion is reduced by a factor of two, as compared to the results shown in Fig. 21). Yellow color shows the zone of fracture face failure.

(36)

f f) where the vertical effective stress σ0(e (overburden stress) at 11 the depth hdepth is determined by the weight of overlying rock: f f) σ0(e = −ρghdepth + pr 11

We now calculate the variation of ∆pcr s with the depth h. According to Fig. 21, the relative critical pressure drop decreases with the depth h. Above the threshold red line, the zone of frac-

(37) 12

ture face failure is located. The decrease in cohesion results in fracture face failure at significantly smaller critical pressure drop (Fig. 22). We note that the expansion of the pressure diffusion front from the fracture face into ambient rock, the zone of plastic deformation will expand as well.

medium with effective moduli K and E (Young’s modulus), depending on the porosity φ. Thus the full deformation ε11 can be represented as ε11 =

2.2.2. Tensile rock failure in 2D formulation Hydraulic fractures are subject to failure while flowback and well startup (Fig. 23), where the pressure is decreased from ph to p f = ph − ∆p s . If the choke is opened too quickly and the bottomhole pressure decreases rapidly, that may lead to proppant pack disintegration and proppant flowback from the fracture into the well. The wells, which are brought to production according to an aggressive scenario of choke management, may produce a large amount of solid particles with the size ranging from fines (0.01mm) to proppant grains (0.1mm-1mm) to pieces of rock (1cm) [60]. For analysis of this process we propose to consider an elastic plane problem of a cavity in the layer of proppant (Fig. 23, b), using the relations for a coupled poroelasticity Biot problem shown in [61].

 
p 1 (σ11 + p) − ν (σ22 + p) + (σ33 + p) − E 3K s

(40)

Using the known relation K=

E 3(1 − 2ν)

(41)

the expression (40) can be rewritten as

ε11 =

! p K 1 [σ11 − ν(σ22 + σ33 )] + 1− . E 3K Ks

(42)

Introducing the notation α=1−

K , Ks

(43)

we eventually write ε11 =

 
1 (σ11 + αp) − ν (σ22 + αp) + (σ33 + αp) E

(44)

Similarly, for the other two components of the strain tensor  
1 (σ22 + αp) − ν (σ11 + αp) + (σ33 + αp) (45) E  
1 = (σ33 + αp) − ν (σ11 + αp) + (σ22 + αp) E

ε22 = ε33

Thus, in order to find the Biot coefficient α one needs to know the effective elastic moduli of dry porous medium. We write the mass balance equation [61]

Figure 23: Pressure decline in the hydraulic fracture during flowback (an aggressive scenario) and well startup (a) and numerical scheme for the cavity problem (b).

We will now illustrate the concept of poroelasticity on an example of the medium element with a single pore containing fluid under pressure p [62]. We decompose the full stress state into two as shown in Fig. 24.

M=

1 ∂p k 2 ∂ε = ∇ p−α M ∂t η ∂t

(46)

Kf ϕr + (α − ϕr )(1 − α)K f /K

(47)

where M is Biot modlus, K f is bulk modulus of the fluid, ϕr is porosity. Equation (46) coincides with equation (5.15) in the book [63]. Variation of the fluid volume fraction in the unit volume of porous medium ζ is related to the bulk strain ε = εkk by the following formulae p = M(ζ − αε)

(48)

If there is no drainage (fluid filtration), then ζ = 0 and pore pressure p depends only on the bulk strain.

Figure 24: Stress state of the element with a single pore.

According to the illustration, the deformation is composed of the deformation of the continuum element with the modulus K s of the matrix material and the deformation of a dry porous

p = −Mαε In this case 13

(49)

Parameter b n L Sv Sh pr ∆p s ϕr φ E ν

Magnitude 0.01[m] 5 10 [m] 89.7 [MPa] 75.9 [MPa] 69 [MPa] 34.5 [MPa] 0.05 30 34.5 [GPa] 0.2

Description Fracture half-width Fracture aspect ratio Size of the computational domain Vertical stresses (13 kpsi) Horizontal stresses (11 kpsi) Initial pore pressure (10 kpsi) Amplitude of the pressure drop Rock porosity Friction angle Young Modulus (5 Mpsi) Poisson’s ratio

Table 4: Main initial simulation parameters [60].

(50)

Figure 25: Distribution of the strength function F MC in drainage conditions. Pore pressure is uniform.

The relations shown above needs to be complemented by the formulas, which relate the elastic moduli of the matrix material (E s and ν s ) and effective properties of the rock (E and ν). For this purpose we will use the relations [64]

decreases rapidly in aggressive scenario. In this case the pressure is increased (Fig. 26, a) and, consequently, the effective strength is decreased, which leads to formation of the failure zone near the narrow end of the cavity, where F MC < 0 (Fig. 26, b).

σi j = 2Gεi j + K ∗ εδi j 2Gν K∗ = + α2 M 1 − 2ν

" #−1 E ϕ 3(1 − ν s )(9 + 5ν s ) = 1+ , Es 1−ϕ 2(7 − 5ν s ) " #−1 E ϕ 3(1 − ν s )(1 + 5ν s ) ν= νs + Es 1−ϕ 2(7 − 5ν s )

(51)

The main initial parameters are shown in Table 4, which is taken from [60]. For simplicity, it is assumed that the properties of the layer do not differ from the embedding medium. Using the relations shown above we will calculate the Biot parameters: α = 0.095, M = 5.3 · 1010 Pa. Note that the coupled problem is related not to the initial pore pressure pr , but to the pressure that is decreased ∆p s . We will run the simulation for two limiting states. First we will assume that the pressure in the fracture is decreased slowly so that pore pressure is uniform (Fig. 25). In order to evaluate if the stresses overcome the tensile strength we will use the Mohr-Coulomb and write the expression:

F MC

   σ1 − σ3 σe1 f f + σe3 f f  = c cos ψ −  + sin φ 2 2

Figure 26: Distribution of pore pressure in conditions without drainage (a) and the strength function F MC (b). Failure zone is highlighted (F MC < 0).

Thus, in conditions without drainage (pressure decrease rate above filtration velocity) are favorable for proppant pack disintegration. In order to determine the safe operating envelope in terms of the ratio of velocities one needs to include the mass balance equation (46). 2.3. Proppant pack compaction The porosity of a proppant pack in a fracture under closure stress depends on the load. Porosity decrease results in the fracture conductivity degradation and, hence, limits the productivity of the well [65]. We consider the variation of the proppant pack porosity under the action of the closure stress, which results in the fracture width decrease, and hence, in the fracture conductivity deterioration. In what follows, we aim at the dependence of the proppant pack inelastic deformation and the closure stress imposed on the proppant pack within the framework of the elasto-plastic model of a proppant pack with internal friction and dilatancy. Consider an elementary volume of a face-centered lattice (Fig. 27, a) and calculate the corresponding porosity φ.

(52)

where σ1 , σ3 are principal stresses, c = 9.54 · 107 Pa is cohesion. If F MC < 0, then the strain state is beyond the strength limit (threshold). According to Fig. 25 the function of strength F MC > 0 and, consequently, with slow pressure decrease the proppant crushing does not occur. The case is very different when the bottom hole pressure is 14

According to Eq. (53) we find the variation of porosity of the proppant pack !# " 1 π (56) ∆ϕ p = ε 1 − ε 1 − ε , 2 6 where ∆ϕ p = ϕ0p − ϕ p . According to the Hertz solution 3(1 − ν2 ) ε= P 4Er2 "

Figure 27: Elementary cell of a face-centered lattice under the closure stress.

(53)

#2/3 " ∆p s . ε = 3(1 − ν2 ) E

we find the expression for the initial proppant pack porosity ϕ0p in the form: ϕ0p = 1 −

π = 0.476 6

,

(57)

where E and ν are the Young’s modulus and Poisson’s ratio of the proppant grain material. Since the force P is related to the variation of pressure in the fracture ∆p s by the formula P = 4r2 ∆p s , we will finally write the expression (56) in the form

Using the expression for the pore volume π 4 V p0 = (2r)3 − πr3 = 8r3 (1 − ) = 3.811r2 3 6

#2/3

(58)

Thus, once the value of ε is determined, we can calculate the new value of porosity as a result of elastic compaction, and, hence, the modified magnitude of the fracture width

(54)

The value of porosity (54) corresponds to the random close packing of non-deformable spheres. With compaction of the pack, the grains come closer at the displacement u and a contact area of the radius a is formed (see Fig. 28). The displacement ε and the radius of the contact area a are related as √ p a = 2r ε(1 − ε), (55)

w = w0 (1 − ε).

(59)

Consider the proppant grains with Young’s modulus E = 1 × 1010 ÷ 3 × 1010 Pa and Poisson’s ratio ν = 0.25. In Fig. 29, we plot the curves of porosity variation ∆ϕ p as a function of the pressure drop ∆p s .

where ε = u/2r is the deformation of the pack under compaction.

Figure 29: Dependence of the variation of porosity ∆ϕ p on the pressure drop ∆p s for various values of Young’s modulus E.

In the first-order approximation, Eqs. (57) and (59) can be applied for elastic compaction of all possible proppant grain packings. The relations obtained allow one to determine the effective Young’s modulus of the proppant pack: p3 E 2 ∆p s ∆p s Ee f f = = p3 (60) ε 9(1 − ν2 )2 Figure 28: Dependence of porosity on strain.

According to Eqs. (58) and (60), the effective modulus Ee f f 15

depends nonlinearly on strain ε (Fig. 30).

Figure 31: Proppant pack in the initial undeformed closely packed state (a) and the final compacted state (b).

We will assume that Poisson’s ratio is zero for elastic strains ν = 0, since the proppant pack is in the conditions of monoaxial stress state σ33 = −∆p s and there is compaction. In this case from the constitutive relations (64) it follows that

Figure 30: The curves Ee f f (ε) for various values of the elasticity modulus E.

The face-centered lattice considered above has the initial porosity ϕ0p = 0.476. Under the loading, proppant particles may be displaced and rearranged with an inelastic shear deformation (see Fig. 31,a) γ shear = tan(ψ) (61)

dσ33 = E3333 dε33 , where

Then the configuration shown in Fig. 31, a) will evolve into the packing with a minimum porosity ϕ p = 0.258 (Fig. 31, b). In this case, the bulk strain occurs θ = ϕ0p − ϕ p .

(66)

" A= 1−

(62)

E3333 = 2GA ! !# 1 1 2 1 − √ Λ 1 − √ αfr . 3 + 2α f r Λ 3 3

(67)

Assuming G = Ee f f /2 and integrating Eq. (66), we obtain finally the relation between the closure stress on the proppant pack and its strain: r ∆p2 3 3 (68) 9(1 − ν2 )2 2s . ε33 = 2A E

The coefficient of dilatancy (compaction) can now be found in the form: θ Λ= = −0.378 (63) γ shear Consider the deformation of a proppant pack within the framework of an effective medium model. For analysis of inelastic behavior, we will use the Prandtl-Reuss equations for material with dilatancy and internal friction [66], which relate the components of the stress tensor differential and the strain tensor differential: dσi j = Ei jkl dεkl . (64)

We will represent the new magnitude of porosity in the form ϕ p = ϕ0p − ε33

(69)

The reduced fracture width in this case is re-calculated by the formula: w = w0 (1 − ε33 ) (70)

The rigidity coefficients of the Prandtl-Reuss constitutive relation (64) have the form (similar to Young’s modulus in a uniaxial tension/compression in a purely elastic case): ! " # K 2 Ei jkl = G (δik δ jl + δil δ jk ) + − δkl δi j − G 3 (65) s s ! G K K ij kl − + Λδi j + α f r δkl , G + α f r ΛK T G T G

We will select the values of Young’s modulus of proppant, following [67]. Consider the interval 3·1010 Pa < E < 7·1010 Pa and plot the curves of deformation of the proppant pack for the friction coefficient α f r = 0.34 and the Poisson’s ratio ν = 0.17 (Fig. 32). In [68], the proppant deformation curves have been obtained using the PFC3D programming code [69] (red curve in Fig. 33). It is clear that simulation and lab test agree in the linear compaction phase, until crushing of proppant begins. In Fig. 33, b we show the comparison of the lab test and the analytical solution (68), for ν = 0.17, E = 5 · 1010 Pa, α f r = 0.34. The plot shows that the theory and experiment agree well until crushing of proppant begins at ε33 > 0.06. The picture of a proppant grain crushed under the action of the closure stress of 50 MPa is shown in Fig. 34 [70]. The con-

where G is the shear modulus, K = 2G(1 + ν)/3(1 − 2ν) is the bulk compression modulus, ν – Poisson’s ratio, si j = σi j − δi j σ – components of the deviatoric stress tensor, σ = σii /3 – mean stresses, T = (si j si j /2)2 – shear stress intensity, and α f r is the internal friction coefficient. Here the parameter Λ is the dilatancy coefficient determined for the proppant pack by relation (63). 16

Figure 34: The photo of a crushed proppant grain, after [70].

Figure 32: Deformation curves of a proppant pack for E = 3 · 107 Pa (curve 1), E = 5 · 107 Pa (curve 2) and E = 7 · 107 Pa (curve 3).

Figure 35: Lab curve of proppant compaction [68] and the curve by analytical formula (74). Figure 33: Numerical (a) and experimental (b) study of proppant pack compaction, after [68].

Substituting (73) into (72), we finally write: tact surface is highly fragmented and the grain is decomposed into several parts. In order to take this phenomenon into account in the present model, we introduce the parameter of damage ωdmg , which takes into account that the surface of the effective cross-section is increased and the effective stress grows as well [53]: σ σe f f = (71) 1 − ωdmg

∆p s =

 E 2A B n  ε 1 − ε 33 n + 1 33 3(1 − ν2 ) 3

!3/2 .

(74)

Assuming n = 2 and B = 56, we will plot the compaction curve with account for proppant crushing, which is shown in Fig. 34 by the grey line. It is clear from the Figure that the analytical formula with account for proppant crushing captured by the damage theory is now in a good agreement with experiment.

Taking into account formula (71), we will transform (68) into ! Z E 2A ε33 ∆p s = (1 − ω (ξ))dξ . (72) dmg 3(1 − ν2 ) 3 0

We have presented an integrated modeling workflow for the coupled fluid mechanics/geomechanics problem of fracture flowback. A family of geomechanics models is proposed for various solid mechanics phenomena that have profound impact on fracture conductivity, namely proppant embedment into the fracture faces with account for rock plasticity, proppant pack compaction (inlcuding elastic compaction, inelastic frictional rearrangement of particles and crushing), tensile rock failure and proppant flowback.

Assume that the parameter of damage is the function of strain [71]: ωdmg = Bεn33 , (73) where B and n are tuning parameters. 17

3. Multiphase flow in a fracture with geomechanics effects   k p = k0p 1 −

In the current section we formulate the mathematical model of a filtration in a propped fracture during flowback with account of geomechanics (non-elastic proppant pack compaction and proppant embedment) and fluid mechanics effects described above (multiphase filtration, fines transport and proppant pack colmatation). The aim of this modelling is to study the fracture conductivity dynamics during flowback, and in particular, to investigate the relative importance of various factors contributing to fracture conductivity degradation. We assume that during flowback the hydraulic fracture is closed on proppant pack, which is filled initially with the fracturing fluid. Before the start of flowback, the pore pressure inside the proppant pack is equal to the pore pressure inside the surrounding rock matrix, so that the “fracture-rock” system is in equilibrium. Then the bottom-hole pressure in the well decreases according to a certain dynamics of the choke opening at the surface, which triggers the influx of reservoir fluid through the fracture walls. For simplicity, we consider the flowback in a single wing of a hydraulic fracture approximated by a 1D model representation, which assumes that the fracture is elongated in the horizontal direction.

=

σtr , − Cmax

ϕ pf l

=

ϕ0p

− σtr .

(76)

∂ f f ∂ f [ρ si (ϕ pf l − Cϕ susp (ρ u w) = qr,i , i = 1, 2 (77) p )w] + ∂t i ∂x i i ∂ susp ∂ (ϕ p Cw) + [(u1p + u2p )w] = −(q1 + q2 )w ∂t ∂x

(78)

∂(σtr w) = (q1 + q2 )w ∂t

(79)

Here qr,i is the flux of fluid ‘i’ from reservoir into the fracture (we assume that reservoir is filled with reservoir fluid only, so that qr,1 = 0 and qr,2 = qr ); qi are colmatation rates of particles transported in the fluid ‘i’, and uip are velocities of suspended particles. Averaging is carried out under the assumption that the fracture aperture is a slowly-varying function of the coordinate dw/dx  1 and all parameters of the filtration are uniform over the fracture aperture, the effect of non-uniform fracture height is ignored at this stage as it is stated above. A note on the implementation of the geomechanics models formulated in sections 2.1.3 and 2.3 into the present multiphase flow model of flowback. In Eqs. (77) – (79), we incrorporate the above-shown geomechanics models as follows: proppant embedment reduces the fracture aperture w according to Eq. (18), while proppant compaction reduces proppant pack porosity φ0p and aperture w according to Eqs. (69), (70) and (74). We assume that while the proppant pack is compacted due to the confinement stress, and proppant grains embed into the rock at the fracture face, the fracture aperture w is reduced by a compaction parameter ε33 and the embedment parameter e0 determined according to Eqs. (18) and (74) as follows:

Two-phase filtration inside the hydraulic fracture is approximated using the standard black-oil model modified due to the presence of suspended fines. The particles are transported with the filtrating fluid through the proppant pack and colmatate the pore space, which leads to a decrease in the proppant pack porosity and permeability. Suspension flow through the proppant pack is described using the modified deep-bed filtration model and two-fluid approach presented in [72, 73]. It is assumed that the suspension consists of the carrier fluid (a mixture of fracturing fluid and reservoir fluid), suspended fines (spherical particles of a certain diameter d s ) with the volume concentration C (calculated with respect to the pore volume in a certain elementary volume of a porous medium) and trapped particles with the concentration σtr (calculated with respect to the total elementary volume of a porous medium). Trapped fines modify porosities ϕ susp (available for suspenp sion flow) and ϕ pf l (available only for a clean fluid) of a proppant pack as follows: ϕ0p

   σtr   k s = ϕ0p k0s  0 ϕ pCmax

Here, k0p is the initial permeability of the proppant pack (in the absence of trapped particles); k0s is the permeability of the close pack of fines. Newtonian and compressible carrier fluids are described by the densities ρif , saturations si and velocities uif (index i = 1, 2 denotes the fracturing and reservoir fluid, respectively). Continuity equations averaged over the fracture aperture w are formulated as follows:

3.1. 1D model for fluid-fluid displacement in a fracture

ϕ susp p

3 σtr   , ϕ0pCmax

w = w0 (1 − ε33 ) − 2e0

(75)

(80)

where w0 is an initial fracture opening. Note that according to Eqs. (18) and (74), parameters e0 and ε33 are functions of the confinement stress and, therefore, depend on the reservoir pore pressure pr . Since proppant pack compaction and proppant embedment are inelastic processes, the fracture aperture can only decrease during the simulations (i.e., the corresponding deformations are irreversible). A decrease in the pack porosity φ0p leads to a decrease in the proppant pack permeability k0p according to the following relation [74]:

Here, ϕ0p is the initial porosity of a “clean” proppant pack and Cmax = 0.65 is the maximum concentration of random close pack of spherical particles. The pore space is divided into two types of channels, namely, large pore throats with permeability k p (in between proppant particles) and small pore throats with the permeability k s (in between trapped particles). The permeabilies related to different types of pore throats are defined as follows: 18

ary p(x, t) (e.g., the fluid pressure in the hydraulic fracture): k0p = 0.204rg2 φ4.58 0

(81) qr (x, t) =

Here, rg is the radius of proppant grains. We considered the following expression for colmatation rate validated on experimental data [73, 75]: qi = γCsi e−βσtr rp γ = 4.2589 rg

rp , β = 29.771 rg

!−1.728 (83)

Here γ and β are free parameters, which were obtained by tuning the model against experiments in [73], r p and rg are radius of suspended particles and grains. We used values rg = 0.675 mm, r p = 0.5r pore , were pore radius r pore can be founded using, for example, relation from [76]. Fluid and particle velocities are defined via suspension velocities ui (in total pore space), ui,l (in large porous channels), and ui,s (in small porous channels) as follows: ui = ui,l + ui,s = −

k p k p,i ∂pi k s k s,i ∂pi − 0 µi (C) ∂x µi ∂x

uif = (1 − C)ui , uip = Culi

(87)

Note that the solution (87) is obtained in case of a steady pressure in the hydraulic fracture. At this stage we use this expression to describe unsteady filtration inside the fracture with the time-dependent pressure field p(x, t). Improved solution to the influx from reservoir taking into account unsteady pressure will be implemented in our follow up studies. The system of governing equations (77) – (79) is solved using the finite-difference method and a uniform staggered grid. We used standard IMPES algorithm for splitting the solution of filtration equations. Quasi-linear parabolic pressure equation is solved implicitly by Picard iterations at each time step, until the convergence is reached. Advection equations for fluid saturations and particle volume fraction are solved using the explicit Euler scheme.

(82) !2.2898

kr (p0r − p(x, t)) √ µ02 κtπ

3.2. Numerical results We carried out numerical simulations of fracture flowback taking into account proppant embedment and compaction as well as colmatation of the proppant pack. The aim of the simulations is two-fold: (i) to analyze a relative impact of various conductivity degradation effects on the well production; (ii) to evaluate the effect of well start-up dynamics on the fracture conductivity degradation and long-term well production. Flowback dynamics is governed by the rate of decrease in the bottom-hole pressure due to opening of the choke at the well head. As mentioned above, we assume that at the beginning of the flowback process, the fracture is filled with fracturing liquid (i = 1) and the pressure inside the fracture is equal to the reservoir pore pressure (p0r ). Parameters of the modeling are summarized in Table 5. The values of rock tress components σ011 , f f) f f) σ022 , σ0(e and σ0(e , as well as hydraulic fracturing pressure 11 22 ph involved into expressions describing geomechanics models (39), (69), (70) and (74) are calculated using the expressions (19), (36) and (37). According to the proppant crushing tests, typical volume concentration of crushed ceramic proppant at the confinement stress of 200 atm (∼ 3000 Psi) is around 1%. In the simulations, we intentionally increased initial volume concentration of fines up to 10% to mimic another sources of fines (spalling rock due to proppant embedment, polymer residue, scale precipitates, clay fines). Accurate description of various sources of fines will be incorporated in our future studies. Well production dynamics is evaluated using several parameters: (i) productivity index PI (non-dimensional); (ii) dimensionless fracture conductivity C FD ; (iii) well flow rate Q0 (m3 /d) and (iv) cumulative production V (m3 ) defined as follows   (u1 + u2 ) / p0r − p x=0, t (88) PI =   0 (u1 + u2 ) / pr − p x=0, t=0

(84)

(85)

Here k p,i is the relative permeability of filtration through the proppant pack, k s,i is the relative permeability of filtration through the pack of fines, µi is the viscosity, µ0i is the viscosity of particle-free fluid, pi is a pressure. In what follows we assume that the effect of capillary pressure on the filtration through the proppant pack can be ignored, which is justified by relatively large pores and large filtration velocities. In general, accurate modelling of flowback requires the solution of a coupled problem of filtration inside a proppant-filled hydraulic fracture and in the surrounding reservoir. To simplify the problem, we decouple the problem into two parts. The the influx from the reservoir into the fracture, we will utilize an analytical solution to pressure diffusion equation in a semi-infinite reservoir with the boundary representing a fracture face. The flow is single-phase, constant boundary condition for velocity at the fracture face is considered [77]: √ 2µ02 u0 κt 0 pr (t) = pr − (86) kr π Here kr is reservoir permeability, p0r is the initial reservoir pore pressure, u0 is constant velocity at the fracture face and κ is piezoconductivity coefficient and µ02 is the viscosity of the reservoir fluid. Value of u0 is a free parameter and is chosen to match the total production rate of the fracture with the field data. The reservoir pressure pr obtained from Eq. (86) is substituted into Eq. (19) for hydraulic fracturing pressure, and then the value of ph is used for calculation of parameters ε33 and e0 . Influx from reservoir qr is described using the Carter formula [77, 78], which is a solution to the pressure diffusion equation in a semi-infinite domain with a prescribed pressure at the bound-

C FD = 19

kav wav kr L f rac

(89)

Parameters Initial reservoir and bottom-hole pressure, p0r Final bottom-hole pressure, pw Fracture half-length, L f rac Fracture height, H f rac Initial fracture opening, w0 Initial fracture porosity, ϕ0p Initial volume concentration of fines, C0 Reservoir permeability, kr Poisson ration of the rock, ν Rock cohesion, c Internal friction angle of the rock, φ Depth of the formation, hdepth

Values 50 MPa 30 MPa 137.55 m 46.82 m 3 · 10−3 m 0.3 0.1 10−3 D 0.3 4.89 MPa 62.2o 3 · 103 m

“smooth” BHP dynamics. Therefore, degradation mechanisms (embedment, compaction) manifest themselves differently in these scenarios, which will have different impact on the longterm production. Before presenting the results of numerical simulations of fracture flowback in the well start-up regimes described above, we calculate the rock failure condition (39). It is found that in the “aggressive” well start-up regime, pressure drop ∆p s increases from 10.7 MPa at t = 0 up to 26.5 MPa at t = 1 day, while critical pressure drop ∆pcr s increases from 11.4 MPa up to 20.0 MPa in the same time interval, both parameters are increasing functions of time. Calculations showed that starting from t = 0.2 days, the condition ∆p s > ∆pcr s is met, so that rock failure at the fracture face occurs. This is not the case for the “smooth” well start-up regime, where rock failure criterion is not met due to significantly smaller difference between the pore pressure and the bottom-hole pressure. Results of the numerical simulations are presented in Figs. 39, 40. After one year of production, the fracture conductivity is 25% higher in the “smooth” well start-up regime as compared to that in “aggressive” one. During the initial stage of flowback (well start-up), the fracture conductivity drops faster in the “aggressive” regime because of higher pressure drop between the reservoir pore pressure and the bottom-hole pressure. The dimensionless fracture conductivity C FD decreases rapidly in both flowback scenarios until the time instant, when the embedment depth is equal to the radius of proppant grains, so that further proppant embedment does not occur. During consequent flowback, the value of parameter C FD decreases only slightly due to the proppant pack compaction. Productivity index is 1.7% less in the “smooth” flowback scenario as compared to that in “aggressive” flowback. In the former flowback regime, the average fracture aperture is larger as compared to that of the latter one, while it is vice-versa for the total influx from reservoir, which leads to smaller value of PI for the “smooth” regime. According to results shown in Fig. 40, well flow rate and cumulative well production in the “smooth” regime are 3.1% and 4.2% higher, respectively, as compared to that of “aggressive” regime after 1 year of production. Note that there is a relatively short time interval after the start of production (well start-up), when the well flow rate and cumulative well production is larger in the “aggressive” well start-up regime as compared to “smooth” due to larger difference between the BHP and pore pressure. To conclude this section, we investigate the effect of irreversible elastic deformations of the proppant pack due to compaction and the rock due to proppant embedement on the flowback in different well start-up regimes. We carried out simulations of flowback scenarios shown in Figs. 41, 42, while the proppant pack compaction is described using the elastic model (58). Proppant embedment is described by the model (18), but in contrast to the simulations shown above, we allowed the embedment depth to restore at decreasing confinement stress to mimic the effect of elasticity. The results of flowback simulations using elastic geomechanic models are shown in Figs. 41, 42. In the absence of ir-

Table 5: Input parameters used for modelling of fracture flowback.

Q0 = (u1 + u2 )H f rac w| x=0

(90)

Here kav and wav are fracture permeability and aperture averaged over the fracture length. To examine the impact of geomechanics effects separately, four test cases of flowback process are considered with the following settings of the model: (1) all fracture degradation effects are switched off, (2) only proppant pack compaction is switched on, (3) both proppant pack compaction and proppant embedement are switched on and (4) all fracture degradation effects are switched on (compaction, embedment and colmatation) (Figs. 36, 37). As shown in Fig. 36, proppant pack compaction provides the largest impact on the production rate. Impact of proppant embedment on the fracture conductivity degradation is comparable with that of proppant pack compaction, while the effect is smaller. After 7 days of production, the dimensionless fracture conductivity C FD drops by 58% in the simulation run with only the proppant pack compaction switched on (case (2)) as compared to the case when all fracture degradation mechanisms are switched off (case (1)). In case (3), where additionally proppant embedement is switched on, the reduction in C FD is 8%, while switching on the colmatation results in only 0.8% reduction in C FD . Well flow rate drops by 16% in case (2), by another 13.7% in case (3), and by 0.85% in case (4). Cumulative production decreases by 26% in case (2), by 14.85% case (3), and by 0.87% case (4). Productivity index PI increases by 70% in the second case, by another 83% in the third case, and decreases by another 3% in the forth case. Next we study the effect of bottom-hole pressure (BHP) dynamics on the fracture degradation and long-term well production. The following two flowback scenarios are considered: (1) the “aggressive” well start-up regime, where BHP decreases linearly from initial pressure p0r down to pw during a short period of time (1 day), and (2) the “smooth” regime, where BHP decreases slowly in sync with the piece-wise constant choke opening strategy during longer period of time (7 days), see Fig. 38. The reservoir pressure is assumed equal in both scenarios, whereas the pressure (and, hence, the longitudinal pressure drop) inside the fracture is different for “aggressive” and 20

Figure 36: Dimensionless fracture conductivity C FD and productivity index PI during fracture flowback. Red curve corresponds to case (1) (fracture conductivity degradation effects are switched off), blue curve – case (2) (only proppant pack compaction is switched on), green curve – case (3) (proppant pack compaction and proppant embedment are switched on), purple curve – case (4) (proppant pack compaction, proppant embedment and colmatation are switched on).

Figure 37: Well flow rate Q0 and cumulative production V corresponding to flowback simulation cases (1) (red curve), (2) (blue curve), (3) (green curve) and (4) (purple curve).

reversible deformations, fracture conductivity, productivity index and well production rate depend only on in-situ difference between BHP and the reservoir pore pressure determining the confinement stress (assuming that the rest of the flowback parameters are fixed). Therefore, after the short transient BHP stage, all this parameters are equal, and there is no impact of the well start-up regime on the long-term well production. The larger cumulative well production in the “aggressive” well startup regime as compared to that in “smooth” regime is reached only due to a larger well production rate in the former regime during the short well start-up stage.

by a non-elastic proppant pack compaction (23% decrease in well cumulative production after 1 year of flowback). Proppant embedment provides comparable impact (additional 15.6% decrease in cumulative production as compared to the pure compaction mechanism only). Note that the effect of latter mechanism on fracture degradation is expected to increase with a decrease in the initial fracture aperture and increase in the diameter of proppant grains. Current modelling workflow allowed us to predict very small effect of proppant colmatation on the well production (around 1% reduction in cumulative well production). Regarding proppant embedment, we would like to note that a typical stress-strain diagram for reservoir rocks has a strain softening interval, which occurs after the material strength threshold is achieved [79]. It was first proposed by [32] to take into account strain softening effect in simulations by specifying the dependence of the internal friction angle on the accumulated plastic deformation. Later on, this suggestion was implemented in the Flac3d numerical code for geomechanics problems [69]. In [31], this approach was generalized for the entire deforma-

4. Discussion 4.1. Interpretation of numerical results Numerical simulations of fracture flowback showed that the impact of geomechanic mechanisms on the fracture conductivity and the well production are significantly more pronounced as compared to colmatation of the proppant pack. We obtained that the fracture conductivity degradation is dominated 21

counted for in the current flowback model, namely, mobilization of fines due to rock spalling around embedded proppant grains, tensile rock failure at the fracture face due to rapid change in pressure distribution along the fracture, and proppant flowback. Therefore we expect that our modeling underestimates the negative effect of rapid pressure drop in fracture conductivity and long-term well production. Another important result we obtained is the need to consider accurate non-elastic models for geomechanic effects accompanying fracture flowback. Our comparison of flowback simulations carried out in the framework of elastic and non-elastic models for proppant pack compaction and proppant embedment showed that a large and highly-transient pressure drop during well start-up promotes irreversible non-elastic deformations of the proppant and rock. Elastic geomechanic models does not allow to describe long-term degradation of the fracture conductivity due to “aggressive” well start-up, which limits the range of application of simplified flowback models based on elastic description of geomechanic effects to predict well production.

Figure 38: Pressure drop dynamics considered in the simulations of “aggressive” and “smooth” well start-up regimes; red and blue curves correspond to BHP for aggressive and smooth regimes, respectively, and reservoir pressure in the vicinity of fracture face is shown by green curve.

4.2. Implications for the field operations The present model could be used for defining an optimum regime (safe operating envelope) of the well startup operation. We note that our consideration is focused on each individual fracture, so the proposed modelling workflow allows one to model each separate stage in a multistage fracturing completion, thereby evaluating the impact on production of the entire well comprising a number of stages. At the same time, the problem of mutual influence of different fractures on each other in a multistage completion (via the flow in the near-horizontal well and via the stress state of the rock) will be the topic of a future work along the lines of the present study. Multivariant numerical experiments of the model allow choosing the best option for each set of the geological and mechanical properties. The results show that on qualitative level the application of “smooth” flowback regime may lead to an increase in the cumulative oil production, as the controlled modelling-based design of the smooth flowback operation allows one to reduce undesired phenomena which reduce fracture conductivity and hence have negative impact on production. We also note that typically the flow rate is increasing (not decreasing) during pumping, and hence fracture is opening (not closing), so embedment does not occur during pumping steps until the pumping is completely stopped. Bridging or screen-out during pumping (which lead to fracture closure on proppant) are undesired phenomena that result in premature stop of the operation, hence embedment does not typically occur until the end of the proppant placement, when the fracture actually closes on proppant. However, the uncertainties in the initial data do not yet allow us to make solid quantitative evaluations of the effect. The following lab and field tests to support the model will allow one to formulate closure relations and geomechanic parameters required for the precise numerical evaluation. The ultimate measure is reality, so we plan for a field testing campaign to evaluate the research hypothesis that the smooth flowback scenario allows to preserve the conductivity and eventually yield

tion curve and a methodology was proposed to determine the variable internal friction angle. Friction angle should be the function of accumulated plastic deformation [32]. As a result, the Mohr-Coulomb approach for elastoplastic reservoir rocks was generalized to all deformation stages up to failure due to plastic strain localization in the strain softening regime [66, 80]. In the present work, this approach is used for numerical modeling of proppant embedment into the fracture face. It is important to note that the lack of initial data (e.g., distribution of pore pressure in the vicinity of the hydraulic fracture at the start of flowback, initial distribution of fines and their volume fraction) might cause overestimation of the proppant pack compaction effect and underestimation of the colmatation. As it is pointed out above, in the current model the source of fines is described in a simplified manner, in which all suspended fines are generated instantly at the start of flowback and there are no suspended fines produced during flowback process due to increase in confinement stress acting on proppant. On the other hand, the initial and final packings of the proppant pack are not certain to be cubic and hexagonal, respectively. Defining the dilatancy coefficient Λ (63) significantly depends on the arrangement of packed proppant grains, so that additional lab tests are required to define it. Simulations of different flowback scenarios showed that the dynamics of bottom-hole pressure affects significantly the fracture conductivity degradation. Rapid decrease in BHP results in a plastic (irreversible) deformation of both proppant (pack compaction) and rock (proppant embedment), which leads to a decrease in the long-term fracture conductivity and loss in well cumulative production. Also we obtained that the “aggressive” well start-up regime is accompanied by the rock failure at the fracture face, which is not the case for the “smooth” well startup regime. Rapid pressure drop in a hydraulic fracture promotes development of fracture degradation mechanisms, which are not ac22

Figure 39: Dynamics of dimensionless fracture conductivity C FD and productivity index PI during “aggressive” (red curve) and “smooth” (blue curve) well start-up regimes with BHP dynamics shown in Fig. 38.

Figure 40: Well flow rate Q0 and cumulative production V dynamics for aggressive (red curve) and smooth (blue curve) regimes with BHP dynamics shown in Fig. 38.

Figure 41: Dynamics of dimensionless fracture conductivity C FD and productivity index PI during “aggressive” (red curve) and “smooth” (blue curve) well start-up regimes in the framework of elastic geomechanics models with BHP dynamics shown in Fig. 38.

higher cumulative production, as compared to the aggressive scenario, which damages conductivity (for details of the field testing campaign, see [29]). It should be noted that the exe-

cution of the “smooth” regime might require extra equipment on the well, as well as longer exploitation of the present one. The additional oil production due to the using “smooth” regime 23

Figure 42: Well flow rate Q0 and cumulative production V dynamics in the framework of elastic geomechanic models for aggressive (red curve) and smooth (blue curve) regimes with BHP dynamics shown in Fig. 38.

should compensate these costs. Current mathematical model should be coupled with economical considerations to evaluate the economic effect and formulate recommendations for well start-up scenarios to be implemented in field.

It is demonstrated that the “smooth” flowback scenario allows one to keep the non-dimensional conductivity of the fracture higher by 7.5%, production rate – by 1.5% and cumulative production – by 1.9% as compared to the “aggressive” scenario, where irreversible (non-elastic) proppant pack compaction and proppant embedment into the rock have significant impact on fracture conductivity degradation. On the basis of rock failure model developed above we obtained that rapid pressure drop in the hydraulic fracture during “aggressive” well start-up regime promotes rock failure at the fracture face. Therefore, the negative impact of “aggressive” well start-up on the fracture conductivity is likely underestimated. Based on this parametric analysis, a safe operating envelope is proposed for evaluation in the upcoming field testing campaign [29].

5. Conclusions In this work, we developed an integrated modeling workflow for transient flowback from a fractured well with account for coupled fluid mechanics/geomechanics phenomena. The focus is on the impact of flowback dynamics on the conductivity of a single fracture. We employed the bottom-up approach to the complex process of fracture flowback during well startup. We developed three models for the key geomechanics factors that impact fracture conductivity: (i) proppant embedment with account for rock plasticity and creep, (ii) tensile rock failure, and (iii) non-elastic proppant pack compaction taking into account rearrangement and crushing of proppant grains. These geomechanics submodels are then incorporated into the 1D transient model of fluid-fluid displacement in a fracture. On the basis of flowback simulations we carried out comparative analysis of the effect of different fracture conductivity degradation effects on the well productivity. It is found that non-elastic proppant pack compaction provides the most pronounced affect on the fracture conductivity. Decrease in the fracture aperture due to proppant embedment provides a comparable yet smaller impact on the well production, while proppant pack colmatation effect is found to be even vanishingly smaller. Tensile rock failure has a binary impact (yes/no): if it is realized, we assume the fracture is hydraulically disconnected from the well due to fracture pinch-out in the near-wellbore. A parametric analysis is conducted to evaluate the impact of geomechanics on fracture conductivity in two scenarios of flowback operation: an “aggressive” regime (widely used in current field operations) characterized by rapid choke opening and abrupt decrease of the bottom-hole pressure, and a “smooth” scenario, where the choke opening increases slowly in the stepwise constant regime, which keeps the bottom-hole pressure curve within a safe operating envelope.

We demonstrated that taking into account only simplified elastic models for geomechanic phenomena peculiar to flowback (proppant pack compaction and proppant embedment) does not allow one to capture the difference between smooth and aggressive well start-up regimes (because elastic model predicts the conductivity to restore when the sharp pressure wave propagates from the well into the fracture and further into the reservoir and the pressure gradient decreases). The effect of transient bottom-hole pressure on fracture conductivity and well production can be captured only if irreversible (nonelastic) deformations of the rock and proppant pack are taken into account, which are due to plasticity and creep. In terms of the recommendations for future work, we would propose the following. On the geomechanics side, we recommend to consider a more accurate description of the stress strain state of the rock near the fracture face (for embedment and for tensile rock failure), the problem of proppant flowback, fracture closure between proppant pillars with account for plastic deformations and creep of the proppant pillars and the fracture face, pinch out of the proppant-free cavern near the wellbore. On the fluid mechanics side, the phenomena to be addressed include non-Newtonian rheology of the fracturing fluid, yield stress of unbroken gel, accurate prediction of dynamic sources of formation fines (due to rock spalling, proppant crushing, rock failure 24

Ei jkl rigidity coefficients in the Prandtl-Reuss relation (64) (i, j, k, l = 1, 2, 3) F vertical force exerted on a particle being embedded into a plastic medium F MC strength function determining the Mohr-Coulomb rock failure condition G Shear modulus Gi, j elasticity parameters of orthotropic medium (see Eq. (25), i, j = 1, 2, 3) hdepth depth of the rock formation H f rac fracture height kr permeability of the rock kp permeability of the proppant pack k p,i relative permeabilities of the proppant pack (i = 1, 2) k0p initial permeability of the proppant pack (in the absence of fines) ks permeability of the “small” pores formed with trapped fines k s,i relative permeabilities of the pack of trapped fines (i = 1, 2) k0s permeability of the close pack of fines kav fracture permeability averaged over its length K bulk modulus Ks bulk modulus of the rock material Kf bulk modulus of the fluid K∗ modified bulk modulus of a rock filled with compressible fluid L spacing between proppant grains filling the hydraulic fracture (see Fig. 13) L f rac fracture half-length M Biot modulus p pressure pr pore pressure p0r initial pore pressure pw well (bottom-hole) pressure pf fluid pressure ph hydraulic fracturing pressure, which is a horizontal stress applied from the rock surface to the proppantfilled hydraulic fracture P applied force Pu limiting traction PI productivity index ∆p s pressure drop inside a proppant-filled hydraulic fracture ∆pini initial pressure drop inside a proppant-filled hydraulic s fracture ∆pcr critical pressure drop inside a proppant-filled hys draulic fracture triggering rock failure qr,i volumetric flow rate (influx) of fluids from reservoir into the fracture (i = 1, 2) qi colmatation rate of fines transported in the fluids (i = 1, 2) Q0 well flow rate r radius of a spherical particle rp radius of fines

and precipitation), and advanced models for inflow from the reservoir (to account for transient pressure inside the fracture). A global optimization problem based on the comprehensive forward modelling of flowback is the ultimate goal, which should provide an optimum choke size sequence for given reservoir, fluids, well and fracture input parameters. Acknowledgements AAO, SAB, IAG, and KIT wish to acknowledge the financial support from the Ministry of Education and Science of the Russian Federation (project No 14.581.21.0027, unique identifier RFMEFI58117X0027). The authors are grateful to R.F. Mutalova for the help with preparing illustrations. This work has been financially supported by LLC ”Gazpromneft Science and Technology Center”. The authors are grateful to M.M. Khasanov, A.A. Pustovskikh, and A.S. Margarit for organizational support of this work; to A.A. Yakovlev, A.V. Shurunov, E.V. Shel, I.G. Faizullin, A.L. Vainshtein, E.F. Saifutdinov, R.P. Uchuev, D.V. Prunov and N.A. Chebykin for stimulating and inspiring discussions on details of potential field implications of this research. Nomenclature Latin Symbols a radius of a contact area (see Figs. 6, 28) b1 , b2 free parameters describing embedment (see Eq. (16)) c cohesion ce f effective cohesion C volume concentration of suspended fines C FD dimensionless fracture conductivity D1 , D2 free parameters describing embedment (see Eq. (15)) e0 depth of embedment of a particle into a plastic medium (Eqs. (18) and (19)) initial depth of embedment eini 0 ep depth of embedment of a particle into elastic medium (Hertz solutions (15) and (16)) E Young modulus Es Young modulus of the rock material Ei elasticity parameters of orthotropic medium (see Eq. (25), i = 1, 2, 3) Ee f f effective elastic modulus of the proppant pack

References [1] J. Adachi, E. Siebrits, A. Peirce, J. Desroches, Computer simulation of hydraulic fractures, Int. J. Rock Mech. Min. Sci. 44 (2007) 739–757. [2] F. G. Strickland, et al., Reasons for production decline in the diatomite, belridge oil field: a rock mechanics view, Journal of petroleum technology 37 (03) (1985) 521–526. [3] B. Robinson, S. Holditch, W. Whitehead, et al., Minimizing damage to a propped fracture by controlled flowback procedures, Journal of petroleum technology 40 (06) (1988) 753–759.

25

rg si j T Tr Ti uif uip ui uli uis V p0 V w w0 wav

σ σs σe f f σ0(e f f ) σ x , σy σOE σE,y σi j σ0i j σtr τ φ ϕ, ϕ0 ϕr ϕp ϕ0p ϕ susp p

radius of proppant grains components of the deviatoric stress tensor (i, j = 1, 2, 3) shear stress intensity time scale of creep rock deformation principal components of the shear stress tensor (i = 1, 2) filtration velocity of the carrier fluids (i = 1, 2, see Eq. (85)) filtration velocity of fines (i = 1, 2, see Eq. (85)) filtration velocity of the suspensions with the carrier fluid i (i = 1, 2) filtration velocity of the suspension through “large” porous channels (i = 1, 2, see Eq. (84)) filtration velocity of the suspension through “small” porous channels (i = 1, 2, see Eq. (84)) pore volume cumulative well production fracture aperture and width of the proppant pack initial fracture aperture and width of the proppant pack fracture aperture averaged over its length

ϕ pf l ∆φ slip ψ

Greek Symbols α Biot coefficient αfr internal friction coefficient β free parameter determining the colmatation rate (see Eq. (79)) γ free parameter determining the colmatation rate (see Eq. (79)) γ shear inelastic shear deformation δi j Kronecker symbols (i, j = 1, 2, 3) ε bulk strain εipl cumulative plastic strain intensity ε0 non-dimensional embedment εini initial non-dimensional embedment εi j components of a strain tensor (i, j = 1, 2, 3) ζ variation of the fluid volume fraction in the unit volume of a porous medium θ, θ0 , θ1 angle (see Figs. 6, 9) κ piezoconductivity coefficient λ lateral pressure coefficient (see Eqs. (28), (29)) λi lateral pressure coefficients (see Eq. (26), i = 1, 2) Λ coefficient of dilatancy (compaction) µi viscosity of the particle-laden fracturing and reservoir fluid, respectively (i = 1, 2) µ0i viscosity of the carrier fracturing and reservoir fluid, respectively (i = 1, 2) ν Poisson ratio νs Poisson ratio of the rock material νi, j elasticity parameters of orthotropic medium (see Eq. (25), i, j = 1, 2, 3) ρ rock density ρi density of the particle-free fracturing and reservoir fluids, respectively (i = 1, 2)

ω ωdmg

mean stress yield strength effective mean stress initial effective mean stress principal components of a stress tensor contact stress vertical component of a contact stress components of the stress tensor (i, j = 1, 2, 3) components of the initial stress tensor (i, j = 1, 2, 3) volume concentration of trapped fines shear stress friction angle angles (see Fig. 9) porosity of the rock porosity of the proppant pack initial porosity of the proppant pack porosity of the proppant pack available for suspension flow porosity of the proppant pack available for fluid flow angle between x-axis and tangential curve to the slip line in Eqs. (8), (9) angle determining the arrangement of proppant in the proppant pack (see Fig. 31) angle calculated clockwise from the x-axis to the main principal stress plane (see Fig. 8) material damage function

[4] A. Reinicke, E. Rybacki, S. Stanchits, E. Huenges, G. Dresen, Hydraulic fracturing stimulation techniques and formation damage mechanisms—implications from laboratory testing of tight sandstone–proppant systems, Chemie der Erde-Geochemistry 70 (2010) 107–117. [5] J. D. Weaver, M. Parker, D. W. van Batenburg, P. D. Nguyen, et al., Fracture-related diagenesis may impact conductivity, SPE Journal 12 (03) (2007) 272–281. [6] R. Duenckel, R. D. Barree, S. Drylie, L. G. O’Connell, K. L. Abney, M. W. Conway, N. Moore, F. Chen, et al., Proppants-what 30 years of study has taught us, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, 2017. [7] B. Robinson, S. Holditch, W. Whitehead, et al., Minimizing damage to a propped fracture by controlled flowback procedures, Journal of petroleum technology 40 (06) (1988) 753–759. [8] J. W. Crafton, D. W. Gunderson, et al., Stimulation flowback management–keeping a good completion good, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, 2007. [9] J. W. Crafton, et al., Modeling flowback behavior or flowback equals” slowback”, in: SPE Shale Gas Production Conference, Society of Petroleum Engineers, 2008. [10] E. Detournay, Mechanics of hydraulic fractures, Annual Review of Fluid Mechanics 48 (2016) 311–339. [11] A. Osiptsov, Fluid mechanics of hydraulic fracturing: a review, J. Petr. Sci. and Eng. 156 (2017) 513–535. [12] J. Wang, D. Elsworth, Role of proppant distribution on the evolution of hydraulic fracture conductivity, Journal of Petroleum Science and Engineering 166 (2018) 249–262. [13] J. Wang, D. Elsworth, M. K. Denison, Propagation, proppant transport and the evolution of transport properties of hydraulic fractures, Journal of Fluid Mechanics 855 (2018) 503–534. [14] H. Zhu, Y.-P. Zhao, Y. Feng, H. Wang, L. Zhang, J. D. McLennan, et al., Modeling of fracture width and conductivity in channel fracturing with nonlinear proppant-pillar deformation, SPE Journal. [15] L. Jia, K. Li, Z. Yan, F. Wan, M. Kaita, et al., A mathematical model for calculating rod-shaped proppant conductivity under the combined effect of compaction and embedment, Journal of Petroleum Science and Engi-

26

neering. [16] J. Zhang, L. Ouyang, D. Zhu, A. Hill, Experimental and numerical studies of reduced fracture conductivity due to proppant embedment in the shale reservoir, Journal of Petroleum Science and Engineering 130 (2015) 37– 45. [17] J. M. Terracina, J. M. Turner, D. H. Collins, S. Spillars, et al., Proppant selection and its effect on the results of fracturing treatments performed in shale formations, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, 2010. [18] X. Ye, N. Tonmukayakul, P. Lord, R. Lebas, et al., Effects of total suspended solids on permeability of proppant pack, in: SPE European Formation Damage Conference & Exhibition, Society of Petroleum Engineers, 2013. [19] L. L. Lacy, A. R. Rickards, S. A. Ali, et al., Embedment and fracture conductivity in soft formations associated with hec, borate and water-based fracture designs, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, 1997. [20] J. D. Weaver, R. D. Rickman, H. Luo, R. Logrhy, et al., A study of proppant formation reactions, in: SPE International Symposium on Oilfield Chemistry, Society of Petroleum Engineers, 2009. [21] N. Rayson, J. Weaver, et al., Improved understanding of proppantformation interactions for sustaining fracture conductivity, in: SPE Saudi Arabia Section Technical Symposium and Exhibition, Society of Petroleum Engineers, 2012. [22] J. A. Ayoub, R. D. Hutchins, F. van der Bas, S. Cobianco, C. N. Emiliani, M. D. Glover, M. Koehler, S. Marino, G. Nitters, D. Norman, et al., New findings in fracture clean-up change common industry perceptions, in: SPE International Symposium and Exhibition on Formation Damage Control, Society of Petroleum Engineers, 2006. [23] J. Frederick, H. Hudson, D. Bilden, et al., The effect of fracture and formation flow variables on proppant pack cleanup: An in-depth study, in: SPE Formation Damage Control Symposium, Society of Petroleum Engineers, 1994. [24] E. May, L. Britt, K. Nolte, et al., The effect of yield stress on fracture fluid cleanup nolte, spe, dowell, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, 1997. [25] M. Balhoff, M. J. Miller, et al., Modeling fracture fluid cleanup in hydraulic fractures, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, 2002. [26] D. I. Potapenko, R. D. Williams, J. Desroches, P. Enkababian, B. Theuveny, D. M. Willberg, K. Moncada, P. Deslandes, N. Wilson, R. Neaton, M. Mikovich, Y. Han, C. G., Securing long-term well productivity of horizontal wells through optimization of postfracturing operations, Vol. SPE-187104-MS, 2017. [27] H. L. Norris III, et al., The use of a transient multiphase simulator to predict and suppress flow instabilities in a horizontal shale oil well, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, 2012. [28] A. Osiptsov, S. Boronin, I. Garagash, K. Tolmacheva, Coupled geomechanics-fluid mechanics phenomena during injection and flowback and their impact on hydraulic fracture conductivity in shales, in: Sixth EAGE Shale Workshop, 2019. [29] A. Osiptsov, A. Vainshtein, S. Boronin, et al., Towards field testing of flowback technology on multistage-fractured horizontal wells, in: SPE Technology conference and Exhibition, Moscow, Russia, SPE 196979, 2019. [30] D. Chen, Z. Ye, Z. Pan, Y. Zhou, J. Zhang, A permeability model for the hydraulic fracture filled with proppant packs under combined effect of compaction and embedment, Journal of Petroleum Science and Engineering 149 (2017) 428–435. [31] I. A. Garagash, N. Dubinya, O. Rusina, S. Tikhotsky, I. Fokin, Estimation of rock strength properties from triaxial test data, Geophysical Research (19) (2018) 57–72. [32] P. A. Vermeer, R. De Borst, Non-associated plasticity for soils, concrete and rock, HERON, 29 (3), 1984. [33] Itasca, Itasca consulting group, inc. flac3d - fast lagrangian analysis of continua in 3 dimensions, ver. 3.1, user’s manual (2006). [34] C. Ming, S. ZHANG, L. Ming, M. Xinfang, Z. Yushi, Z. Tong, L. Ning, L. Sihai, Calculation method of proppant embedment depth in hydraulic fracturing, Petroleum Exploration and Development 45 (1) (2018) 159– 166.

[35] Y. Liu, J. Guo, C. Lu, Experimental analysis of proppant embedment mechanism, Chemistry and Technology of Fuels and Oils 54 (2) (2018) 204–210. [36] W. C. Oliver, G. M. Pharr, An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments, Journal of materials research 7 (6) (1992) 1564–1583. [37] A. Ishlinsky, The axial-symmetrical problem in plasticity and the brinell test, Theoretical Research Translation 2 (1944) 47. [38] R. L. Jackson, H. Ghaednia, S. Pope, A solution of rigid–perfectly plastic deep spherical indentation based on slip-line theory, Tribology Letters 58 (3) (2015) 47. [39] D. Tabor, The hardness of metals, Clarendon Press, Oxford, UK, 1951. [40] J. Alcal´a, D. Esqu´e-de los Ojos, Reassessing spherical indentation: Contact regimes and mechanical property extractions, International Journal of solids and structures 47 (20) (2010) 2714–2732. [41] E. Olsson, P.-L. Larsson, A unified model for the contact behaviour between equal and dissimilar elastic–plastic spherical bodies, International Journal of Solids and Structures 81 (2016) 23–32. [42] F.-P. Gao, N. Wang, B. Zhao, A general slip-line field solution for the ultimate bearing capacity of a pipeline on drained soils, Ocean Engineering 104 (2015) 405–413. [43] K. Li, Y. Gao, Y. Lyu, M. Wang, et al., New mathematical models for calculating proppant embedment and fracture conductivity, SPE Journal 20 (03) (2015) 496–507. [44] A. Dinnik, Rock pressure and the design of supports for a circular shaft, Inzh. rabotnik (7). [45] J. Desroches, T. Bratton, Formation characterization: Well logs, Economides, M. & Nolte, K.(2000). Reservoir Stimulation, Third Edition. Wiley. [46] J. Guo, Y. Liu, et al., Modeling of proppant embedment: elastic deformation and creep deformation, in: SPE International Production and Operations Conference & Exhibition, Society of Petroleum Engineers, 2012. [47] I. FLAC3D, Fast lagrangian analysis of continua in 3 dimensions (2002). [48] S. Nakagawa, Sustainability of hydraulic fracture conductivity in ductile and expanding shales, in: Annual Status Report FY2017, FWP ESD 14084, Lawrence Berkeley National Laboratory, 2017. [49] M. Vincent, Proppant types: A review of proppant characteristics, in: PowerPoint Slides – Advanced Well Stimulation class (PEGN 522), 2009. [50] L. Kachanov, Fundamentals of fracture mechanics (1974). [51] L. Kachanov, Rupture time under creeping conditions, Problems of continuum mechanics.- Ed. by JRM Radok.- Philadelphia: SIAM (1961) 202–218. [52] Y. N. Rabotnov, A mechanism of a long time failure [in russian], Creep problems in structural members (1959) 5–7. [53] Y. N. Rabotnov, Creep rupture, in: Applied mechanics, Springer, 1969, pp. 342–349. [54] Y. N. Rabotnov, Creep problems in structural members. [55] N. S. Bulychev, Mechanics of underground structures [in Russian], 1982. [56] M. D. Zoback, Reservoir geomechanics, Cambridge University Press, 2010. [57] I. A. Garagash, Simulation of geomechanical processes during the exploitation of gas fields in permafrost and hydrate-bearing rocks, Electronic Journal ”Georesources, Geo-Energy, Geopolitics” 2 (06). [58] S. Lehnitsky, Theory of Elasticity of an Anisotropic Body [in Russian], Nauka, Moscow, 1977. [59] M. R. Gillard, O. O. Medvedev, P. R. Hosein, A. Medvedev, F. Pe˜nacorada, E. d’Huteau, et al., A new approach to generating fracture conductivity, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, 2010. [60] M. Chertov, et al., Numerical modeling of failure in poroelastic rocks sensitive to pressure drop rate, in: 51st US Rock Mechanics/Geomechanics Symposium, American Rock Mechanics Association, 2017. [61] E. Detournay, A. Cheng, Fundamentals of poroelasticity, comprehensive rock engineering, vol 2: Analysis and design methods, eds, Brown, ET, Fairhurst, CH and Hoek, E. [62] I. Zinchenko, Y. Moiseev, I. Garagash, N. Baransky, K. Smirnov, E. Kozlov, et al., Enhanced prediction of wellbore stability using seismic-based geomechanical modeling (russian), in: SPE Russian Oil and Gas Technical Conference and Exhibition, Society of Petroleum Engineers, 2006. [63] O. Coussy, Poromechanics. hoboken (2004). [64] M. Kachanov, I. Tsukrov, B. Shafiro, Effective moduli of solids with cavi-

27

[65]

[66]

[67] [68]

[69] [70] [71] [72]

[73]

[74]

[75] [76] [77] [78]

[79] [80]

ties of various shapes, Applied Mechanics Reviews 47 (1S) (1994) S151– S174. I. Garagash, Numerical modeling of a hydraulic fracture in dilatant weak rocks and evolution of proppant porosity, in: Book of Abstracts: 10th Hydraulic Fracturing Summit, 15–18 July 2010, Limassol, Cyprus, 2010. J. W. Rudnicki, J. Rice, Conditions for the localization of deformation in pressure-sensitive dilatant materials, Journal of the Mechanics and Physics of Solids 23 (6) (1975) 371–394. S. Man, Compression and flow behavior of proppants in hydraulically induced fracture, Ph.D. thesis, University of Calgary (2016). W. Zheng, Laboratory and discrete element study of proppant crushing and embedment and their influence on fracture conductivity, Ph.D. thesis, Doctoral dissertation, University of British Columbia (2017). C. Itasca, Pfc3d (particle flow code in 3 dimensions), version 4.0, Minneapolis: ICG 3. A. Reinicke, Mechanical and hydraulic aspects of rock-proppant systems: laboratory experiments and modelling approaches. I. Garagash, On the internal instability of materials with damage plasticity, Proc. Kazach. Acad. Sci 3 (1981) 2–7. S. A. Boronin, A. A. Osiptsov, K. I. Tolmacheva, Multi-fluid model of suspension filtration in a porous medium, Fluid Dynamics 50 (6) (2015) 759–768. K. Tolmacheva, S. Boronin, A. Osiptsov, Formation damage and cleanup in the vicinity of flooding wells: Multi-fluid suspension flow model and calibration on lab data, Journal of Petroleum Science and Engineering. A. A. Osiptsov, Hydraulic fracture conductivity: effects of rod-shaped proppant from lattice-boltzmann simulations and lab tests, Advances in water resources 104 (2017) 293–303. S. Xu, B. Gao, J. E. Saiers, Straining of colloidal particles in saturated porous media, Water Resources Research 42 (12). P. Glover, E. Walker, Grain-size to effective pore-size transformation derived from electrokinetic theory, Geophysics 74 (1) (2008) E17–E29. G. Barenblatt, V. Entov, V. Ryzhik, Theory of fluid flows through natural rocks. E. Carter, Optimum fluid characteristics for fracture extension, in: Howard GC, Fast CR, editors. Drilling and production practices, 1957, pp. 261–270. J. C. Jaeger, N. G. Cook, R. Zimmerman, Fundamentals of rock mechanics, John Wiley & Sons, 2009. I. Garagash, V. Nikolaevskiy, Non-associated laws of plastic flow and localization of deformation, Adv. Mech 12 (1) (1989) 131–183.

28

 



Coupled fluid mechanics/geomechanics modeling workflow for transient flowback is developed; Proppant embedment into fracture face with rock plasticity and creep, tensile rock failure and propant pack compaction with proppant rearrangement and crushing are taken into account Parametric study demonstrates that a smooth flowback scenario is preferable over an aggressive regime to preserve the fracture conductivity and to increase production