Multi-Size Proppant Pumping Schedule of Hydraulic Fracturing: Application to a MP-PIC Model of Unconventional Reservoir for Enhanced Gas Production

Abstract

Slickwater hydraulic fracturing is becoming a prevalent approach to economically recovering shale hydrocarbon. It is very important to understand the proppant’s transport behavior during slickwater hydraulic fracturing treatment for effective creation of a desired propped fracture geometry. The currently available models are either oversimplified or have been performed at limited length scales to avoid high computational requirements. Another limitation is that the currently available hydraulic fracturing simulators are developed using only single-sized proppant particles. Motivated by this, in this work, a computationally efficient, three-dimensional, multiphase particle-in-cell (MP-PIC) model was employed to simulate the multi-size proppant transport in a field-scale geometry using the Eulerian–Lagrangian framework. Instead of tracking each particle, groups of particles (called parcels) are tracked, which allows one to simulate the proppant transport in field-scale geometries at an affordable computational cost. Then, we found from our sensitivity study that pumping schedules significantly affect propped fracture surface area and average fracture conductivity, thereby influencing shale gas production. Motivated by these results, we propose an optimization framework using the MP-PIC model to design the multi-size proppant pumping schedule that maximizes shale gas production from unconventional reservoirs for given fracturing resources.

1. Introduction

Hydraulic fracturing is a well stimulation technique to enhance shale gas production from unconventional reservoirs in the petroleum industry [1,2,3]. During hydraulic fracturing, highly pressurized fracturing fluids are pumped to create fractures which provide highly conductive paths for shale gas to flow from reservoirs to the wellbores. Proppant (e.g., sand, alumina, or zirconia) carried by the injected fracturing fluid will prevent the fracture from closing and maintain the high conductivity of propped fracture during shale gas production [4]. Because of extensive use of hydraulic fracturing over the last decade in the United States, production of shale gas from unconventional reservoirs has increased from 5.7 trillion cubic feet in 2000 to over 26 trillion cubic feet in 2019.

In hydraulic fracturing, proppant transport takes place along the horizontal direction by the fracturing fluid, and along the vertical direction by the gravitational settling. In unconventional reservoirs, slickwater is widely used as the fracturing fluid due to its advantages such as low cost and reduced gel damage [5]. However, due to the low-viscosity of slickwater, the proppant settles quickly and forms a proppant bank, which results in a short propped fracture length. Therefore, it is very important to understand the proppant transport behavior during slickwater hydraulic fracturing treatment to effectively create a desired propped fracture geometry.

Many experimental studies have been conducted to understand the proppant transport process with the aim of achieving improved fracture conductivity. Specifically, slots have been used to investigate the effects of fracturing fluid and fracture geometry on the proppant distribution in a single fracture [6,7,8]. Wang et al. (2003) [9] conducted slot experiments to derive bi-power law correlations for bed load proppant transport during hydraulic fracturing using slickwater. In some of the unconventional reservoirs, naturally present fractures will interact with hydraulic fractures and result in complex fracture geometries. In order to understand the proppant transport in a complex fracture network, many experimental studies have been conducted recently using complex slot configurations [10,11,12]. However, all these experimental studies are performed at lab-scales, and thus, these results cannot be directly applied to field-scale hydraulic fracturing.

In order to simulate field-scale proppant transport, the Eulerian–Lagrangian method and Eulerian–Eulerian method have been widely considered. In the Eulerian–Eulerian method, both the fracturing fluid and proppant are solved in Eulerian grids using continuum governing equations [13,14,15]. These continuum models are developed to quickly solve proppant transport, but they cannot capture the discrete characteristics of proppants. In addition, for multi-size proppant pumping schedules, one has to treat particles with different sizes or materials as separate phases, which may require solving additional continuity and momentum equations [16].

