A Markov-Chain-Based Method to Characterize Anomalous Diffusion Phenomenon in Unconventional Reservoir

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 socalled "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 ultralow permeability and complex fracture system.In turn, academic progress can shed light on the flow mechanism 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/dualpermeability 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 singlephase 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 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.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. 5.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 halflength   .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.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, as Eq.1 shows.where  = ,  and   is the pressure in domain  .Since the properties are constant in each domain, Eq.1 can be rearranged into the form as is shown in Eq.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 m 2 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.
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 singlephase 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.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 welldeveloped 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.

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).
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.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 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 Both definitions in Eq.9 and Eq.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 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.
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. 1 is the th entry of this vector.
⃑ (0) = [0, 0, ⋯ , 0, 1, 0, ⋯ , 0] ……………………( 16) 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. .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 ( ~10 2 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 preprocessing 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°; 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.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 vs.  .The plots and parameters of the fracture networks, and the results are shown below.

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 straightforward 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 ft 2 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 ft 2 MSD is reached due to the totally homogeneous fracture distribution.In other cases, most deviations happen within the range 10 ft 2 ~ 10 3 ft 2 MSD, 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.
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.

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;

FigureFigure 3 -
Figure 2-Illustration of the displacements with respect to different coordinates

Figure 4 -
Figure 4-Map view plot of the fracture network and 10 sampled nodes in Case 0

Figure 6 -FigureFigure 8
Figure 6-Map view plot of the fracture network and 10 sampled nodes in Case 1

Figure 9 -Figure
Figure 9-Map view plot of the fracture network and 10 sampled nodes in Case 2

Figure 12 -Figure
Figure 12-Map view plot of the fracture network and 10 sampled nodes in Case 3

Table 1 -
Statistics for the fracture network in various Cases 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.