Copyright © the author(s). This work is licensed under a Creative Commons Attribution 4.0 International License. DOI: 10.14800/IOGR.427 Received February 25, 2018; revised March 17, 2018; accepted March 24, 2018. *Corresponding author: shuailiu15@gmail.com 1 A Markov-Chain-Based Method To Characterize Anomalous Diffusion Phenomenon In Unconventional Reservoir Shuai Liu*, Texas A&M University, College Station, USA; Han Li, CGG, Houston, USA; and Peter P. Valkó, Texas A&M University, College Station, USA Abstract The recent success in developing unconventional reservoirs has aroused many new challenges to the theory of reservoir engineering. In this paper, we try to investigate the anomalous diffusion phenomenon caused by the heterogeneity due to the fracture network on the reservoir scale. Firstly, we revisit the physical background of the single-phase flow diffusivity equation by discussing the equivalent single particle diffusion. Combining the characteristics of single particle diffusion with complex fracture geometry, it is indicated that anomalous diffusion phenomenon will be dominant on the reservoir scale, even for single phase production behavior. Then a model based on Markov chain is presented to demonstrate the proposed anomalous diffusion by simulating the particle normal diffusion on a geometric graph and then calculating the relation of the mean square displacement vs. time in the embedding Euclidean space. Based on the simulation results, in consequence, we make discussions on the characteristic size of the heterogeneity due to the fracture network on the reservoir scale, summarize two types of pattern for the anomalous diffusion, and accordingly provide a supportive argument for using the fractional diffusivity equation, in place of the classical one, to model the flow and production behavior in highly fractured unconventional reservoirs. Introduction In the last decade, one of the most prominent progresses for the world’s energy industry was undoubtedly the so-called “shale revolution”. The economic development of unconventional reservoirs, such as the shale gas and tight oil, was made possible by the combination of horizontal drilling and hydraulic fracturing. This practice has been motivating more and more interests in the academic community to investigate the nature of fluid flow and production characteristics in porous media with ultra-low permeability and complex fracture system. In turn, academic progress can shed light on the flow mechanism within unconventional reservoirs and improve the understanding and development of these types of plays. Recent successful development of unconventional reservoirs, especially shales and tight sands with highly complex fracture systems, yields a bunch of new challenges. Many researchers have tried to overcome these challenges by analyzing special phase behavior (Nojabaei et al. 2012; Luo et al. 2016; Jin and Firoozabadi 2016; Luo et al. 2017) and transport phenomena (Riewchotisakul and Akkutlu 2016; Wu and Chen 2016) in nano-size pores, and by considering geomechanics (Kim and Moridis 2012; Yu et al. 2017) such as stress-sensitive permeability. These are all promising avenues, and some informative and inspiring works have been presented or published as has been referred above. On the other hand, despite these progresses, some fundamental concepts and tools for the unconventional reservoir engineering are directly inherited from conventional reservoirs as special situations with ultra-low permeability. Generally speaking, in the unconventional reservoir routine methodology used for modelling flow on the reservoir scale and analyzing production is still along the conventional route, which, for the simplest single-phase case, applies the diffusivity equation with some averaged or upscaled parameters. Mainly based on the mailto:shuailiu15@gmail.com 2 solutions to the single-phase diffusivity equations, Blasingame et al. (1986a, 1986b, 1988, 1989) rigorously derived the rate-decline relation and developed a method to extract the formation properties, such as permeability, drainage area, and hydraulic fracture length by analyzing the flowrate-time data. Due to the formation characteristics of unconventional reservoirs, where the traditional pressure transient analysis (PTA) becomes infeasible, the rate transient analysis (RTA) or production analysis, which has been developed from the above referred works, is the popular way to analyze the reservoir performance. Clarkson and Pedersen (2010) examined the use of classic RTA for analysis of tight oil reservoirs and proposed an integrated RTA approach, which can provide reasonable estimates of hydraulic fracture and reservoir properties. As a typical application of RTA in the shale reservoirs, Belyadi et al. (2015) use this tool to estimate the productivity of Marcellus shale wells and identify the most successful completion methods in the investigated play. To take into account the existence of complex fractures in the formations, the dual- porosity model (Warren and Root 1963; Kazemi 1969; De Swaan 1978) and its variants (Abdassah and Ershaghi 1986; Liu et al. 1987) are widely employed in the study of unconventional reservoirs. Bello and Wattenbarger (2008) combined the slab matrix transient dual porosity model with RTA to study the usage of various shape factor formulations and the effect of matrix geometry on the transient linear response. Fuentes-Cruz and Valko (2015) applied the dual-porosity/dual-permeability model to study the proposed concept of variable matrix-block size and their resulting mathematical model is fundamentally different compared with the standard dual-porosity model due to the interporosity flow. Given the ultra-low permeability of the formation matrix and the extensively elongated transient regime, trilinear model (Lee and Brockenbrough 1983) has been introduced as a simplified but useful “asymptotic” model to study the flow in unconventional reservoirs. Ozkan et al. (2009) applied an analytical trilinear solution to describe the performance of the multiple fractured well and concluded that smaller fracture spacing corresponds to better productivity. Though it hasn’t aroused wide attention in the community of petroleum engineering, anomalous diffusion has been observed in abundant experiments (Adams and Gelhar 1992; Berkowitz and Scher 2001) and well- studied analytically (Metzler et al. 1994; Berkowitz and Scher 1997) and numerically (Zhang et al. 2008; Vilaseca et al. 2011) by hydrologists, biologists, and physicists. Usually, it is the tracers or the small protein molecules that undertake the anomalous diffusion in a flow field that has very high heterogeneity, such as complex fracture systems or massive large molecular as obstacles. However, since the diffusivity equation describes both the tracers’ diffusion and the flow through porous media, by analogy and similarity the flow through porous media is very likely to show characteristics of anomalous diffusion in highly heterogeneous formations, for example shales. Actually, some predecessors in petroleum engineering have done some works on this topic. Raghavan (2011) generalized the concepts of classical diffusion by the fractional derivatives to explain features of anomalous diffusion. Based on fractional derivatives, he described a non- Darcy constitutive equation for the flux. Raghavan and Chen (2017) replaced by the fractional constitutive flux laws the Darcy’s law in their mature dual-porosity model for the multiple fracture horizontal well to obtain the pressure distribution in the drainage region. They provided the asymptotic solutions for the long- term reservoir responses and found that power-law behaviors reflect the heterogeneity in the system. Albinali and Ozkan (2016a) used anomalous diffusion to represent flow in the naturally fractured region between hydraulic fractures for the transient, single-phase production. Based on the sensitivity analyses, they suggested to use this model on a wide range of flow heterogeneity even without the intrinsic details of the formation properties. Albinali and Ozkan (2016b) discussed the basis of anomalous diffusion and combined this new concept with the dual-porosity model to interpret the flow in the heterogeneous formation. Subdiffusion exponents and other coefficients can be extracted from the anomalous diffusion model to help us learn more about the reservoir-rock quality and stimulation efficiency. Holy and Ozkan (2017) recently developed a 1-D numerical model for the linear, single-phase flow undertaking anomalous diffusion. This work provides the foundation for more general multi-dimensional numerical models in the future. The outline of this paper is listed as follows. In the following section, the physical background of the single-phase flow diffusivity equation is revisited by discussing its equivalent form describing individual fluid particle motion. Then, a method is developed by simulating the continuous-time Markov chain (CTMC) on a geometric graph to model the single-phase flow in the fracture network, which is embedded 3 in the formation. Using the simulation results we can calculate the relation of transferred mean square displacement (MSD) vs. time to demonstrate the proposed anomalous diffusion phenomenon. Then by the models described above, the simulation results for a highly fractured formation are given in the Simulated Result section. The plots of MSD with respect to time are displayed to show apparently the characteristics of anomalous diffusion with a non-unit slope. Some further discussions about the simulation results are provided. As a result, the fractional diffusivity equation is suggested to account for the anomalous transport phenomenon in the highly fractured unconventional reservoirs. Revisiting Diffusivity Equation Single-phase flow of slightly compressible fluid in porous media is always the primary and most essential problem for petroleum industry (Dake 1983). Although it may be called as one of the thoroughly investigated topics in this field, we would like to propose a new perspective for this equation when the porous media has a micro-structure in the form of complex fracture network. For simplicity and without losing generality, we will concentrate on the problem in two dimensions. The assumptions for the discussion and the model in this paper are listed as follows. 1. The fluid is of single phase and slight compressibility. The formation volume factor 𝐵 and viscosity 𝜇 can be considered essentially constant. The total compressibility of fracture rock and matrix rock saturated with the fluid are 𝑐𝑡𝑓 and 𝑐𝑡𝑚, respectively. In addition, only isothermal flow is considered. 2. The formation is horizontal, and its thickness ℎ is constant. Its upper and lower boundary is of no flow conditions. 3. The porous media consists of two main parts: the fractures and the matrix. The schematic graph in Figure 1a illustrates the domains of our model. • Although, in reality, the unconventional reservoir formations have the strong heterogeneity due to fractures in nearly all scales (Gale et al. 2014), for the purpose of this work we only take into account the macroscopic fractures (larger than ~10−3 ft). The fracture system consists of natural fracture sets, induced fractures and hydraulic fractures, all of which are interconnected to form a “fracture network” and capable of directly contributing to the flow into the wellbore. The natural fractures and induced fractures have the same formation properties: constant fracture permeability 𝑘𝑓, constant fracture porosity 𝜙𝑓 and constant aperture 𝑏. The hydraulic fractures are assumed to be of planar shape and have infinite conductivity. All the fractures have the same height as the formation thickness and are perpendicular to the horizontal bedding plane. • The matrix is isotropic and homogeneous. This means that we neglect the small-scale fractures, and isolated fractures that have no connections to the interconnected fracture system as mentioned previously. The matrix has constant parameters: permeability 𝑘𝑚, porosity 𝜙𝑚. • Based on the nature of the unconventional reservoir that the fractures have pretty much higher permeability than matrix, it is assumed that the fluid initially located in fractures only transports in the fracture system, and that the fluid initially located in the matrix firstly transports slowly through the matrix into the fracture system, after which it continues transporting only in the fracture system. 4. Based on the above assumptions, the flow in both domains is dictated by Darcy’s law. Due to the pretty small size of aperture 𝑏 compared to ℎ we can assume laminar flow in the fracture system. The hydraulic fractures are of uniform spacings, as shown in Figure 1b, so that each one is located in the center of its own drainage volume. By the geometric and physical symmetry, we regard the symmetrical element, a quarter of the drainage area for one hydraulic fracture, as our problem domain. As shown in Figure 1c, it has the length 𝑙𝑑, the width 𝑤𝑑, and the hydraulic fracture half-length 𝑥𝑓. Furthermore, the drainage area’s outer boundaries are of no flow conditions, except for the part of infinite-conductivity fracture that is at constant pressure condition. For simplicity, we also neglect the flow in fractures across the drainage area boundary, if any. 4 (a) (b) (c) Figure 1—Illustrations of the fracture system and matrix, the drainage area of a hydraulic fracture, and the problem domain. By the above assumptions, the problem has been reduced to a 2-D single phase flow problem with the gravity being neglected. In either the homogeneous and isotropic domains, the matrix or the fracture system, the pressure distribution of slightly compressible flow can be mathematically modeled by the diffusivity equation (Eq. 1). 𝑘𝑖 𝜇 ∆𝑝𝑖 = 𝜙𝑖𝑐𝑡𝑖 𝜕𝑝𝑖 𝜕𝑡 ,…………………………………………………………………………….…...……(1) where 𝑖 = 𝑓,𝑚 and 𝑝𝑖 is the pressure in domain 𝑖. Since the properties are constant in each domain, Eq.1 can be rearranged into the form as shown in Eq. 2. 𝑘𝑖 𝜇𝜙𝑖𝑐𝑡𝑖 ∆𝑝𝑖 = 𝜕𝑝𝑖 𝜕𝑡 ...…………………… ……………………………..…………………….……………(2) As a regular step, we wrap up into a single parameter the parameters on the left-hand-side before the Laplace operator. In Eq. 3 the parameter 𝜂𝑖 is usually called the hydraulic diffusivity coefficient with the unit m2 s⁄ in SI unit. Obviously, 𝜂𝑖 is also a constant parameter in either the fracture system or the matrix. Substituting 𝜂𝑖 back into Eq. 1 gives us a result which has the form of the equation dictating the general diffusion process (Eq. 4), as shown in Eq. 5. 𝜂𝑖 = 𝑘𝑖 𝜇𝜙𝑖𝑐𝑡𝑖 ,……………..……………..………………………..…..…………………………………. (3) 𝐷∆𝐶 = 𝜕𝐶 𝜕𝑡 ,………………………………………………………….…………………………………..(4) Horizontal wellbore Hydraulic fracture Natural fractures & Induced fractures Drainage area wd ld xf 5 𝜂𝑖∆𝑝𝑖 = 𝜕𝑝𝑖 𝜕𝑡 .…………………………………………………….………………………………………(5) Since 𝐷 and 𝜂𝑖 have the same dimension and the above two equations share the same form, there should be some physical analogy between the pressure 𝑝𝑖 and the concentration 𝐶. With the assumption of slightly compressible fluid, the relation between the pressure and density is linear, which means that pressure is equivalent to density in this case. Thus, pressure can be taken as some type of “density” or “concentration” for the single-phase fluid particles. In this perspective, Eq.5 describes the aggregate behavior of the huge number of fluid particles in the porous media. By the theory of normal diffusion (Vlahos 2008), which has been investigated since Einstein (1905)’s work on Brownian motion, the MSD, 〈𝑟2〉, of the fluid particles is related to the time 𝑡 by the diffusivity coefficient 𝐷, as shown in Eq. 6. Eq. 5 can be derived from Eq. 6 (Vlahos 2008), which means Eq. 5 is only valid for the aggregate behavior of those particles whose MSD vs. 𝑡 relationship is dictated by Eq. 6. That is, the validity for using Eq. 5 for a flow on the domain with a specified scale depends on the validity of Eq. 6 for the fluid particles in the same domain. 〈𝑟2〉 ~ 𝑡.………………………………………………………………………………….……………..(6) According to the above analysis, the reason why Eq. 5 can be successfully applied to the conventional reservoir’s flow in multiple scales is that the fluid particles moves approximately in a Euclidean space due to the relatively homogeneous porous media, and that the average motion of the particles is dictated by Eq.6. And for the same reason, Eq. 5 is also valid for the flow in some tight sands with very low permeability but no well-developed fracture networks, except for a quite small 𝜂𝑖. However, this isn’t the case for the unconventional reservoirs with complex fracture systems. By the assumptions, all the produced fluid comes directly from the fracture system. And the high production rate during the early period (several months to years) after the beginning of production all comes from the fluid initially located in the fracture system because of the ultra-low matrix permeability. Thus, the major fluid flow and the production in this period can be readily modeled by solving Eq. 5 with the fracture parameters on the fracture domain, only if the span, shape, and other details of the fracture network are available, which is basically impossible by the current state-of-the-art technology. Consequently, we are forced to model the fracture flow based on the combined domains, since we have much more information and confidence to determine the drainage area. Using the perspective of diffusing fluid particles, it is obvious that the particles only move in the fracture network instead of the full Euclidean space. Therefore, although the particle motion is described by Eq.6 using the 1-D coordinate attached to the fracture, this relation of MSD vs. 𝑡 needs to be transferred to the 2-D coordinate attached to the whole drainage area. This is illustrated in Figure 2. In this figure, a particle diffuses from point 1 to point 2 along the yellow path in the fracture. Its displacement with respect to the fracture coordinate is 𝑟1 + 𝑟2 + 𝑟3 + 𝑟4, while that with respect to drainage area coordinate is 𝑑, which is much smaller than the previous one. Apparently, this transferring should take into account the geometric characteristics of the fracture network, which means the transferred relation of MSD vs. 𝑡 will not follow the linear relation as Eq.6. According to some prior works (Berkowitz and Scher 2001; Vlahos 2008), it can be intuitively proposed that the modified relation should have the power law form as shown in Eq.7. 〈𝑟2〉 ~ 𝑡𝛼.………………………………………………………..…………………..…………………(7) If 𝛼, the diffusivity exponent, doesn’t keep unit, the corresponding process is named as the anomalous diffusion. To demonstrate this phenomenon in highly fractured formation is the main task of the rest of this paper. As a consequence, due to the possible invalidity of Eq.6 when 𝛼 isn’t the unit anymore, using Eq.5 to model the flow and production based on the whole drainage area may be only a very rough simplification and fails to capture some features of the flow through the unconventional reservoir with complex fracture networks. 6 Model Based on Markov Chain To simulate the fluid particle’s motion in the fracture network, which is embedded in a 2-D Euclidean space, we take advantage of its normal diffusion. As has been studied in many publications (Itô 1974; Rogers and Willianms 1994), the fluid particles under normal diffusion can be mathematically modeled to have Markov property. In more details, denoting the location of a single particle at time 𝑡 as 𝑋(𝑡), we have a continuous- time stochastic process {𝑋(𝑡): 𝑡 ≥ 0}, which is considered to have Markov property only if the conditional probability satisfies Eq.8 (Itô 1974; Rogers and Willianms 1994). 𝑃[𝑋(𝑡) = 𝑗|𝑋(𝑠) = 𝑖, 𝑋(𝑡𝑛−1) = 𝑖𝑛−1, 𝑋(𝑡𝑛−2) = 𝑖𝑛−2, ⋯ , 𝑋(𝑡1) = 𝑖1 ] = 𝑃[𝑋(𝑡) = 𝑗|𝑋(𝑠) = 𝑖],.....(8) where 0 ≤ 𝑡1 ≤ 𝑡2 ≤ ⋯ ≤ 𝑡𝑛−2 ≤ 𝑡𝑛−1 ≤ 𝑠 ≤ 𝑡 is any non-decreasing sequence of n+1 times and 𝑖1, 𝑖2, ⋯ , 𝑖𝑛−2, 𝑖𝑛−1, 𝑖, 𝑗 are any n+1 states in the state space of Markov chain. It means that each step of stochastic “jump” only depends on the current states, and the particle acts as it “forgets” the states it has previously experienced. Since in this work the particles only move in the fracture network, the state space of the Markov chain only contains the points belonging to the fractures. For simplicity and the limitedness of our computational resources, we only take the endpoints and the intersection points of the fracture segments as the states, as illustrated in Figure 3. When a particle occupies a state at a given time, it will “jump” after some “waiting time” at the next step to one of the neighbor states. The target of the “jumping” is chosen randomly according to the probability distribution determined by the diffusivity coefficient, the length, and the aperture of all the fracture segments directly connected to the current state, as shown in Eq.9. Figure 2—Illustration of the displacements with respect to different coordinates. Figure 3—Illustration of endpoints and intersection points as states. r1 r2r3 r4 d 1 2 1 2 5 74 6 3 7 𝑃[𝑋(𝑡) = 𝑗|𝑋(𝑠) = 𝑖 ] = { 𝜂𝑗𝑏𝑗 𝑙𝑗⁄ ∑ 𝜂𝑘𝑏𝑘 𝑙𝑘⁄𝑘∈𝑁𝑖 , 𝑗 ≠ 𝑖 0, 𝑗 = 𝑖 ,……………………………………………………..(9) where 𝑃[𝑋(𝑡) = 𝑗|𝑋(𝑠) = 𝑖 ] is the probability for jumping from the current state 𝑖 to the neighbor state 𝑗, 𝑁𝑖 is the set of neighbor states of 𝑖, and 𝜂𝑘, 𝑏𝑘 and 𝑙𝑘 are respectively the diffusivity coefficient, aperture and length of the fracture segment connecting 𝑖 and one of its neighbor state 𝑘 . Only if this kind of probability distribution is applied will the resulting expectation of the stochastic process hold consistency with the solution of the diffusivity equation. However, the above probability distribution only applies to the states corresponding to regular points, not to those points located on the hydraulic fracture. As stated in the assumptions, the hydraulic fracture has infinite conductivity, which is translated into total absorbing states for the particle stochastic motion. It means that whenever a particle jumps into these states, it will stay there forever and will not jump to other states again. Apparently, the probability distribution of these absorbing states is 𝑃[𝑋(𝑡) = 𝑗|𝑋(𝑠) = 𝑖 ] = { 0, 𝑗 ≠ 𝑖 1, 𝑗 = 𝑖 ……………………..………………………..…………….(10) This definition of absorbing state is consistent with the constant pressure condition on the hydraulic fracture boundary. No matter what type the state is, its probability distribution satisfies the normalization condition that ∑ 𝑃[𝑋(𝑡) = 𝑗|𝑋(𝑠) = 𝑖 ]𝑗∈𝑁𝑖 = 1....………………………………………………..…..…………….(11) Both definitions in Eqs.9 and 10 use only the information of current state 𝑖, which manifests this process as a Markov chain. Besides the spatial increment, the temporal increment (the “waiting time”) should also keep the Markov property, which means that only the “memoryless” Poisson distribution should be applied to the temporal increment. While the Markov property of time will later be considered implicitly by the Kolmogorov forward equation, we state the temporal increment’s Poisson distribution here for completeness. By the above concepts, we have implicitly described the fracture network as a graph object, which has the endpoints and the intersection points of the fracture segments as the nodes and the fracture segments between these nodes as edges. This graph is certainly related to the corresponding fracture network, so the graph is classified to be a geometric graph with Euclidean coordinates as one of the node attributes and length as one of the edge attributes. Besides, another critical attribute for each edge is the weight, which defines the transition rate of the CTMC and is related to the probability distribution previously discussed. The edge weight is naturally defined as 𝜂𝑏 𝑙⁄ . Briefly speaking, introducing the graph definition and terminologies explicitly will make our simulation easier to be implemented by programming. So far, we have reduced our problem to a feasible task that a CTMC is to be simulated on a geometric graph object, whose set of nodes are taken as the finite state space, and the weight, 𝜂𝑏 𝑙⁄ , of the edge between node 𝑖 and 𝑗 as the major part of transition rate 𝑞𝑖𝑗. By the Kolmogorov forward equation (Itô 1974), for the node 𝑖, we have 𝑑𝑃𝑖𝑗(𝑡) 𝑑𝑡 = ∑ 𝑃𝑖𝑘(𝑡)𝑞𝑘𝑗𝑘∈𝑆 ...………………………………….…………………..…………….………(12) In Eq.12, 𝑃𝑖𝑗(𝑡) = 𝑃[𝑋(𝑡) = 𝑗|𝑋(0) = 𝑖] is the probability for the particle occurring at node 𝑗 at time 𝑡 after it starts moving from the initial node 𝑖. 𝑞𝑘𝑗 as the transition rate between node 𝑘 and node 𝑗 has the definition as shown in Eq.13. 𝑞𝑘𝑗 = { 𝜂𝑗𝑏𝑗 (𝑙𝑗𝑎 2)⁄ , 𝑘 ≠ 𝑗 𝑎𝑛𝑑 𝑗 ∈ 𝑁𝑘 −∑ 𝜂𝑗𝑏𝑗 (𝑙𝑗𝑎 2)⁄𝑗∈𝑁𝑘 , 𝑘 = 𝑗 0, 𝑜𝑡ℎ𝑒𝑟𝑠 ,……………………………………………(13) 8 where 𝑎 is only a characteristic length for the system to make the unit consistent. And 𝑆 is the set of all the nodes in the graph. Writing the ordinary differential equation (ODE) Eq.12 into the matrix form, we have 𝑑�⃑� (𝑡) 𝑑𝑡 = �⃑� (𝑡)𝐐,………………………………………………………………..……………………….(14) where �⃑� (𝑡) is a stochastic vector, and 𝐐 = {𝑞𝑖𝑗} is the transition rate matrix. From graph theory, we know (Wikipedia 2017) that for a well-defined weighted graph, each entry of 𝐐 is the opposite of the corresponding entry in the graph’s Laplacian matrix, which can be readily obtained. With the help of matrix exponential, the solution of Eq.14 can be symbolically expressed as �⃑� (𝑡) = �⃑� (0) exp(𝐐𝑡),………………………………………………………………………………..(15) where �⃑� (0) is the particle’s initial distribution among the nodes. If we temporarily assume that the matrix exponential in Eq.15 can be successfully evaluated, we can obtain the probabilities for a particle occurring at each node at any time 𝑡 after it starts to move from node 𝑖. The initial stochastic vector for this case is given by Eq.16. �⃑� (0) = [0, 0,⋯ , 0, 1, 0,⋯ , 0] ,………………………………………………….……..………………(16) where 1 is the 𝑖th entry of this vector. These probabilities are equivalent to the particle distribution at time 𝑡 when massive particles all start from 𝑖 to undertake the stochastic process. On the other hand, since it is a geometric graph with Euclidean coordinates as node attributes, we can calculate the square of distance between node 𝑖 and all other nodes to obtain a “vector of distance square”, 𝑣 𝑑𝑠, as shown in Eq.17. 𝑣 𝑑𝑠 = [𝑑1,𝑖 2 , 𝑑2,𝑖 2 , ⋯ , 𝑑𝑖−1,𝑖 2 , 0, 𝑑𝑖+1,𝑖 2 , ⋯ , 𝑑𝑛,𝑖 2 ] ,………………………………………………………….(17) where 𝑑𝑗,𝑖 is the Euclidean distant between node 𝑖 and 𝑗, 𝑑𝑗,𝑖 = [(𝑥𝑖 − 𝑥𝑗) 2 + (𝑦𝑖 − 𝑦𝑗) 2 ] 1 2⁄ . Recalling that our main task is to calculate the relation of transferred MSD vs. 𝑡 for the formation embedding fracture network (or the Euclidean space embedding geometric graph), this task can be done by taking the dot product between �⃑� (𝑡) and 𝑣 𝑑𝑠 as shown in Eq.18. 〈𝑟2〉(𝑡) = �⃑� (𝑡) ∙ 𝑣 𝑑𝑠..……………………..…………………………………………….……………(18) Theoretically, we have provided the solution to our problem, but in practice the stable evaluation of the matrix exponential is problematic. Since usually the longest fractures (~102 ft) are several orders of magnitude longer than the shortest ones (~10-3 ft), the value range of the non-zero entries in the matrix 𝐐 is quite large, which makes 𝐐 a very stiff matrix. Some common methods (Padé approximation, Krylov space methods (Moler and Van Loan 2003)) and packages (expokit (Sidje 1998), MatrixExp in Mathematica (Wolfram Research, Inc. 2017)) for evaluating matrix exponential do not work very well for the cases encountered in this paper. So, we had to resort to numerical methods solving the ODE Eq.14 directly (Basically any “stiff” ODE solver can be used). We have used Mathematica’s NDSolve function (Wolfram Research, Inc. 2017), which integrates many “stiff” ODE solvers into its options, to solve Eq.14 numerically for getting �⃑� (𝑡). Besides the last step of solving ODE, some other packages and algorithms have been used to do some pre-processing work. The 2-D discrete fracture network (DFN) is created using the FRACGEN; The DFN is processed by a modified python package (Splichte 2013) based on the Bentley-Ottmann algorithm; The graph object is readily created using a python package, networkx (Hagberg et al. 2008). Simulation Results Using the described model in the last section, we conduct the simulation to the flow in the fracture networks with different scenarios (Table 1): Case 0: uses the regular meshes with uniform grid size in both 𝑥 and 𝑦 axis; Case 1: uses the fracture network with two orthogonal sets which are randomly generated by FRACGEN; Case 2: also uses the fracture network with two orthogonal sets, but the overall network rotates 45°; 9 Case 3: uses the fracture network with two sets having 40° between them, and all the upper domain; Case 4: uses a very complex fracture network, which has 4 sets of fractures and is the sample, “MWX4f”, in FRACGEN. Table 1—Statistics for the fracture network in various Cases. Case Name Max (ft) Min (ft) Mean (ft) Median (ft) Number of fractures Number of nodes Case 0 10 10 10 10 2525 1285 Case 1 92.4 1.70E-03 11 7.6 1969 1753 Case 2 83.8 1.00E-03 10.5 7.5 1940 1690 Case 3 123.7 1.00E-03 13.1 8.5 753 698 Case 4 81.6 1.00E-03 10 5.9 2572 2214 The domain of the first 4 cases all have the dimensions of 500 ft ×250 ft coming from the well spacing 1000 ft and fracture spacing 500 ft. And all cases have a hydraulic fracture half-length of 400 ft (𝐼𝑥 = 0.8) which is represented by the blue line on the boundary. Since each case contains hundreds and even thousands of nodes, the comprehensive investigation is impossible without the help of some statistical methods, which is out of the scope of this paper. So, we only select 10 nodes randomly in each case to do the simulation on them and to show the relation of MSD versus 𝑡. The plots and parameters of the fracture networks, and the results are shown below. Case 0. Figure 4—Map view plot of the fracture network and 10 sampled nodes in Case 0. 10 Node 34 Node 141 Node 155 Node 181 Node 440 Node 839 Node 844 Node 1074 0.10 1 10 100 1000 104 t 1 100 104 MSD, ft2 0.10 1 10 100 1000 104 t 1 100 104 MSD, ft 2 0.10 1 10 100 1000 104 t 1 100 104 MSD, ft2 0.10 1 10 100 1000 104 t 1 100 104 MSD, ft2 0.10 1 10 100 1000 104 t 1 100 104 MSD, ft2 0.10 1 10 100 1000 104 t 1 100 104 MSD, ft 2 0.10 1 10 100 1000 104 t 1 100 104 MSD, ft 2 0.10 1 10 100 1000 104 t 1 100 104 MSD, ft 2 11 Node 1093 Node 1126 Figure 5—log-log plots of MSD vs. time for 10 sampled nodes in Case 0. Case 1. Figure 6—Map view plot of the fracture network and 10 sampled nodes in Case 1. Figure 7—Histogram of the fracture length in Case 1. 0.10 1 10 100 1000 104 t 1 100 104 MSD, ft 2 0.10 1 10 100 1000 104 t 1 100 104 MSD, ft 2 12 Node 72 Node 150 Node 236 Node 357 Node 542 Node 801 Node 964 Node 1104 1 100 104 t 1 100 104 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft2 1 100 104 t 1 100 104 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft2 1 100 104 t 1 100 104 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft 2 13 Node 1526 Node 1662 Figure 8—log-log plots of MSD vs. time for 10 sampled nodes in Case 1. Case 2. Figure 9—Map view plot of the fracture network and 10 sampled nodes in Case 2. Figure 10—Histogram of the fracture length in Case 2. 1 100 104 t 1 100 104 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft 2 14 Node 139 Node 428 Node 731 Node 965 Node 1031 Node 1127 Node 1194 Node 1247 1 100 104 t 1 100 104 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft 2 0.1 100 105 t 0.001 1 1000 MSD, ft2 1 100 104 t 0.1 10 1000 105 MSD, ft 2 1 100 104 t 0.1 10 1000 105 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft2 15 Node 1391 Node 1650 Figure 11—log-log plots of MSD vs. time for 10 sampled nodes in Case 2. Case 3. Figure 12—Map view plot of the fracture network and 10 sampled nodes in Case 3. Figure 13—Histogram of the fracture length in Case 3. 1 100 104 t 1 100 104 MSD, ft 2 0.10 1 10 100 1000 104 t 0.10 1 10 100 1000 104 MSD, ft2 16 Node 7 Node 10 Node 187 Node 225 Node 358 Node 377 Node 471 Node 576 1 100 104 t 1 100 104 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft 2 0.1 100 105 t 0.01 10 104 MSD, ft 2 1 100 104 t 0.1 10 1000 105 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft 2 1 100 104 t 1 100 104 MSD, ft 2 17 Node 678 Node 690 Figure 14—log-log plots of MSD vs. time for 10 sampled nodes in Case 3. Case 4. Figure 15—Map view plot of the fracture network and 10 sampled nodes in Case 4. Figure 16—Histogram of the fracture length in Case 4. 0.10 1 10 100 1000 104 t 1 100 104 MSD, ft 2 0.010 0.100 1 10 t 0.100 10 1000 MSD, ft 2 18 Node 340 Node 488 Node 631 Node 648 Node 1137 Node 1152 Node 1714 Node 1958 0.1 10 1000 105 t 0.1 10 1000 105 MSD, ft 2 0.1 10 1000 105 t 0.1 10 1000 105 MSD, ft 2 0.1 10 1000 105 t 0.1 10 1000 105 MSD, ft 2 0.1 10 1000 105 t 0.1 10 1000 105 MSD, ft 2 0.1 10 1000 105 t 0.1 10 1000 105 MSD, ft 2 0.1 10 1000 105 t 0.1 10 1000 105 MSD, ft 2 0.1 10 1000 105 t 0.1 10 1000 105 MSD, ft 2 0.1 10 1000 105 t 0.1 10 1000 105 MSD, ft 2 19 Node 2128 Node 2182 Figure 17—log-log plots of MSD vs. time for 10 sampled nodes in Case 4. Discussions of the Results When performing the above 5 case studies, for simplicity we assign unit value to the diffusivity coefficient 𝜂𝑗, the aperture 𝑏𝑗, the fracture height ℎ, and the characteristic length 𝑎 and only use 1 𝑙𝑗⁄ as the weight of the edges to construct the weighted graph’s Laplacian matrix. Since by our assumptions 𝜂𝑗, 𝑏𝑗, ℎ and 𝑎 are all constants, only a constant factor is needed to transfer the “time” 𝑡 in the above plots to the real time, and the phenomena displayed aren’t affected by the values of these parameters as far as they fall into the ranges guaranteeing the basic assumptions of this work. In each of the above plots, besides the resulting relations of MSD vs. 𝑡 presented in log-log plot, some dash lines with different colors are added to accommodate our analysis. Firstly, the most straight-forward one is the orange dash line, which has exactly the unit slope. Because Case 0 uses the regular grid to do the simulation, it is exactly the same thing as conducting numerical simulation to the homogeneous reservoirs with the most common spatial discretization. Thus, Case 0’s result should be expected to show the nature of normal diffusion, as we have discussed in the previous section. It is exactly what we see in Case 0’s plots (Figure 5). The unit-slope orange line tracks firmly the resulting curve until deviation happens at the late time, when most of the particles have been “locked” in the absorbing nodes and MSD levels out. Some very small deviations occur at the early or intermediate time of some nodes, Node 155, Node1074, and Node 1126. This is due to the boundary effect of our finite graph on these nodes that are relatively close to the boundaries. Then, using Case 0 as a reference, we can immediately notice the discrepancies between other cases (Figure 8, Figure 11, Figure 14, and Figure 17) and the normal diffusion. Although it seems that every node in from Case 1 to Case 4 has its unique relation due the various local characteristics of the random created DFN, an obvious common phenomenon is that most of the curves deviate from the unit- slope orange in quite early time. Many of them are concave downwards displaying subdiffusion, with some of them being concave upwards corresponding to superdiffusion. The pair of green dash lines in each plot give us some sense of the scale of heterogeneity. Roughly speaking, the average fracture length in all 5 cases is about 10 ft, which corresponds to 100 ft2 for distance square. So, the green line pair shows the point where the MSD reaches a value that can be taken as characteristic size of heterogeneity. In Case 0 where the fracture segments are all 10 ft long, the slope keeps the unit value before and after 100 ft2 MSD is reached due to the totally homogeneous fracture distribution. In other cases, most deviations happen within the range 10 ft2~1000 ft2MSD, or 1 ft ~ 10 ft, again roughly speaking. Moreover, we can see that each plot has the exact unit-slope section in the very early time when MSD is very small. Recall that only fractures longer than 10-3 ft are modeled in our work, and the number of fractures with length less than 10-1 ft is very small from all the four histograms. This means that the flow domain with its characteristic dimensions less than the average fracture length can be considered homogeneous to some extent. However, since only 10 nodes have been sampled out of thousands for each case, this point can only be taken as a reasonable suggestion. It should also be noted that this observation has been made while neglecting all other types of heterogeneity. 0.10 1 10 100 1000 104 t 0.10 1 10 100 1000 104 MSD, ft 2 0.10 1 10 100 1000 104 t 0.10 1 10 100 1000 104 MSD, ft 2 20 Although currently we cannot connect the fracture network characteristics (fracture density, fracture length distribution, fracture sets and clusters, and so on) to the properties of the diffusion due to the limited number of sampled nodes, we are capable to summarize some common diffusion patterns when looking at the results from various DFNs comprehensively. We have already noted the early time concave period. Similarly, the absorbing boundary domination at the late time is quite apparent. From the homogeneous Case 0, it is apparent that the absorbing boundary has the effect of greatly leveling out the curve. Bearing this in mind, we roughly draw the red straight dash line by hand to try to catch up its average deviating trends in early or intermediate times, and to try to rule out the absorbing boundary effects. So, these handmade trendlines are to some extent secant lines of the resulting curves. Consequently, the slope of these trendlines, which is exactly the diffusivity exponent 𝛼 in Eq. (7), tells the types of anomalous diffusion in a given range of scale: subdiffusion for 𝛼 < 1 and superdiffusion for 𝛼 > 1 . Generally speaking, there are two diffusion patterns for the flow in the fracture network, if the anomalous diffusion does happen. Type 1 is that after the early normal diffusion section, the curve directly concaves downward, and the flow begins to undertake subdiffusion until the absorbing boundary dominates. The typical examples for this type are Node 964, 1104 in Case 1 (Figure 8), Node 139, 428, 731, 965, 1650 in Case 2 (Figure 11), Node 10, 471 in Case 3 (Figure 14), and Node 1137, 1714, 2128 in Case4 (Figure 17). Type 2 is characterized with a hump-like shape, which means the curve concaves upward to undertake superdiffusion firstly, and then becomes subdiffusion in a larger MSD range. The typical examples for this type are Node 150, 357, 801 in Case 1 (Figure 8), Node 1031 in Case 2 (Figure 11), Node 7, 377 in Case 3 (Figure 14), and Node 340, 1958, 2182 in Case 4 (Figure 17). Only if the absorbing boundary begins to dominate, the undertaking diffusion is overwhelmingly affected and even totally hidden. These two patterns are kind of intuitive for a flow into hydraulic fractures (the absorbing nodes). We feel that more complex patterns can be expected in some highly heterogeneous formations, possibly alternating subdiffusion and superdiffusion time periods. Based on our prior discussions, we now try to answer our major question: Is it valid enough to apply the classical diffusivity equation on the reservoir scale to model single-phase flow and to analyze the production? In our opinion, the diffusion pattern classification provides a perspective for the answer. Looking back again at Case 0, we make a tangential line (the horizontal black dash line) from the final plateau of the curve, and the tangential point and the intersection point with the unite-slope trendline both correspond to a time (the vertical black dash lines). Not so surprisingly, we find that in the normal diffusion case, the differences between these two times for different nodes are always around 1 log cycle. So, this can serve as a standard to tell the validity of the averaged model (or homogenization) with the classical diffusivity equation. The type 1 pattern, when only subdiffusion happens, obviously enlarges this time difference to around 1 1 2⁄ or even 2 log cycles, and hence disproves the validity. On the other hand, the superdiffusion and subdiffusion in type 2 pattern counteract with each other to some extent and can lead to the time difference around 1 log cycle, such as Node 150 in Case 1 (Figure 8), Node 7, Node 377 in Case 3 (Figure 14), and Node 340 in Case 4 (Figure 17). So, if both type of nodes exist in a situation, we should determine or estimate which type dominates. It would require massive simulation conducted, or some advanced statistical tool employed. However, type 2 pattern can also lead to confusion by substantially shortening the time difference, which displays a superdiffusion on average, like the Node 357, 542, 801 in Case 1 (Figure 8), Node 1031 in Case 2 (Figure 11), Node 2181 in Case 4 (Figure 17). Actually, depending on the relation between superdiffusion, subdiffusion and absorbing boundary effects, type 2 pattern can manifest itself as “apparent” superdiffusion, subdiffusion or normal diffusion on average. The above observations provide further supportive argument for using the fractional diffusivity equation for modeling flow in porous media and production behavior in highly fractured unconventional reservoirs. In the end, we can also make some discussions about the common dual-porosity model in the perspective of anomalous diffusions. Although the dual-porosity model also bears the two flow domain assumptions (fracture and matrix), it implicitly makes further assumption that the flow particles in the fracture domain undertake the diffusion in a Euclidean space. Thus, although it considers matrix supplementing fluid into fracture correctly, it still fails to capture the heterogeneity due to the geometry of the fracture network, which may limit the model’s application in highly fractured unconventional reservoirs. 21 Conclusions In summary, based on our simulation work and the relevant discussions, we can reach the following conclusions: 1. The physical background of the single-phase diffusivity equation is revisited. Combining the characteristics of single particle diffusion with complex fracture geometry, it is indicated that anomalous diffusion phenomenon will be dominant on the reservoir scale, even for single phase production behavior; 2. A Markov-chain-based model is presented to model single particle diffusion and it is demonstrated that anomalous diffusion characteristics emerge from considering normal diffusion on a graph that is embedded in a 2-D Euclidean space; 3. According to the simulation results, two known types of diffusion pattern (subdiffusion and superdiffusion) may rise from the same graph geometry (possibly alternating in time); 4. The fractional diffusivity equation have advantages in characterizing flow and production in highly fractured unconventional reservoirs. Conflicts of Interest The author(s) declare that they have no conflicting interests. Nomenclature 𝛼 = diffusivity exponent 𝜂 = hydraulic diffusivity coefficient, ft2/sec 𝜇 = viscosity, cp 𝜙 = porosity 𝑎 = characteristic length, ft 𝐵 = formation volume factor, RB/STB 𝑏 = fracture aperture, ft 𝐶 = concentration, mol/ft3 𝑐𝑡 = total compressibility, psi-1 𝐷 = diffusivity coefficient, m2/s 𝑑 = euclidean distance between to nodes in the geometric graph, ft ℎ = formation thickness, ft 𝑘 = permeability, mD 𝑙 = fracture length, ft 𝑙𝑑 = length of the problem domain, ft 𝑁𝑖 = the set of the neighbor states of state 𝑖 𝑃 = probability 𝑝 = pressure, psi 𝐐 = transition rate matrix 𝑞 = transition rate, 1/sec 𝑟 = partical displacement, ft 𝑆 = the set of all the nodes in a graph 𝑡 = time, sec 𝑣 𝑑𝑠 = vector of distant square, ft2 𝑤𝑑 = width of the problem domain, ft 𝑋 = location/state of a particle Subscripts 𝑓 = fracture domain 𝑚 = matrix domain 22 References Abdassah, D. and Ershaghi, I. 1986. Triple-Porosity Systems for Representing Naturally Fractured Reservoirs. SPE Form Eval 1(2): 273-280. SPE-13409-PA. Adams, E. E. and Gelhar, L. W. 1992. Field Study of Dispersion in a Heterogeneous Aquifer: 2. Spatial Moments Analysis. Water Resources Research 28(12): 3293-3307. Albinali, A. and Ozkan, E. 2016a. Analytical Modeling of Flow in Highly Disordered, Fractured Nano-Porous Reservoirs. Paper presented at the SPE Western Regional Meeting, Anchorage, Alaska, 23-26 May. SPE-180440- MS. Albinali, A. and Ozkan, E. 2016b. Anomalous Diffusion Approach and Field Application for Fractured Nano-Porous Reservoirs. Paper presented at the SPE Annual Technical Conference and Exhibition, Dubai, UAE, 26-28 September. SPE-181255-MS. Bello, R.O. and Wattenbarger, R.A. 2008. Rate Transient Analysis in Naturally Fractured Shale Gas Reservoirs. Paper presented at the CIPC/SPE Gas Technology Symposium 2008 Joint Conference, Calgary, Alberta, Canada, 16-19 June. SPE-114591-MS. Belyadi, H., Yuyi, S., and Junca-Laplace, J. P. 2015. Production Analysis Using Rate Transient Analysis. Paper presented at the SPE Eastern Regional Meeting, Morgantown, West Virginia, 13-15 October. SPE-177293-MS. Berkowitz, B. and Scher, H. 1997. Anomalous transport in random fracture networks. Physical review letters 79 (20): 4038. Berkowitz, B. and Scher, H. 2001. The Role of Probabilistic Approaches to Transport Theory in Heterogeneous Media. Transport in Porous Media 42(1-2): 241-263. Blasingame, T. A. and Lee, W. J. 1986a. Properties of Homogeneous Reservoirs, Naturally Fractured Reservoirs, and Hydraulically Fractured Reservoirs From Decline Curve Analysis. Paper presented at the Permian Basin Oil and Gas Recovery Conference, Midland, Texas, 13-15 March. SPE-15018-MS. Blasingame, T. A. and Lee, W. J. 1986b. Variable-Rate Reservoir Limits Testing. Paper presented at the Permian Basin Oil and Gas Recovery Conference, Midland, Texas, 13-15 March. SPE-15028-MS. Blasingame, T. A. and Lee, W. J. 1988. The Variable-Rate Reservoir Limits Testing of Gas Wells. Paper presented at the SPE Gas Technology Symposium, Dallas, Texas, 13-15 June. SPE-17708-MS. Blasingame, T. A., Johnston, J. L., and Lee, W. J. 1989. Type-Curve Analysis Using the Pressure Integral Method. Paper presented at the SPE California Regional Meeting, Bakersfield, California, 5-7 April. SPE-18799-MS. Clarkson, C. R. and Pedersen, P. K. 2010. Tight Oil Production Analysis: Adaptation of Existing Rate-Transient Analysis Techniques. Paper presented at the Canadian Unconventional Resources and International Petroleum Conference, Calgary, Alberta, Canada, 19-21 October. SPE-137352-MS. Dake, L. P. 1983. Fundamentals of reservoir engineering, 1st edition. Amsterdam, Holland: Elsevier. De Swaan, A. 1978. Theory of Waterflooding in Fractured Reservoirs. SPE J 18 (2): 117-122. SPE-5892-PA. Einstein, A. 1905. Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Annalen der physik 322(8): 549-560. FRACGEN and NFFLOW version 14.9, A DOE-sponsored project. https://edx.netl.doe.gov/dataset/fracgen-and- nfflow-version-14-9. Accessed December 18, 2017. Fuentes-Cruz, G. and Valko, P. P. 2015. Revisiting the Dual-Porosity/Dual-Permeability Modeling of Unconventional Reservoirs: The Induced-Interporosity Flow Field. SPE J 20(1): 124-141. SPE-173895-PA. Gale, J. F., Laubach, S. E., Olson, J. E., et al. 2014. Natural Fractures in Shale: A Review and New Observations. AAPG bulletin 98(11): 2165-2216. Hagberg, A. A., Schult, D. A., and Swart, P. J. 2008. Exploring Network Structure, Dynamics, and Function using NetworkX. Paper presented at the Proceedings of the 7th Python in Science Conference, Pasadena, CA, 2-5 August. Holy, R. and Ozkan, E. 2017. Numerical Modeling of 1D Anomalous Diffusion in Unconventional Wells Using a NonUniform Mesh. Paper presented at the SPE/AAPG/SEG Unconventional Resources Technology Conference, Austin, Texas, 24-26 July. URTEC-2695593-MS. Itô, K. 1974. Diffusion processes. London, UK: John Wiley and Sons, Inc. Jin, Z. and Firoozabadi, A. 2016. Thermodynamic Modeling of Phase Behavior in Shale Media. SPE J 21(1): 190- 207. SPE-176015-PA. Kazemi, H. 1969. Pressure Transient Analysis of Naturally Fractured Reservoirs with Uniform Fracture Distribution. SPE J 9(4): 451-462. SPE-2156-PA. Kim, J. and Moridis, G. J. 2012. Numerical Studies on Coupled Flow and Geomechanics with the Multiple Porosity Model for Naturally Fractured Tight and Shale Gas Reservoirs. Paper presented at the 46th U.S. Rock Mechanics/Geomechanics Symposium, Chicago, Illinois, 24-27 June. https://edx.netl.doe.gov/dataset/fracgen-and-nfflow-version-14-9 https://edx.netl.doe.gov/dataset/fracgen-and-nfflow-version-14-9 23 Lee, S.T. and Brockenbrough, J. 1983. A New Analytic Solution for Finite Conductivity Vertical Fractures With Real Time and Laplace Space Parameter Estimation. Paper presented at the SPE Annual Technical Conference and Exhibition, San Francisco, California, 5-8 October. SPE-12013-MS. Liu, X., Chen, Z., and Jiang, L. 1987. Exact Solution of Double-Porosity, Double-Permeability Systems Including Wellbore Storage and Skin Effect. Paper presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 27-30 September. SPE-16849-MS. Luo, S., Lutkenhaus, J. L., and Nasrabadi, H. 2016. Confinement-Induced Supercriticality and Phase Equilibria of Hydrocarbons in Nanopores. Langmuir 32(44): 11506-11513. Luo, S., Lutkenhaus, J. L., and Nasrabadi, H. 2017. Multi-Scale Fluid Phase Behavior Simulation in Shale Reservoirs by a Pore-Size-Dependent Equation of State. Paper presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 9-11 October. SPE-187422-MS. Metzler, R., Glöckle, W.G., and Nonnenmacher, T.F. 1994. Fractional Model Equation for Anomalous Diffusion. Physica A: Statistical Mechanics and its Applications 211(1): 13-24. Moler, C. and Van Loan, C. 2003. Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-five Years Later. SIAM review 45(1): 3-49. Nojabaei, B., Johns, R.T., and Chu, L. 2012. Effect of Capillary Pressure on Fluid Density and Phase Behavior in Tight Rocks and Shales. Paper presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 8-10 October. SPE-159258-MS. Ozkan, E., Brown, M.L., Raghavan, R.S., et al. 2009. Comparison of Fractured Horizontal-Well Performance in Conventional and Unconventional Reservoirs. Paper presented at the SPE Western Regional Meeting, San Jose, California, 24-26 March. SPE-121290-MS. Raghavan, R. 2011. Fractional derivatives: application to transient flow. Journal of Petroleum Science and Engineering 80(1): 7-13. Raghavan, R. and Chen, C. 2017. Addressing the Influence of a Heterogeneous Matrix on Well Performance in Fractured Rocks. Transport in Porous Media 117(1), 69-102. Riewchotisakul, S. and Akkutlu, I. Y. 2016. Adsorption-Enhanced Transport of Hydrocarbons in Organic Nanopores. SPE J 21(6): 1960-1969. SPE-175107-PA. Rogers, L.C.G. and Williams, D. 1994. Diffusions, Markov processes and martingales: Volume 2. London, UK: Cambridge University Press. Sidje, R. B. 1998. A Software Package for Computing Matrix Exponentials. ACM Transactions on Mathematical Software (TOMS) 24(1): 130-156. Splichte. 2013, GitHub repository, https://github.com/splichte/lsi. Accessed December 18, 2017. Vilaseca, E., Isvoran, A., Madurga, S., et al. 2011. New Insights into Diffusion in 3D Crowded Media by Monte Carlo Simulations: Effect of Size, Mobility and Spatial Distribution of Obstacles. Physical Chemistry Chemical Physics 13(16):7396-7407. Vlahos, L., Isliker, H., Kominis, Y., et al. 2008. Normal and Anomalous Diffusion: A Tutorial. Warren, J. E. and Root, P. J. 1963. The Behavior of Naturally Fractured Reservoirs. SPE J 3(3): 245-255. SPE-426- PA. Wikipedia. 2017. Transition Rate Matrix. https://en.wikipedia.org. Accessed December 18, 2017. Wu, K. and Chen, Z. 2016. Real Gas Transport through Complex Nanopores of Shale Gas Reservoirs. Paper SPE- 180086-MS presented at the SPE Europe featured at 78th EAGE Conference and Exhibition, Vienna, Austria, 30 May-2 June. Yu, W., Xu, Y., Weijermars, R., et al. 2017. Impact of Well Interference on Shale Oil Production Performance: A Numerical Model for Analyzing Pressure Response of Fracture Hits with Complex Geometries. Paper presented at the SPE Hydraulic Fracturing Technology Conference and Exhibition, The Woodlands, Texas, 24–26 January. SPE-184825-MS. Zhang, Y., Meerschaert, M.M., and Baeumer, B. 2008. Particle tracking for time-fractional diffusion. Physical Review E 78(3): 36-42. Shuai Liu is a PhD-degree candidate in the department of petroleum engineering at Texas A&M University. His research is about the fracture design and optimization in unconventional reservoirs. Liu holds a master’s and bachelor’s degree in petroleum engineering from China University of Petroleum, Beijing. Han Li is a seismic imaging geophysicist working at CGG in Houston, Texas. He got his Ph.D degree in petroleum engineering department at Texas A&M University in 2017. He co-authored more than 12 articles and conference papers. His expertise is hydraulic fracturing, rock mechanics, numerical simulation and geophysics. https://github.com/splichte/lsi https://en.wikipedia.org/w/index.php?title=Transition_rate_matrix&oldid=815062240 24 Peter P. Valkó is the holder of the R. Whiting Chair in Petroleum Engineering at Texas A&M University. He holds B.S. and M.S. degrees from Veszprém University, Hungary, and Ph.D. degree from the Institute of Catalysis, Novosibirsk, Russia. Previously he taught at academic Institutions in Austria and Hungary and worked for the Hungarian Oil Company (MOL). His research interests include hydraulic fracturing, performance of stimulated wells and numerical inversion methods.