In the Eulerian–Lagrangian method, the Lagrangian method is used to solve proppant transport, and Eulerian grids are used to solve the fracturing fluid flow [17,18,19,20]. Compared to the Eulerian–Eulerian framework, the Eulerian–Lagrangian framework can consider multi-size proppant pumping schedules; the latter is computationally more demanding than the former. Within the Eulerian–Lagrangian framework, two models are widely used for proppant transport simulation: the computational fluid dynamics-discrete element method (CFD-DEM) and the multiphase particle-in-cell (MP-PIC) method. The main difference between these two methods is how to deal with the dispersed phase (i.e., proppant particles). In CFD-DEM models, each individual particle is treated separately, and the velocity and position of every particle are calculated by integrating the forces acting on them over time. While this approach is effective in capturing proppant transport at relatively small scales, it is difficult to apply this method to field-scale simulations because of its extremely high requirement of computation power. In MP-PIC models, the dispersed phase is approximated by many parcels; each parcel contains a group of proppant particles. Because of these characteristics, MP-PIC models are able to capture the important features of proppant transport in field-scale geometries at a greatly reduced computational cost.

While the MP-PIC approach has been widely adopted in other particulate systems [16,21,22,23,24,25], its application to multiphase oil/gas systems has not received much attention. Motivated by recent developments [26,27,28], a three-dimensional MP-PIC model has been employed to simulate the proppant transport behavior in a field-scale geometry. In hydraulic fracturing, it is very important to achieve the required fracture geometry, as it will maximize shale gas production rate from unconventional reservoirs [29]. Recently, several studies have been conducted to compute pumping schedules to achieve an optimal fracture geometry [14,29,30,31,32,33,34,35,36,37,38]. However, these pumping schedules are developed using single-sized proppant particles, and the effect of multi-size proppant pumping schedules on the performance of hydraulic fracturing still remains unclear. In the case of unconventional reservoirs, shale gas production depends on propped fracture surface area (PFSA) and average fracture conductivity (AFC). Therefore, it is very important to achieve desired PFSA and AFC values by manipulating the pumping schedule (i.e., fracturing fluid flow rate, proppant concentration, and proppant size injected at the wellbore). In this regard, we use a MP-PIC model within an optimization framework to compute the optimal multi-size proppant pumping schedule that maximizes shale gas production from unconventional reservoirs for given fracturing resources (i.e., proppant amount, etc.).

This paper is organized as follows: In Section 2, we describe a MP-PIC model. In Section 3, we apply the MP-PIC model to a field-scale fracture geometry and demonstrate the importance of considering multi-size proppant pumping schedules. In Section 4, we apply the MP-PIC model within an optimization framework to determine the multi-size proppant pumping schedule that maximizes shale gas production from unconventional reservoirs for given fracturing resources. In Section 5, we present several simulation results to show increased productivity obtained by the computed pumping schedule.

2. High-Fidelity Model Formulation Using the MP-PIC Model

If the MP-PIC method is applied to hydraulic fracturing, mass and momentum equations are solved to describe the fracturing fluid behavior, treating it as a continuous phase, and the Lagrangian description is used to solve Liouville equation to describe the proppant transport behavior. Further, instead of tracking each particle like in CFD-DEM, a group of particles (called a parcel) is tracked, and the particles within each parcel are presumed to have an identical density, volume, velocity, and position (Figure 1). The governing equations for the fluid and particle phases are described below.

2.1. Fluid Phase

The following mass and momentum equations are used to describe the fluid phase.

$$\frac{\partial \left({\alpha }_f{\rho }_f\right)}{\partial t}+\nabla ·\left({\alpha }_f{\rho }_f{u}_f\right)=0$$

(1)

$$\frac{\partial \left({\alpha }_f{\rho }_f{u}_f\right)}{\partial t}+\nabla ·\left({\alpha }_f{\rho }_f{u}_f{u}_f\right)=-\nabla p-F_{fp}+{\alpha }_f{\rho }_{f}g+\nabla ·{\alpha }_f{\tau }_f$$

(2)

where

$$\begin{array}{c}\hfill S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\end{array}$$

(3a)

$$\begin{array}{c}\hfill \varphi (x,u_p,{\rho }_p,V_p,t)\end{array}$$

(3b)

$$\begin{array}{c}\hfill {\alpha }_p=\int \int \int \varphi V_{p}d{V}_{p}d{\rho }_{p}d{u}_p\end{array}$$

(3c)

$$\begin{array}{c}\hfill {\alpha }_f=1-\int \int \int \varphi V_{p}d{V}_{p}d{\rho }_{p}d{u}_p\end{array}$$

(3d)

$$\begin{array}{c}\hfill {\tau }_{ij}=2{\mu }_f{S}_{ij}-\frac{2}{3}{\mu }_f{\delta }_{ij}\frac{\partial u_i}{\partial x_j}\end{array}$$

(3e)

where subscripts f and p are used to denote the fluid and particle phases, respectively; $u_f$ is the fluid velocity; ∇ is the vector differential operator; ${\rho }_f$ is the fluid density; $F_{fp}$ is the rate of momentum exchange between the fluid phase and particle phase; p is the fluid pressure; g is the gravitational acceleration constant; ${\alpha }_p$ and ${\alpha }_f$ are the proppant particle and fluid volume fractions, respectively; ${\tau }_{ij}$ is the fluid phase stress tensor; and ${\mu }_f$ is the fluid viscosity [26,27].

2.2. Particle Phase

Particle distribution function, $\varphi (x,u_p,{\rho }_p,V_p,t)$, is used to describe the dynamics of the particle phase, where x is the particle position, $u_p$ is the particle velocity, ${\rho }_p$ is the particle density, $V_p$ is the particle volume, and t is the current time. The time evolution of $\varphi$ is obtained by solving the Liouville equation while simultaneously considering drag force, pressure gradient, gravity, and inter-particle stress, which is given by [21]:

$$\frac{\partial \varphi }{\partial t}+\nabla ·\left(\varphi u_p\right)+{\nabla }_{u_p}·\left(\varphi A\right)=0$$

(4)

where

$$A=\frac{d{u}_p}{dt}=D_p(u_f-u_p)-\frac{1}{{\rho }_p}\nabla p+g-\frac{1}{{\alpha }_p{\rho }_p}\nabla {\tau }_p$$

(5)

where A is the proppant parcel acceleration, $D_p$ is the drag coefficient, and ${\tau }_p$ is the inter-particle stress. On the right hand side in Equation (5), the first term is the interface drag force, the second term is the dynamic pressure gradient, the third term is the acceleration due to gravity, and the fourth term is the acceleration due to inter-particle stress. $D_p$ is determined by using the drag force model developed by [39], which is given below:

$$D_p=C_d\frac{3}{8}\frac{{\rho }_f}{{\rho }_p}\frac{\left|u_f-u_p\right|}{R_p}$$

(6)

where

$$\begin{array}{c}\hfill C_d=\left\{\begin{array}{cc}\frac{24}{Re}\left({\alpha }_f^{-2.65}+\frac{R{e}^{2/3}}{6}{\alpha }_f^{-1.78}\right)\hfill & Re<1000\hfill \\ 0.44{\alpha }_f^{-2.65}\hfill & Re>1000\hfill \end{array}\right.\end{array}$$

(7a)

$$\begin{array}{c}\hfill Re=\frac{2{\rho }_f\left|u_f-u_p\right|r_p}{{\mu }_f}\end{array}$$

(7b)

$$\begin{array}{c}\hfill r_p={\left(\frac{3V_p}{4\pi }\right)}^{1/3}\end{array}$$

(7c)

where $C_d$ is the dimensionless drag force coefficient, $Re$ is the Reynolds number, and $r_p$ is the proppant particle radius. The inter-particle stress for particle-particle collision is calculated using the following equation [40]:

$${\tau }_p=\frac{P_s{\alpha }_p^{\beta }}{{\alpha }_{cp}-{\alpha }_p}$$

(8)

where $\beta$ is a constant between 2 and 5, $P_s$ is another constant, and ${\alpha }_{cp}$ is the proppant parcel volume fraction at closed packing. Mao et al. (2020) [27] conducted slot experiments to validate the MP-PIC model. The values of $\beta$ and $P_s$ were chosen as 2 and 4 Pa, respectively, to achieve a good fit between the experimental data and MP-PIC model. More details regarding the MP-PIC model validation can be found in Mao et al. (2020) [27].

2.3. Propped Fracture Conductivity

We can obtain the spatio-temporal profile of proppant concentration, $c(x,y,t)={\alpha }_p(x,y,t)$, by solving Equations (1)–(8). However, in order to measure the permeability of a fracture, we need to determine the fracture conductivity, which is defined as the product of proppant bank permeability and fracture width. Fracture permeability, k, is determined by the Kozeny–Carman relation [41], which is given below:

$$k=\frac{{\alpha }_f^3}{C{S}^2{(1-{\alpha }_f)}^2}$$

(9)

where C is the Kozeny–Carman constant which is taken as 5 for unconsolidated sand pack, and S is the proppant particle surface per unit volume of proppant bank. Fracture conductivity, $C_{fd}$, is determined for a given fracture width $W_f$ and S, using the following equation:

$$C_{fd}=1.1\times {10}^{10}\frac{W_f}{S^2}\frac{{\alpha }_f^3}{{(1-{\alpha }_f)}^2}$$

(10)

The coefficient in Equation (10) is the unit conversion factor between the fracture conductivity, $C_{fd}$, which is in md·ft; and fracture width ($W_f$) and the proppant particle surface per unit volume of proppant bank, S, both of which are in inches.

2.4. Reservoir Simulator

The reservoir simulation software from Computer Modeling Group Ltd. (CMG, Calgary, AB, Canada) is used to simulate shale gas production from an unconventional reservoir. The spatial profiles of fracture conductivity and the propped fracture geometry (i.e., the propped average fracture width, propped fracture length, and propped fracture height) obtained from the MP-PIC model at the end of hydraulic fracturing are fed as inputs to CMG software to simulate shale gas production from a well in a unconventional reservoir using a single fracture. A flowchart describing the computational procedures between the MP-PIC model and CMG is given in Figure 2.

Remark 1.

Zeng et al. (2019) [28] developed a model by coupling a two-dimensional MP-PIC method with the Perkins–Kern–Nordgren (PKN) fracture propagation model to study the proppant transport process. In our work, a three-dimensional MP-PIC model developed by Mao et al. (2019) [26] was employed to simulate the proppant transport behavior in a field-scale geometry. More details regarding the implementation of the MP-PIC model can be found in Mao et al. (2019) [26]. We use the MP-PIC model as a virtual experiment to generate data and use the data to develop a mathematical formulation to find the pumping schedule that maximizes the shale gas production from a shale gas well. The major contribution of this work is not developing the MP-PIC model but using the model to develop a multi-size proppant pumping schedule technique for optimal operation of a shale gas well.

Remark 2.

In the MP-PIC model, we considered several particles of identical density, volume, velocity, and specific position as a computational parcel. The particle motion is governed by Newton’s second law. The forces affecting the particles’ motion are fluid-particle drag, gravity, dynamic pressure gradient, and particle stress. When calculating these terms, we consider the velocity and acceleration of the particulate phase. The particle size will play an important role when we calculate the fluid particle drag because the drag coefficient is related to the particle size. Additionally, when we calculate the particle volume fraction, we need to sum up all the physical particle volumes within the Eulerian cell.

Remark 3.

Mao et al. (2020) [27] conducted slot experiments to validate the MP-PIC model. The predicted proppant concentration and proppant bank height profiles were very close to experimental data, which verified that the theory was implemented correctly. More details regarding the MP-PIC model validation can be found in Mao et al. (2020) [27].

3. Field-Scale Simulation Results Using the MP-PIC Model

In this section, by varying the proppant pumping schedule, we study its effect on the proppant transport behavior in a field-scale fracture geometry using the MP-PIC model. The fluid flow and proppant transport were fully coupled, and we used the open source computational fluid dynamics software, OpenFOAM, to simulate fluid flow and proppant transport. The parameters used in the MP-PIC model are given in

Table 1. Model parameters used for the MP-PIC simulation.

Parameter	Value
$\beta$	2
$P_s$	4 Pa
Kozeny–Carman constant (C)	5
Fracture half length	180 m
Fracture height	30 m
Fracture width	0.00762 m
Proppant particle density (${\rho }_p$)	2648 $\mathrm{kg}/{\mathrm{m}}^3$
Pure fluid density ${\rho }_f$	1000 $\mathrm{kg}/{\mathrm{m}}^3$
Particle volume fraction at close packing (${\alpha }_{cp}$)	$0.65$
Parcels injected per second	2500

. Since the fracture is symmetrical, we only modeled a half-fracture using the computation cell configuration shown in Figure 3. The considered fracture geometry is comparable to a field-scale; in particular, a fracture of 180 m in half-length, 30 m in height, and 0.00762 m in width is considered. In the fracture width direction, we considered only one cell because of the high aspect ratio (length over width) of the fracture. The fracturing fluid is injected at the beginning of the half fracture (x = 0) and through the center along the vertical direction (y = 13.5 m to 16.5 m).

The MP-PIC model was simulated using the Ada and Terra clusters provided by Texas A&M high performance research computing facility. The MP-PIC method required around 96 h of calculation via 28 CPUs to simulate proppant transport in the described field-scale fracture geometry. In the literature, available CFD-DEM simulations only considered the fracture length to be up to 10 m, which is much smaller than field-scale fracture geometries; and fracturing fluid injection time to be up to 200 s, which is also not a field-scale injection time [19]. Because of the gain in computational efficiency by the use of the MP-PIC method, we were able to simulate proppant transport in the field-scale fracture geometry.

Then, we studied the effect of proppant size on proppant transport using the fracture configuration shown in Figure 3. For this analysis, we used two different proppant sizes: 100 mesh (150 μm) and $40/70$ mesh (320 μm). This study was performed through three field-scale case studies (denoted as A, B, and C) using three different pumping schedules with a fixed fracturing fluid flow rate of 10 bbl/min and total proppant amount of $M_{prop}$ = 31,500 lb, which are presented in Figure 4. In Case A, we considered the pumping schedule in which the mass percentage of 150 μm proppant is 30%. In Case B and Case C, we varied the mass percentage of 150 μm proppant to 50% and 70%, respectively, by adjusting the injection time of each pumping stage. The overall operation of hydraulic fracturing is for 1 h, and the proppant density is 2650 kg/m${}^3$.

These three pumping schedules were applied to the MP-PIC model. The final proppant concentration distributions are presented in Mao et al. (2020) [27], and it is clear that the proppant transport behavior is highly sensitive to the fracturing fluid flow rate, proppant concentration, and proppant size used during injection. More specifically, it was observed that in Case C, the proppant particles travel farther than the other two cases, wherein the mass percentages of smaller size proppant particles are higher. This is because the gravitational settling of smaller proppant particles is less substantial and can be easily carried deeper into the fracture by the fracturing fluid. In Case A, the proppant particles do not travel much due to their large size, which results in significant gravitational settling, and thereby a short proppant bank.

After obtaining the final proppant concentration distribution, we then calculated the fracture conductivity profile across the fracture for these cases using Equations (9) and (10), as shown in Mao et al. (2020) [27]. It was observed that in Case A, the high fracture conductivity region was mainly due to the high mass percentage of large proppant particles (i.e., high-permeability proppant). On the other hand, in Case C, the fracture conductivity was relatively low due to the presence of small proppant particles (i.e., low-permeability proppant).

After obtaining the fracture conductivity profile and the propped fracture geometry, we then simulated shale gas production using the CMG reservoir simulator for these three cases.

Table 2. Propped fracture surface area (PFSA), average fracture conductivity (AFC), and the cumulative shale gas production volume for the three pumping schedules.

Case	AFC (md·ft)	PFSA (m2)	Cumulative Shale Gas Production Volume for 10 Years (MMft3)
Case A	516.96	2282.5	14.87
Case B	398.42	2301.5	15.02
Case C	284.11	2364.5	14.79

summarizes the AFC, PFSA, and the cumulative shale gas production volume obtained using these pumping schedules. Our results indicated that pumping schedules significantly affect PFSA and AFC, and thereby the cumulative shale gas production volume from unconventional reservoirs. It is observed that using a pumping schedule with a high mass percentage of small size proppant particles (Case C) results in a large PFSA because proppants can penetrate deep inside the fracture due to insignificant gravitational settling. On the other hand, using a pumping schedule with a high mass percentage of large size proppant particles (Case A) results in a high AFC due to the high permeability of large proppant particles. The cumulative shale gas production volume at the end of 10 years is the greatest in Case B, wherein equal mass percentages of small and large proppant particles are used.

The above results show that it is very important to achieve desired PFSA and AFC values to maximize the cumulative shale gas extraction volume from unconventional reservoirs through the manipulation of the pumping schedule (i.e., fracturing fluid flow rate, proppant concentration, and proppant size injected at the wellbore). Motivated by this, we propose an optimization framework in the following section to obtain the pumping schedule that maximizes the cumulative shale gas production volume from unconventional reservoirs by achieving desired PFSA and AFC values for given fracturing resources (i.e., proppant amount, etc.).

4. Optimal Operation of Hydraulic Fracturing Using the MP-PIC Model

Before we compute an optimal multi-size proppant pumping schedule using the MP-PIC model, we present the currently available pumping schedule design techniques and their drawbacks.

4.1. Background on Existing Pumping Schedules

A power-law type pumping schedule for proppant concentration with a fixed fracturing fluid flow rate was proposed by Nolte (1986) [42], which is one of the most commonly used techniques in the oil and gas industry. However, Nolte’s pumping schedule has some limitations: (1) it does not consider all the phenomena that happen during hydraulic fracturing, particularly proppant settling due to gravity; (2) it only considers single fracture propagation; (3) the form of pumping schedule is predefined and applied in an open-loop manner to create a single hydraulic fracture, which may result in early termination of hydraulic fracturing when there is a plant-model mismatch.

To deal with the limitations of Nolte’s pumping schedule, Gu and Desroches (2003) [43] developed a pumping schedule design technique using a hydraulic fracturing simulator. Specifically, an iterative method is used to solve an inverse problem to achieve a desired fracture length and distribution of proppant by manipulating the proppant concentration in pumping schedules. Dontsov and Peirce (2014) [44] developed another pumping schedule design technique to improve Nolte’s pumping schedule; however, they neglected the relationship between proppant concentration and fracture propagation. Recently, Siddhamshetty et al. (2019) [37] developed a new model-based control system to compute a fracturing fluid pumping schedule by mitigating the stress-shadow effects in simultaneously-growing fractures to achieve a desired fracture geometry. Siddhamshetty et al. (2020) [45] further extended their previous work to maximize the total fracture surface area in naturally fractured unconventional reservoirs. However, all of these pumping schedules are developed using single-sized proppant particles, and the effect of multi-size proppant pumping schedules on the performance of hydraulic fracturing was not investigated. In summary, while it is very important to achieve desired PFSA and AFC values to maximize shale gas extraction from unconventional reservoirs, as shown in Section 3, they were not considered in these pumping schedule design techniques.

To address this limitation, we have formulated an optimization problem in the following subsection. The MP-PIC method presented in Section 2 requires running 28 CPUs for 96 h to simulate proppant transport in a field-scale fracture geometry. Because of the large computational requirement, although it is much more affordable than other CFD-DEM models, it cannot be directly used in this optimization problem for the computation of an optimal multi-size proppant pumping schedule. After hydraulic fracturing, we used CMG reservoir simulator to predict shale gas production at the end of 10 years, which again requires huge computational resources. Therefore, it is necessary to develop computationally efficient models to approximate the MP-PIC model and CMG reservoir simulator, which is presented in the following subsection.

4.2. Handling the Computational Requirement of the MP-PIC Model

In this subsection, we developed two computationally efficient models which will be used in the optimization framework to reduce the computational requirements of the MP-PIC model and CMG simulator (Figure 5). The first model is a multi-variable output error state-space (MOESP)-based reduced-order model (ROM) that describes the relationships between the manipulated input variables (i.e., proppant concentration, and proppant size injected at the wellbore) and the output variables (i.e., PFSA and AFC) of the MP-PIC model. In the following subsections, we first describe how this ROM is developed by utilizing data generated from the MP-PIC model.

4.2.1. ROM to Simulate the Hydraulic Fracturing Process

In this subsection, we develop a ROM using MOESP algorithm to describe the proppant transport phenomena during hydraulic fracturing, which is given below:

$$\begin{array}{c}x\left(t_{k+1}\right)=Ax\left(t_k\right)+Bu\left(t_k\right)\hfill \end{array}$$

(11a)

$$\begin{array}{c}y\left(t_k\right)=Cx\left(t_k\right)\hfill \end{array}$$

(11b)

where $u\left(t_k\right)={[C_{x0}\left(t_k\right),d_p\left(t_k\right)]}^{\mathrm{T}}$ are the input variables, where $C_{x0}\left(t_k\right)$ and $d_p\left(t_k\right)$ are the proppant concentration and proppant size injected at the wellbore, respectively; $y\left(t_k\right)=[PFSA\left(t_k\right),AFC\left(t_k\right)]$ are the output variables; and $x\left(t_k\right)$ are the ROM states. We have used the data generated from the MP-PIC model to determine the model parameters A, B, and C, and the initial state estimate, $x\left(0\right)$.

We carried out open-loop simulations by applying a training input (Figure 6) to the MP-PIC model. Utilizing these simulation data, we obtained an ROM, which is a 3rd-order, linear, time-invariant, state-space model. The training input is chosen by considering the minimum and maximum allowable fracturing proppant concentration and proppant size. We have used two different proppant sizes, 100 mesh (150 μm) and 40/70 mesh (320 μm), and considered a non-decreasing trend with respect to the injected proppant concentration, as shown in Figure 6. We observed in Figure 7 that PFSA and AFC obtained from the ROM are close to the values obtained from the MP-PIC model.

We then validated the ROM by using different pumping schedules (Figure 8). The performance of the ROM for the validation pumping schedules is shown in Figure 9, where it is observable that PFSA and AFC obtained from the ROM are close to the true values from the MP-PIC model.

4.2.2. Map to Correlate the Fracture Properties to the Cumulative Gas Production Volume

In this work, we have used CMG reservoir simulator to compute cumulative shale gas production volumes using final fracture properties and geometry at the end of hydraulic fracturing process. Specifically, we have used a package called GEM available in CMG simulator to model shale gas production by providing the component properties and the rock-fluid properties of a shale reservoir. The inputs to CMG-GEM are AFC and PFSA. We used a variety of AFC and PFSA values as inputs to CMG-GEM to generate the corresponding shale gas production profiles. Using the input (PFSA and AFC) and output (cumulative shale production at the end of 10 years) data, we developed a map (Figure 10) where the x-axis is AFC, y-axis is PFSA, and z-axis is the cumulative shale gas production volume at the end of 10 years. There is a general trend that shale gas production increases with PFSA and AFC.

4.3. Optimization Framework to Maximize Cumulative Shale Gas Production Volume from Unconventional Reservoirs

In this subsection, we propose an optimization framework to compute the optimal multi-size proppant pumping schedule that maximizes the cumulative shale gas production volume from unconventional reservoirs for given fracturing resources (i.e., proppant amount, etc.). We have utilized the MOESP-based ROM and map within the following optimization problem to compute the optimal pumping schedule:

$$\begin{array}{cc}\hfill \underset{\begin{array}{c}{C}_{stage,1},\dots ,C_{stage,8}\\ d_{p,stage,1},\dots ,d_{p,stage,8}\end{array}}{max}& Ga{s}_f(AFC\left(t_f\right),PFSA\left(t_f\right))\hfill \end{array}$$

(12a)

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& \mathrm{MOESP}-\mathrm{based}\mathrm{ROM}\left(\mathrm{Equation}\right(11\left)\right)\hfill \end{array}$$

(12b)

$$\begin{array}{cc}\hfill & \mathrm{map}\left(\mathrm{Figure}10\right)\hfill \end{array}$$

(12c)

$$\begin{array}{c}{C}_{stage,m-1}\le C_{stage,m}\le 5\mathrm{PPGA}\hfill \end{array}$$

(12d)

$$\begin{array}{cc}\hfill & d_{p,stage,m}\in \{150\mathsf{\mu }\mathrm{m},320\mathsf{\mu }\mathrm{m}\}\hfill \end{array}$$

(12e)

$$\begin{array}{cc}\hfill & m=1,\dots ,8\hfill \end{array}$$

(12f)

$$\begin{array}{c}\Delta \left({\displaystyle \sum _{m=1}^8}Q{C}_{stage,m}\right)=M_{prop}\hfill \end{array}$$

(12g)

where $\Delta$ is each pumping stage duration; $t_f$ is total hydraulic fracturing operation time; $C_{stage,m}$ and $d_{p,stage,m}$ are the proppant concentration and proppant size injected at the mth pumping stage, respectively; $Ga{s}_f$ is the cumulative shale gas production at the end of 10 years; Q is the fracturing fluid flow rate, which is taken as 10 bbl/min during the entire hydraulic fracturing operation; and $M_{prop}$ is the total amount of proppant available for injection.

In the optimization problem of Equation (12), we used the MOESP-based ROM of hydraulic fracturing process (Equation (12b)) and the map approximating CMG reservoir simulator (Equation (12c)) to reduce the computational requirement. Equations (12d) and (12e) are the constraints on the input variables (e.g., proppant concentration and proppant size injected at the wellbore). The proppant diameters are chosen in non-decreasing order. The units of fracturing fluid flow rate and proppant concentration are bbl/min and PPGA (i.e., by definition, one pound of proppant added to 1 gallon of fracturing fluid), respectively. We considered the total amount of proppant injected during hydraulic fracturing to be equal to 31,500 lb, which is introduced as a constraint using Equation (12g).

5. Proposed Optimization Framework Results

In this section, we present the results to demonstrate the performance of the proposed optimization framework. Please note that in this work we do not consider the fracture propagation, and we only modeled the proppant transport within an half fracture of 180 m in length, 30 m in height, and 0.00762 m in width. Within the optimization framework, we used the MOESP-based ROM and the map to reduce the computational requirement. In the proposed optimization framework, $\Delta$ was chosen to be 450 s with the total operation time of 3600 s. Correspondingly, the pumping schedule has eight pumping stages. The total amount of proppant to be considered is $M_{prop}$ = 31,500 lb.

The pumping schedule obtained using the proposed optimization framework is shown in Figure 11. When we applied this pumping schedule to the MP-PIC model, the AFC and PFSA values were 457.69 md·ft and 2292 ${\mathrm{m}}^2$, respectively. Using these properties at the end of hydraulic fracturing, the cumulative shale gas production volume at the end of 10 years was 15.47 ${\mathrm{MMft}}^3$. The final proppant distribution obtained using the optimal pumping schedule is shown in Figure 12.

The mass percentage of 150 μm proppant in the pumping schedule obtained from the proposed optimization framework was 43%. We compared it with the other cases considered in Section 3. In Case A, we considered the pumping schedule in which the mass percentage of 150 μm proppant was 30%. In Case B and Case C, we varied the mass percentage of 150 μm proppant to 50% and 70%, respectively.

Table 3. AFC, PFSA, and the cumulative shale gas production volume for pumping schedules with different mass percentages of 150 μm and 320 μm proppant sizes.

Case	AFC (md·ft)	PFSA (m2)	Cumulative Shale Gas Production for 10 Years (MMft3)
Case A	516.96	2282.5	14.87
Case B	398.42	2301.5	15.02
Case C	284.11	2364.5	14.79
Optimal	457.69	2292	15.47

shows the AFC, PFSA, and shale gas production volume obtained using these pumping schedules. It was observed that the cumulative shale gas production obtained by using the proposed optimization framework was higher than the other cases.

All the calculations of model training, validation, and deployment in the optimization problem were facilitated using MATLAB on a Dell workstation, powered by Intel(R) Core(TM) i7-4790 [email protected], running the Windows 8 operating system.

6. Conclusions

In this work, we employed a computationally efficient, three-dimensional MP-PIC model to simulate multi-size proppant transport in a field-scale geometry using the Lagrangian description for the proppant particles, and the Eulerian continuum description for the fracturing fluid. In the MP-PIC model, a group of particles (called parcel) are tracked instead of tracking each particle like in CFD-DEM, which reduces the computational cost, and thus allows one to simulate proppant transport in field-scale geometries. We used the MP-PIC model to simulate proppant transport at a field-scale using different multi-size proppant pumping schedules and found that shale gas production from unconventional reservoirs depends on the PFSA and AFC values at the end of hydraulic fracturing. We also found that a smaller size proppant results in a larger PFSA because proppants can penetrate deeper inside the fracture due to less gravitational settling. On the other hand, using large proppant particles results in a high AFC due to their high permeability.

Motivated by these sensitivity studies, we used the MP-PIC model for optimal operation of hydraulic fracturing. To reduce the computational requirement of using the MP-PIC model and CMG within the optimization framework, a MOESP-based ROM was developed to describe proppant transport during hydraulic fracturing, and a map was developed to predict the cumulative shale gas production volume as a function of PFSA and AFC at the end of hydraulic fracturing. Then, we used the MOESP-based ROM and map in the optimization framework to design an optimal multi-size proppant pumping schedule. Simulation results presented in this work show that the AFC and PFSA values obtained by the proposed pumping schedule were able to maximize the cumulative shale gas production volume from an unconventional reservoir for given fracturing resources.