Patent Publication Number: US-11643924-B2

Title: Determining matrix permeability of subsurface formations

Description:
TECHNICAL FIELD 
     The present disclosure relates to characteristics of subsurface formations, particularly to matrix permeability of subsurface formations. 
     BACKGROUND 
     One of the most important parameters for characterizing unconventional or source-rock reservoirs is permeability. The permeability includes: fracture permeability and matrix permeability. Fracture permeability characterizes the ease with which fluid can flow through the natural fractures in the rock, as well as those fractures created through hydrofracturing. Matrix permeability characterizes the ease with which a fluid can flow within the intact portion of the rock. Fracture permeability usually falls within the millidarcy (mD) range, whereas matrix permeability is in the microdarcy (μD) to nanodarcy (nD) range. Although matrix permeability is much lower than fracture permeability, under some circumstances matrix permeability may control the production of a rock reservoir. 
     Unconventional rocks are highly laminated, and many weakness planes can exist approximately parallel to the lamination due to the orientation of mineral components and distribution of pores and organic matters. When rock samples are retrieved to the earth surface from a reservoir, the rock samples tend to have more fractures due to the removal of the overburden pressure. Fractures generated through this process are called induced fractures relative to in situ fractures existing in rocks naturally. Thus, measuring matrix permeability of reservoir rocks can be challenging due to the presence of artificially generated fractures and microfractures. 
     SUMMARY 
     The present specification describes methods, apparatus, and systems for determining matrix permeability of subsurface formations, e.g., fractured reservoir rocks. 
     One aspect of the present disclosure features a method of determining permeability of a subsurface formation, including: positioning a sample of the subsurface formation in a measurement cell, fluidly connecting an inlet and an outlet of the sample to an upstream reservoir and a downstream reservoir, respectively, flowing a fluid through the sample from the upstream reservoir to the downstream reservoir, measuring changes of an upstream pressure associated with the upstream reservoir and a downstream pressure associated with the downstream reservoir in a measurement time period, and determining a matrix permeability of the subsurface formation based on measurement data before the upstream pressure and the downstream pressure merge at a merging time point. 
     The subsurface formation can include at least one of shale, limestone, siltstone, or sandstone. The fluid can include one or more gases from a group of gases including methane, argon, nitrogen, carbon dioxide, helium, ethane, and propane. 
     In some implementations, determining the matrix permeability of the subsurface formation includes: determining an initial time point when the upstream pressure and the downstream pressure are larger than a central pressure at an axial center of the sample, e.g., the middle point of the axis of a cylindrical sample, extracting measurement data within a particular time range from the initial time point to an ending time point, and determining the matrix permeability of the subsurface formation using the extracted measurement data within the particular time range. The initial time point is earlier than the merging time point. The ending time point can be earlier than the merging time point. The ending time point can be later than the merging time point. 
     In some examples, determining the initial time point includes: determining the initial time point to be a time point no earlier than a declining time point when both the upstream pressure and the downstream pressure start to decline. 
     In some examples, determining the starting time point includes: determining the central pressure at the axial center of the sample using the upstream pressure and the downstream pressure according to a formula: 
                 P   ms     =         -     (           V   ur       V   p       ⁢     3   2       +     1     4             )       ⁢     P   ur       -       (           V     d   ⁢   r         V   p       ⁢     3   2       +     1   4       )     ⁢     P   dr       +       3   2     ⁢     (         V   p     +     V   ur     +     V   dr         V   p       )     ⁢     P   ∞           ,         
where P ms , P ur , and P dr  represent the central pressure, the upstream pressure, and the downstream pressure, respectively, V ur  and V dr  represent volumes of the upstream reservoir and the downstream reservoir, respectively, V p  represents a total volume of pores in the sample, and P ∞  represents an equilibrium pressure at an end of the measurement time period when P ur  and P dr  are stable after merging.
 
     In some implementations, determining the matrix permeability of the subsurface formation using the extracted measurement data within the particular time range includes: determining the matrix permeability based on data arrays G(t) and f(t) obtained from the extracted measurement data within the particular time range, where 
                 G   ⁡   (   t   )     =         (       5   ⁢   AA     +     1   2       )     ⁢     (         P   ur     (   t   )     -       P   ur     (     t   0     )       )       +       (       5   ⁢   BB     +     1   2       )     ⁢     (         P   dr     (   t   )     -       P   dr     (     t   0     )       )           ,         
and
 
                 f   ⁡   (   t   )     =       ∫     t   0     t         1     μ   ⁢   Φ   ⁢     C   g     ⁢     L   2         ⁢   4   ⁢   8   ⁢     (         (     AA   +     1   2       )     ⁢     P   ur       +       (     BB   +     1   2       )     ⁢     P   dr       -       (     1   +     A   ⁢   A     +     B   ⁢   B       )     ⁢     P   ∞         )     ⁢   d   ⁢   t         ,         
where AA=V ur /V p , BB=V dr /V p  V ur  and V dr  represent volumes of the upstream reservoir and the downstream reservoir, respectively, V p  represents a total volume of pores in the sample, t 0  represents the initial time point, t represents a time point within the particular time range, μ represents a viscosity of the fluid (e.g., a gas), Φ represents a porosity of the sample, L represents a length of the sample, and C g  represents a compressibility of the fluid (e.g., a gas), P ur  and P dr  represent the upstream pressure and the downstream pressure, respectively, and P ∞  represents an equilibrium pressure at an end of the measurement time period when P ur  and P dr  are stable after P ur  and P dr  merge.
 
     In some implementations, determining the matrix permeability based on the data arrays includes determining the matrix permeability according to a formula:
 
 G ( t )=− kf ( t )+ C,  
 
where k represents the matrix permeability, and C represents a constant.
 
     Determining the matrix permeability based on the data arrays can include: performing linear regression or linear fitting on the data arrays to get a slope, and determining the matrix permeability to be a negative of the slope. 
     In some implementations, the method can further include: before flowing the fluid through the sample, confining the measurement cell with a confining fluid to cause a confining pressure in the measurement cell where the confining fluid can be separated from the flowing fluid by one or more layers of sleeves, and applying an effective stress pressure to the sample, e.g., from all directions. Flowing the fluid through the sample can include: flowing the fluid to apply a pore pressure to the sample, where the pore pressure is smaller than the confining pressure. The method can further include: maintaining a pressure difference between the pore pressure and the confining pressure to be the effective stress pressure while flowing the fluid through the sample from the upstream reservoir to the downstream reservoir. 
     In some implementations, measuring the changes of the upstream pressure and the downstream pressure in the measurement time period includes: carrying out a pressure pulse decay permeability (PDP) measurement. The method can further include: preparing the sample of the subsurface formation by sealing a core plug sample of the subsurface formation with a sleeve. The method can also include: determining a fracture permeability of the subsurface formation based on the measurement data before the upstream pressure and the downstream pressure merge. 
     Another aspect of the present disclosure features a system for determining permeability of a subsurface formation. The system includes a measuring system and a computing system. The measuring system includes a measurement cell configured to contain a sample of the surface formation and upstream and downstream reservoirs fluidly connecting to an inlet and an outlet of the sample, respectively. The measuring system is configured to: flow a fluid through the sample from the upstream reservoir to the downstream reservoir, and measure changes of an upstream pressure associated with the upstream reservoir and a downstream pressure associated with the downstream reservoir in a measurement time period. The computing system is configured to determine a matrix permeability of the subsurface formation based on measurement data before the upstream pressure and the downstream pressure merge at a merging time point. 
     In some implementations, the computing system is separate from the measuring system. In some implementations, the measuring system includes the computing system. 
     In some implementations, the measuring system is configured to: before flowing the fluid through the sample, confine the measurement cell with a confining fluid to cause a confining pressure in the measurement cell where the confining fluid can be separated from the flowing fluid by one or more layers of sleeves, apply an effective stress pressure to the sample, e.g., from all directions, flow the fluid to apply a pore pressure to the sample, where the pore pressure is smaller than the confining pressure, and maintain a pressure difference between the confining pressure and the pore pressure to be the effective stress pressure while flowing the fluid through the sample from the upstream reservoir to the downstream reservoir. The measuring system can be configured to measure the changes of the upstream pressure and the downstream pressure in the measurement time period by carrying out a pressure pulse decay permeability (PDP) measurement. 
     In some implementations, the computing system is configured to: determine an initial time point when the upstream pressure and the downstream pressure are declining and are larger than a central pressure at an axial center of the sample (e.g., a middle point of the axis of the sample), extract measurement data within a particular time range from the initial time point to an ending time point, e.g., earlier or later than the merging time point, and determine the matrix permeability of the subsurface formation using the extracted measurement data within the particular time range. 
     In some examples, the computing system is configured to: determine data arrays G(t) and f(t) obtained from the extracted measurement data within the particular time range, where 
                 G   ⁡   (   t   )     =         (       5   ⁢   AA     +     1   2       )     ⁢     (         P   ur     (   t   )     -       P   ur     (     t   0     )       )       +       (       5   ⁢   BB     +     1   2       )     ⁢     (         P   dr     (   t   )     -       P   dr     (     t   0     )       )           ,         
and
 
                 f   ⁡   (   t   )     =       ∫     t   0     t         1     μ   ⁢   Φ   ⁢     C   g     ⁢     L   2         ⁢   4   ⁢   8   ⁢     (         (     AA   +     1   2       )     ⁢     P   ur       +       (     BB   +     1   2       )     ⁢     P   dr       -       (     1   +     A   ⁢   A     +     B   ⁢   B       )     ⁢     P   ∞         )     ⁢   dt         ;         
and determine the matrix permeability based on the data arrays according to a formula:
 
 G ( t )=− kf ( t )+ C,  
 
where k represents the matrix permeability, C represents a constant, AA=V ur /V p , BB=V dr /V p , V ur  and V dr  represent volumes of the upstream reservoir and the downstream reservoir, respectively, V p  represents a total volume of pores in the sample, t 0  represents the initial time point, t represents a time point within the particular time range, μ represents gas viscosity, Φ represents a porosity of the sample, L represents a length of the sample, and C g  represents gas compressibility, P ur  and P dr  represent the upstream pressure and the downstream pressure, respectively, and P ∞  represents an equilibrium pressure at an end of the measurement time period when P ur  and P dr  are stable after P ur  and P dr  merge.
 
     Implementations of the above techniques include methods, systems, computer program products and computer-readable media. In one example, a method can be performed by at least one processor coupled to at least one non-volatile memory and the methods can include the above-described actions. In another example, one such computer program product is suitably embodied in a non-transitory machine-readable medium that stores instructions executable by one or more processors. The instructions are configured to cause the one or more processors to perform the above-described actions. One such computer-readable medium stores instructions that, when executed by one or more processors, are configured to cause the one or more processors to perform the above-described actions. 
     The details of one or more implementations of the subject matter of this specification are set forth in the accompanying drawings and associated description. Other features, aspects, and advantages of the subject matter will become apparent from the description, the drawings, and the claims. 
    
    
     
       DESCRIPTION OF DRAWINGS 
         FIG.  1    is a schematic diagram illustrating an example system for determining matrix permeability of a subsurface formation, according to one or more example embodiments of the present disclosure. 
         FIG.  2 A  illustrates a cylindrical core sample used in the systems and methods, according to one or more example embodiments of the present disclosure. 
         FIG.  2 B  illustrates a data analysis model for the sample of  FIG.  2 A , according to one example embodiment of the disclosure. 
         FIG.  3 A  shows a diagram illustrating transient pressures in upstream and downstream sides of a sample during an example measurement, according to one or more example embodiments of the present disclosure. 
         FIG.  3 B  shows a plot between data arrays G(t) and f(t) based on data selected from the measurement data of  FIG.  3 A  before the upstream pressure and the downstream pressure merge, according to one or more example embodiments of the present disclosure. 
         FIG.  4    is a flowchart of an example process of determining a matrix permeability of a subsurface formation, according to one or more example embodiments of the present disclosure. 
     
    
    
     DETAILED DESCRIPTION 
     As the presence of fractures may enlarge the permeability of a source rock sample with induced fractures, matrix permeability measurement results from a steady-state permeability (SSP) method and a transient pulse-decay permeability (PDP) method can be much larger than a true matrix permeability of reservoir rocks. The permeability thus measured hereinafter can be referred to fracture permeability. 
     In some cases, a Gas Research Institute (GRI) method resorts to crushing rocks such that the crushed particles are small enough so that no fractures cut into or cut through the particles and permeability from pressure decay on the particles represents the matrix permeability. However, the GRI method cannot measure permeability under different stress conditions. In some cases, a PDP method is applied to a single continuum sample, and measures fracture permeability if the sample is fractured and measures matrix permeability if the sample is not fractured. 
     In some cases, a modified PDP method measures both the fracture permeability and matrix permeability under different stress conditions from fractured unconventional rocks using data obtained from a late time period after the upstream pressure and downstream pressure merge into the same pressure for a given time. However, the sensitivity of the evolution of merged pressure measurements to matrix permeability may be not high. The merged pressure can change very slowly with time, which may also create problems for accurately measuring the rock matrix permeability in practice. The flow into the matrix may occur even before pressures merge. 
     Implementations of the present disclosure provide methods, systems and techniques for determining matrix permeability from fractured source rock samples using data obtained before upstream pressure and downstream pressure merge, which can resolve the above issues and provide quicker, higher sensitivity, and higher accuracy for measurements of matrix permeability. 
     In some implementations, for data analysis, a source rock sample is prepared. The sample can be represented into upstream and downstream halves along an axis of the sample. A pore pressure at a central of the sample can be calculated based on mass balance using pressures in upstream and downstream reservoirs and a final merged pressure. Pressure data after both the upstream pressure and downstream pressure decline but before the upstream pressure and downstream pressure merge can be obtained. Volume ratios between the volumes of the upstream and downstream reservoirs and the pore volume of the sample can be also obtained. Then, based on Darcy law and the summation of gas mass fluxes, two arrays G(t) and f(t) that are related by matrix permeability can be obtained. By performing linear regression or plotting, matrix permeability of the sample can be determined. The obtained matrix permeability can be used in reservoir simulation models to predict well productions from corresponding unconventional reservoirs. 
     The techniques can address existing challenges for measuring permeability of unconventional formations. For examples, the techniques can address a challenge of coring and core handling of heterogeneous rock samples that create extensive microcracking. The presence of these microcracks may directly affect the permeability measured, the lower the rock permeability, the greater the effect of the induced microcracks. This effect is most prevalent for laminated, low permeability, organic-rich, mudstones, where the organic to mineral contact and the interfaces associated with the laminated fabric are weak contacts that are prone to part during unloading. The methods implemented in the present disclosure are performed on source rock samples with induced fractures. The effect of the induced microcracks may be reduced or eliminated. 
     The techniques can also address the heterogeneity challenge in measuring permeability of unconventional formations, low permeability rocks. These rocks possess intrinsic variability in texture and composition that results from geologic processes of deposition and diagenesis. As a result, these rocks exhibit a broad distribution of permeabilities. Conventional permeability measurements developed for homogeneous media, have focused on the evaluation of a single representative value of permeability, without accounting for the distribution of permeabilities. The resulting consequences are that the “single permeability” is ill-defined and not necessarily representative of the rock containing the distribution of permeabilities. In comparison, the methods implemented in the present disclosure can be performed on multiple continua samples with induced fractures to get both fracture permeability and matrix permeability in the same measurement. 
     The techniques can also address the possibility of underestimating the matrix permeability, on samples of standard size, for example, 1 to 1.5 inches in diameter and 1 to 2 inches in length. The method using crushed fragments of sample tends to be the standard method most often used for measuring permeability in ultra-low permeability rocks. However, the measured permeabilities of the crushed sample fragments do not represent the mean value of the whole permeability distribution of the rock before it was crushed, unless a further calibration or correction is made to these measurements. If tiny fractures do exist in the reservoir conditions, crushing the rock to remove the effect of fractures and performing the permeability pressure decay measurements may generate much smaller permeability values then should be. In comparison, the methods implemented in the present disclosure can be performed with a whole plug sample with or likely with fractures or microfractures, instead of crushed samples having different particle size fractions. Thus, the data can be obtained at a high pore pressure (e.g., larger than 100 pounds per square inch—psi) and requires no or little corrections. 
     Example Measurement System 
       FIG.  1    illustrates an example system  100  for determining matrix permeability of a subsurface formation, according to one or more example embodiments of the present disclosure. The system  100  may include an upstream reservoir  110 , a downstream reservoir  120 , and a measurement cell  130  (or a core holder) connecting to the upstream and downstream reservoirs  110 ,  120  with valves  122 ,  124  through a fluidic channel  106  (e.g., a tube). The fluid in  110 , 120 , and  106  can be marked as a pore fluid. 
     The measurement cell  130  can be a vessel containing a confining fluid  140  (e.g., a gas or a liquid). The measurement cell  130  can be connected to a pump to control a confining pressure caused by the fluid  140  in the measurement cell  130 . A sample  150  of subsurface formation can be placed in the measurement cell  130  and can be subject to the confining pressure. In some embodiments, the subsurface formation includes at least one of shale, limestone, siltstone, or sandstone. The confining fluid  140  and the pore fluid can be separated by one or more layers of sleeves surrounding the sample  150  and metal piece at the ends of the sample  150 . 
     The upstream reservoir  110  and the downstream reservoir  120  can be fluid reservoirs and be connected to two different pumps to control pressures, respectively. The fluid (i.e., the pore fluid) can include one or more gases including at least one of methane, argon, nitrogen, carbon dioxide, helium, ethane, or propane. The pumps connected to the upstream reservoir  110  and the downstream reservoir  120  can include flow meters to measure mass flow rates at an inlet and outlet of the sample  150 . A number of pressure and temperature transducers can be used to monitor pressure and temperature conditions at different locations in the system  100 . For example, a pressure gauge  160  can be connected to the upstream reservoir  110 , a pressure gauge  162  can be connected between the valve  122  and an inlet of the measurement cell  130 , and a pressure gauge  164  can be connected between the value  124  and an outlet of the measurement cell  130 . 
     The sample  150  may be fluidly connected to the supply of fluid from the upstream and downstream reservoirs  110 ,  120  via opposite ends of the sample  150 , as shown in  FIG.  1   . In some examples, the sample  150  is a core plug sample of unconventional shale rock, e.g., a horizontal sample plug with a bedding parallel to an axis of a cylindrical shape of the sample  150 . The sample  150  can be with one or more layers of sleeves for sealing. For example, the sample  150  can be wrapped in layers of plastic and metal sheet and put between two end pieces with metal lines to connect the sample  150  and the upstream and downstream reservoirs  110 ,  120  while keeping the sample  150  well-sealed. 
     Example Data Analysis Model 
       FIG.  2 A  illustrates a cylindrical core plug sample  150  used in the systems and methods, according to one or more example embodiments of the present disclosure. The sample  150  has a length L and a central axis  151 .  FIG.  2 B  illustrates a data analysis model  200  for the sample  150  of  FIG.  2 A , according to one or more example embodiments of the present disclosure. 
     The determination of rock matrix permeability can rely on analysis of measurement data with a relationship between the data and the permeability. As fluid flow in a rock sample is a one-dimensional process, the analysis of fluid mass flow for upstream side and downstream side of the sample can be done separately in one dimensional coordinate system shown in  FIG.  2 B . For illustration purposes only, gas is used as an example of the flowing fluid, i.e., the pore fluid. The data analysis is done on each half of the sample  150  (upstream half  152  and downstream half  154 ). The one dimensional system for both halves  152 ,  154  are shown with the original point (or center point) (i.e., x=0) at an axial center of the sample  150 . x is a distance from the original point (x=0) to a position within the sample  150 . For example, L/2 is the distance from the original point (or axial center  155 ) to end faces  153 ,  157  of the sample  150 . 
     P us  and P ds  respectively represent pressures in the upstream half  152  and the downstream half  154  of the sample  150 . The pressures at position x from the axial center (x=0) and time t can be expressed as P us  (x,t) and P ds (x,t) for the upstream half  152  and the downstream half  154 , respectively. For simplicity, P us  (x,t) and P ds (x,t) are expressed as P us  and P ds , respectively. P ms  represents the pressure at the axial center  155  of the sample  150 , i.e., at the original point x=0. The gas density ρ(x,t) is expressed as ρ. l is the half length of the sample  150  with the full length of L. x is in the range of (0, L/2) and z is defined as x/l. V p  stands for a total volume of pores in the sample  150 . A represents a cross-sectional area of the sample  150 . Φ represents a porosity of the sample  150 , which can be expressed as the fraction of the pore volume Vp relative to the total rock volume V bulk  of the sample  150 . That is, V bulk =LA; and 
             ϕ   =         V   p       V   bulk       .           
V ur  and V dr  represent volumes of the upstream and downstream reservoirs  110 ,  120 , respectively. P ur (t) and P dr (t) respectively represent pressures in the upstream reservoir  110  and the downstream reservoir  120 , which are changing with time t and are expressed as P ur  and P dr , respectively.
 
     At a particular time period, both the pressures P ur  and P dr  in the upstream and downstream reservoirs  110 ,  120  are larger than the pressure P ms  in the center of the sample  150 . Gas flows into the sample  150  from both the upstream reservoir  110  and the downstream reservoir  120 . The pore volume V p  of the sample  150  can be expressed as:
 
 V   p   =ΦAL= 2Φ Al   (1).
 
     Considering that there can be an approximate quadratic distribution of pressure from the end surface  153 ,  157  of the sample to the center  155  of the sample  150 , P us  and P ds  in the upstream half  152  and the downstream half  154  can be expressed as: 
     
       
         
           
             
               
                 
                   
                     
                       P 
                       
                         u 
                         ⁢ 
                         s 
                       
                     
                     = 
                     
                       
                         
                           
                             ( 
                             
                               
                                 P 
                                 ur 
                               
                               - 
                               
                                 P 
                                 
                                   m 
                                   ⁢ 
                                   s 
                                 
                               
                             
                             ) 
                           
                           ⁢ 
                           
                             
                               ( 
                               
                                 x 
                                 
                                   L 
                                   2 
                                 
                               
                               ) 
                             
                             2 
                           
                         
                         + 
                         
                           P 
                           
                             m 
                             ⁢ 
                             s 
                           
                         
                       
                       = 
                       
                         
                           
                             ( 
                             
                               
                                 P 
                                 ur 
                               
                               - 
                               
                                 P 
                                 
                                   m 
                                   ⁢ 
                                   s 
                                 
                               
                             
                             ) 
                           
                           ⁢ 
                           
                             z 
                             2 
                           
                         
                         + 
                         
                           P 
                           
                             m 
                             ⁢ 
                             s 
                           
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   2 
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     P 
                     
                       d 
                       ⁢ 
                       s 
                     
                   
                   = 
                   
                     
                       
                         
                           ( 
                           
                             
                               P 
                               dr 
                             
                             - 
                             
                               P 
                               
                                 m 
                                 ⁢ 
                                 s 
                               
                             
                           
                           ) 
                         
                         ⁢ 
                         
                           
                             ( 
                             
                               x 
                               
                                 L 
                                 2 
                               
                             
                             ) 
                           
                           2 
                         
                       
                       + 
                       
                         P 
                         
                           m 
                           ⁢ 
                           s 
                         
                       
                     
                     = 
                     
                       
                         
                           ( 
                           
                             
                               P 
                               dr 
                             
                             - 
                             
                               P 
                               
                                 m 
                                 ⁢ 
                                 s 
                               
                             
                           
                           ) 
                         
                         ⁢ 
                         
                           z 
                           2 
                         
                       
                       + 
                       
                         
                           P 
                           
                             m 
                             ⁢ 
                             s 
                           
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     A mass balance equation for either half  152 ,  154  of the sample can be written as: 
                         ∂     (     ρ   ⁢   Φ     )         ∂   t       =       -       ∂   q       ∂   x         =       -       ∂   q       ∂   z         ⁢     1   l           ,           (   4   )               
where q is the gas mass flux. If the change of porosity c during the measurement can be neglected, there can be:
 
                         ∂     (       ρ     u   ⁢   s       ⁢   Φ     )         ∂   t       =       Φ   ⁢       ∂     ρ     u   ⁢   s           ∂   t         =       Φ   ⁢       ∂     ρ     u   ⁢   s           ∂     P     u   ⁢   s           ⁢       ∂     P     u   ⁢   s           ∂   t         =       Φρ   ⁢     C   g     ⁢       ∂     P     u   ⁢   s           ∂   t         =     Φ   ⁢   ρ   ⁢     C   g     ⁢       d   ⁡   (         (       P     u   ⁢   r       -     P     m   ⁢   s         )     ⁢     z   2       +     P     m   ⁢   s         )       d   ⁢   t                 ,           (   5   )                                 ∂     (       ρ     d   ⁢   s       ⁢   Φ     )         ∂   t       =       Φ   ⁢       ∂     ρ     d   ⁢   s           ∂   t         =       Φ   ⁢       ∂     ρ     d   ⁢   s           ∂     P     d   ⁢   s           ⁢       ∂     P     d   ⁢   s           ∂   t         =       Φ   ⁢   ρ   ⁢     C   g     ⁢       ∂     P     d   ⁢   s           ∂   t         =     Φ   ⁢   ρ   ⁢     C   g     ⁢       d   ⁡   (         (       P     d   ⁢   r       -     P     m   ⁢   s         )     ⁢     z   2       +     P     m   ⁢   s         )       d   ⁢   t                 ,           (   6   )               
where
 
                 ρ   ⁢     C   g       =         ∂     ρ     u   ⁢   s           ∂     P     u   ⁢   s           =       ∂     ρ     d   ⁢   s           ∂     P     d   ⁢   s               ,         
ρ represents a density of the pore fluid such as a gas, and C g  represents a compressibility of the pore fluid such as a gas. By definition,
 
                 C   ⁢   g     =       ∂   ρ       ρ   ⁢     ∂   P           ,         
thus
 
     
       
         
           
             
               
                 ρ 
                 ⁢ 
                 
                   C 
                   g 
                 
               
               = 
               
                 
                   ∂ 
                   ρ 
                 
                 
                   ∂ 
                   P 
                 
               
             
             . 
           
         
       
     
     The total mass change Q1 in the sample  150  (including both halves  152 ,  154 ) can be expressed as: 
     
       
         
           
             
               
                 
                   
                     Q 
                     ⁢ 
                     1 
                   
                   = 
                   
                     
                       
                         - 
                         Φρ 
                       
                       ⁢ 
                       
                         C 
                         g 
                       
                       ⁢ 
                       l 
                       ⁢ 
                       A 
                       ⁢ 
                       
                         
                           ∫ 
                           0 
                           1 
                         
                         
                           
                             
                               d 
                               ⁡ 
                               ( 
                               
                                 
                                   
                                     ( 
                                     
                                       
                                         P 
                                         dr 
                                       
                                       - 
                                       
                                         P 
                                         
                                           m 
                                           ⁢ 
                                           s 
                                         
                                       
                                     
                                     ) 
                                   
                                   ⁢ 
                                   
                                     z 
                                     2 
                                   
                                 
                                 + 
                                 
                                   P 
                                   
                                     m 
                                     ⁢ 
                                     s 
                                   
                                 
                                 + 
                                 
                                   
                                     ( 
                                     
                                       
                                         P 
                                         ur 
                                       
                                       - 
                                       
                                         P 
                                         
                                           m 
                                           ⁢ 
                                           s 
                                         
                                       
                                     
                                     ) 
                                   
                                   ⁢ 
                                   
                                     z 
                                     2 
                                   
                                 
                                 + 
                                 
                                   P 
                                   
                                     m 
                                     ⁢ 
                                     s 
                                   
                                 
                               
                               ) 
                             
                             dt 
                           
                           ⁢ 
                           dz 
                         
                       
                     
                     = 
                     
                       
                         - 
                         
                           
                             V 
                             p 
                           
                           2 
                         
                       
                       ⁢ 
                       ρ 
                       ⁢ 
                       
                         C 
                         g 
                       
                       ⁢ 
                       
                         
                           
                             d 
                             ⁡ 
                             ( 
                             
                               
                                 
                                   4 
                                   3 
                                 
                                 ⁢ 
                                 
                                   P 
                                   ms 
                                 
                               
                               + 
                               
                                 
                                   1 
                                   3 
                                 
                                 ⁢ 
                                 
                                   P 
                                   ur 
                                 
                               
                               + 
                               
                                 
                                   1 
                                   3 
                                 
                                 ⁢ 
                                 
                                   P 
                                   dr 
                                 
                               
                             
                             ) 
                           
                           dt 
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   7 
                   ) 
                 
               
             
           
         
       
     
     The total mass change Q2 in both upstream and downstream reservoirs  110 ,  120  can be expressed as: 
     
       
         
           
             
               
                 
                   
                     Q 
                     ⁢ 
                     2 
                   
                   = 
                   
                     
                       
                         - 
                         
                           V 
                           ur 
                         
                       
                       ⁢ 
                       ρ 
                       ⁢ 
                       
                         C 
                         g 
                       
                       ⁢ 
                       
                         
                           d 
                           ⁢ 
                           
                             P 
                             
                               u 
                               ⁢ 
                               r 
                             
                           
                         
                         
                           d 
                           ⁢ 
                           t 
                         
                       
                     
                     - 
                     
                       
                         V 
                         dr 
                       
                       ⁢ 
                       ρ 
                       ⁢ 
                       
                         C 
                         g 
                       
                       ⁢ 
                       
                         
                           
                             d 
                             ⁢ 
                             
                               P 
                               
                                 d 
                                 ⁢ 
                                 r 
                               
                             
                           
                           
                             d 
                             ⁢ 
                             t 
                           
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   8 
                   ) 
                 
               
             
           
         
       
     
     As a result of mass balance, Q2=−Q1, thus: 
     
       
         
           
             
               
                 
                   
                     
                       
                         - 
                         
                           V 
                           ur 
                         
                       
                       ⁢ 
                       ρ 
                       ⁢ 
                       
                         C 
                         g 
                       
                       ⁢ 
                       
                         
                           d 
                           ⁢ 
                           
                             P 
                             
                               u 
                               ⁢ 
                               r 
                             
                           
                         
                         
                           d 
                           ⁢ 
                           t 
                         
                       
                     
                     - 
                     
                       
                         V 
                         d 
                       
                       ⁢ 
                       ρ 
                       ⁢ 
                       
                         C 
                         g 
                       
                       ⁢ 
                       
                         
                           d 
                           ⁢ 
                           
                             P 
                             
                               d 
                               ⁢ 
                               r 
                             
                           
                         
                         
                           d 
                           ⁢ 
                           t 
                         
                       
                     
                   
                   = 
                   
                     
                       
                         V 
                         ⁢ 
                         p 
                       
                       2 
                     
                     ⁢ 
                     ρ 
                     ⁢ 
                     
                       C 
                       g 
                     
                     ⁢ 
                     
                       
                         
                           d 
                           ⁡ 
                           ( 
                           
                             
                               
                                 4 
                                 3 
                               
                               ⁢ 
                               
                                 P 
                                 
                                   m 
                                   ⁢ 
                                   s 
                                 
                               
                             
                             + 
                             
                               
                                 1 
                                 3 
                               
                               ⁢ 
                               
                                 P 
                                 
                                   u 
                                   ⁢ 
                                   r 
                                 
                               
                             
                             + 
                             
                               
                                 1 
                                 3 
                               
                               ⁢ 
                               
                                 P 
                                 
                                   d 
                                   ⁢ 
                                   r 
                                 
                               
                             
                           
                           ) 
                         
                         
                           d 
                           ⁢ 
                           t 
                         
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
     Equation 9 can be integrated as: 
     
       
         
           
             
               
                 
                   
                     
                       
                         V 
                         ur 
                       
                       ⁢ 
                       
                         P 
                         ur 
                       
                     
                     - 
                     
                       
                         V 
                         dr 
                       
                       ⁢ 
                       
                         P 
                         dr 
                       
                     
                   
                   = 
                   
                     
                       
                         V 
                         p 
                       
                       ( 
                       
                         
                           
                             4 
                             6 
                           
                           ⁢ 
                           
                             P 
                             ms 
                           
                         
                         + 
                         
                           
                             1 
                             6 
                           
                           ⁢ 
                           
                             P 
                             ur 
                           
                         
                         + 
                         
                           
                             1 
                             6 
                           
                           ⁢ 
                           
                             P 
                             dr 
                           
                         
                       
                       ) 
                     
                     + 
                     
                       C 
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   10 
                   ) 
                 
               
             
           
         
       
     
     When t=∞, P ur =P dr =P ms =P ∞ , where P ∞  is the pressure when t=∞ at the end of the measurement. Thus, C can be solved from Equation 10:
 
 C =(− V   ur   −V   dr   −V   p ) P   ∞   (11).
 
     Combining Equations 10 and 11 gives: 
     
       
         
           
             
               
                 
                   
                     P 
                     
                       m 
                       ⁢ 
                       s 
                     
                   
                   = 
                   
                     
                       
                         - 
                         
                           ( 
                           
                             
                               
                                 
                                   V 
                                   ur 
                                 
                                 
                                   V 
                                   p 
                                 
                               
                               ⁢ 
                               
                                 3 
                                 2 
                               
                             
                             + 
                             
                               1 
                               4 
                             
                           
                           ) 
                         
                       
                       ⁢ 
                       
                         P 
                         
                           u 
                           ⁢ 
                           r 
                         
                       
                     
                     - 
                     
                       
                         ( 
                         
                           
                             
                               
                                 V 
                                 dr 
                               
                               
                                 V 
                                 p 
                               
                             
                             ⁢ 
                             
                               3 
                               2 
                             
                           
                           + 
                           
                             1 
                             4 
                           
                         
                         ) 
                       
                       ⁢ 
                       
                         P 
                         dr 
                       
                     
                     + 
                     
                       
                         3 
                         2 
                       
                       ⁢ 
                       
                         ( 
                         
                           
                             
                               V 
                               p 
                             
                             + 
                             
                               V 
                               ur 
                             
                             + 
                             
                               V 
                               dr 
                             
                           
                           
                             V 
                             p 
                           
                         
                         ) 
                       
                       ⁢ 
                       
                         
                           P 
                           ∞ 
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
           
         
       
     
     From Equation 12, the following equation can be obtained: 
     
       
         
           
             
               
                 
                   
                     
                       d 
                       ⁢ 
                       
                         P 
                         
                           m 
                           ⁢ 
                           s 
                         
                       
                     
                     
                       d 
                       ⁢ 
                       t 
                     
                   
                   = 
                   
                     
                       
                         - 
                         
                           ( 
                           
                             
                               
                                 
                                   V 
                                   
                                     u 
                                     ⁢ 
                                     r 
                                   
                                 
                                 
                                   V 
                                   p 
                                 
                               
                               ⁢ 
                               
                                 3 
                                 2 
                               
                             
                             + 
                             
                               1 
                               4 
                             
                           
                           ) 
                         
                       
                       ⁢ 
                       
                         
                           d 
                           ⁢ 
                           
                             P 
                             
                               u 
                               ⁢ 
                               r 
                             
                           
                         
                         
                           d 
                           ⁢ 
                           t 
                         
                       
                     
                     - 
                     
                       
                         ( 
                         
                           
                             
                               
                                 V 
                                 
                                   d 
                                   ⁢ 
                                   r 
                                 
                               
                               
                                 V 
                                 p 
                               
                             
                             ⁢ 
                             
                               3 
                               2 
                             
                           
                           + 
                           
                             1 
                             4 
                           
                         
                         ) 
                       
                       ⁢ 
                       
                         
                           
                             d 
                             ⁢ 
                             
                               P 
                               
                                 d 
                                 ⁢ 
                                 r 
                               
                             
                           
                           
                             d 
                             ⁢ 
                             t 
                           
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
       
     
     Based on the mass balance, the summation of gas mass fluxes at positions x is given in the following equation: 
                     q   ⁢   1     =       -       ∫   0   x           ∂     (       ρ     u   ⁢   s       ⁢   Φ     )         ∂   t       ⁢   d   ⁢   x         -       ∫   0   x           ∂     (       ρ   ds     ⁢   Φ     )         ∂   t       ⁢     dx   .                   (   14   )               
Note that there are an x position in the upstream half  152  of the sample  150  and also an x position in the downstream half  154  of the sample  150 , as illustrated in  FIG.  2 B .
 
     Inserting Equations 5 and 6 into Equation 14 results in: 
     
       
         
           
             
               
                 
                   
                     
                       q 
                       ⁢ 
                       1 
                     
                     = 
                     
                       
                         - 
                         Φ 
                       
                       ⁢ 
                       ρ 
                       ⁢ 
                       
                         C 
                         g 
                       
                       ⁢ 
                       
                         l 
                         ⁡ 
                         ( 
                         
                           
                             2 
                             ⁢ 
                             z 
                             ⁢ 
                             
                               
                                 d 
                                 ⁢ 
                                 
                                   P 
                                   
                                     m 
                                     ⁢ 
                                     s 
                                   
                                 
                               
                               
                                 d 
                                 ⁢ 
                                 t 
                               
                             
                           
                           + 
                           
                             
                               1 
                               3 
                             
                             ⁢ 
                             
                               
                                 z 
                                 3 
                               
                               ( 
                               
                                 
                                   
                                     d 
                                     ⁢ 
                                     
                                       P 
                                       
                                         u 
                                         ⁢ 
                                         r 
                                       
                                     
                                   
                                   
                                     d 
                                     ⁢ 
                                     t 
                                   
                                 
                                 - 
                                 
                                   
                                     d 
                                     ⁢ 
                                     
                                       P 
                                       
                                         m 
                                         ⁢ 
                                         r 
                                       
                                     
                                   
                                   
                                     d 
                                     ⁢ 
                                     t 
                                   
                                 
                               
                               ) 
                             
                           
                           + 
                           
                             
                               1 
                               3 
                             
                             ⁢ 
                             
                               
                                 z 
                                 3 
                               
                               ( 
                               
                                 
                                   
                                     d 
                                     ⁢ 
                                     
                                       P 
                                       
                                         d 
                                         ⁢ 
                                         r 
                                       
                                     
                                   
                                   
                                     d 
                                     ⁢ 
                                     t 
                                   
                                 
                                 - 
                                 
                                   
                                     d 
                                     ⁢ 
                                     
                                       P 
                                       
                                         m 
                                         ⁢ 
                                         r 
                                       
                                     
                                   
                                   
                                     d 
                                     ⁢ 
                                     t 
                                   
                                 
                               
                               ) 
                             
                           
                         
                         ) 
                       
                     
                   
                   . 
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
           
         
       
     
     Based on Darcy law, the summation of the gas flow fluxes can also be expressed as: 
                       q   ⁢   2     =         -   k     ⁢     ρ   μ     ⁢       ∂     P     u   ⁢   s           ∂   z       ⁢     1   l       -     k   ⁢     ρ   μ     ⁢       ∂     P     d   ⁢   s           ∂   z       ⁢     1   l           ,           (   16   )               
where k represents the matrix permeability and μ represents gas viscosity.
 
     Based on q2=q1 and Equations 15 and 16, the following equation can be obtained: 
     
       
         
           
             
               
                 
                   
                     
                       k 
                       ⁢ 
                       
                         1 
                         μ 
                       
                       ⁢ 
                       
                         
                           ∂ 
                           
                             P 
                             
                               u 
                               ⁢ 
                               s 
                             
                           
                         
                         
                           ∂ 
                           z 
                         
                       
                       ⁢ 
                       
                         1 
                         l 
                       
                       ⁢ 
                       
                         1 
                         
                           Φ 
                           ⁢ 
                           
                             C 
                             g 
                           
                           ⁢ 
                           l 
                         
                       
                     
                     + 
                     
                       k 
                       ⁢ 
                       
                         1 
                         μ 
                       
                       ⁢ 
                       
                         
                           ∂ 
                           
                             P 
                             
                               d 
                               ⁢ 
                               s 
                             
                           
                         
                         
                           ∂ 
                           z 
                         
                       
                       ⁢ 
                       
                         1 
                         l 
                       
                       ⁢ 
                       
                         1 
                         
                           Φ 
                           ⁢ 
                           
                             C 
                             g 
                           
                           ⁢ 
                           l 
                         
                       
                     
                   
                   = 
                   
                     
                       2 
                       ⁢ 
                       z 
                       ⁢ 
                       
                         
                           d 
                           ⁢ 
                           
                             p 
                             
                               m 
                               ⁢ 
                               s 
                             
                           
                         
                         
                           d 
                           ⁢ 
                           t 
                         
                       
                     
                     + 
                     
                       
                         1 
                         3 
                       
                       ⁢ 
                       
                         
                           z 
                           3 
                         
                         ( 
                         
                           
                             
                               d 
                               ⁢ 
                               
                                 p 
                                 
                                   u 
                                   ⁢ 
                                   r 
                                 
                               
                             
                             
                               d 
                               ⁢ 
                               t 
                             
                           
                           - 
                           
                             
                               d 
                               ⁢ 
                               
                                 p 
                                 
                                   m 
                                   ⁢ 
                                   s 
                                 
                               
                             
                             
                               d 
                               ⁢ 
                               t 
                             
                           
                         
                         ) 
                       
                     
                     + 
                     
                       
                         1 
                         3 
                       
                       ⁢ 
                       
                         
                           
                             z 
                             3 
                           
                           ( 
                           
                             
                               
                                 d 
                                 ⁢ 
                                 
                                   p 
                                   
                                     d 
                                     ⁢ 
                                     r 
                                   
                                 
                               
                               
                                 d 
                                 ⁢ 
                                 t 
                               
                             
                             - 
                             
                               
                                 d 
                                 ⁢ 
                                 
                                   p 
                                   
                                     m 
                                     ⁢ 
                                     s 
                                   
                                 
                               
                               
                                 d 
                                 ⁢ 
                                 t 
                               
                             
                           
                           ) 
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   17 
                   ) 
                 
               
             
           
         
       
     
     Integrating the above equation 17 on z from 0 to 1 results in 
                         k       µΦC   g     ⁢     L   2         ⁢   4   ⁢   8   ⁢     (         (       A   ⁢   A     +     1   2       )     ⁢     P   ur       +       (       B   ⁢   B     +     1   2       )     ⁢     P   dr       -       (     1   +     A   ⁢   A     +     B   ⁢   B       )     ⁢     P   ∞         )       =         -     (       5   ⁢   AA     +     1   2       )       ⁢       d   ⁢     P     u   ⁢   r           d   ⁢   t         -       (       5   ⁢   BB     +     1   2       )     ⁢       d   ⁢     P     d   ⁢   r           d   ⁢   t             ,           (   18   )               
where AA=V ur /V p , and BB=V dr /V p .
 
     Integrating both sides of the above equation 18 from t 0  to t gives
 
 G ( t )=− kf ( t )+contestant  (19),
 
     
       
         
           
             
               
                 
                   
                     
                       G 
                       ⁡ 
                       ( 
                       t 
                       ) 
                     
                     = 
                     
                       
                         
                           ( 
                           
                             
                               5 
                               ⁢ 
                               AA 
                             
                             + 
                             
                               1 
                               2 
                             
                           
                           ) 
                         
                         ⁢ 
                         
                           ( 
                           
                             
                               
                                 P 
                                 ur 
                               
                               ( 
                               t 
                               ) 
                             
                             - 
                             
                               
                                 P 
                                 
                                   u 
                                   ⁢ 
                                   τ 
                                 
                               
                               ( 
                               
                                 t 
                                 0 
                               
                               ) 
                             
                           
                           ) 
                         
                       
                       + 
                       
                         
                           ( 
                           
                             
                               5 
                               ⁢ 
                               BB 
                             
                             + 
                             
                               1 
                               2 
                             
                           
                           ) 
                         
                         ⁢ 
                         
                           ( 
                           
                             
                               
                                 P 
                                 dr 
                               
                               ( 
                               t 
                               ) 
                             
                             - 
                             
                               
                                 P 
                                 dr 
                               
                               ( 
                               
                                 t 
                                 0 
                               
                               ) 
                             
                           
                           ) 
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   20 
                   ) 
                 
               
             
           
         
       
     
     
       
         
           
             
               
                 
                   
                     f 
                     ⁡ 
                     ( 
                     t 
                     ) 
                   
                   = 
                   
                     
                       ∫ 
                       
                         t 
                         0 
                       
                       t 
                     
                     
                       
                         1 
                         
                           μ 
                           ⁢ 
                           Φ 
                           ⁢ 
                           
                             C 
                             g 
                           
                           ⁢ 
                           
                             L 
                             2 
                           
                         
                       
                       ⁢ 
                       4 
                       ⁢ 
                       8 
                       ⁢ 
                       
                         ( 
                         
                           
                             
                               ( 
                               
                                 AA 
                                 + 
                                 
                                   1 
                                   2 
                                 
                               
                               ) 
                             
                             ⁢ 
                             
                               P 
                               ur 
                             
                           
                           + 
                           
                             
                               ( 
                               
                                 BB 
                                 + 
                                 
                                   1 
                                   2 
                                 
                               
                               ) 
                             
                             ⁢ 
                             
                               P 
                               dr 
                             
                           
                           - 
                           
                             
                               ( 
                               
                                 1 
                                 + 
                                 
                                   A 
                                   ⁢ 
                                   A 
                                 
                                 + 
                                 
                                   B 
                                   ⁢ 
                                   B 
                                 
                               
                               ) 
                             
                             ⁢ 
                             
                               P 
                               ∞ 
                             
                           
                         
                         ) 
                       
                       ⁢ 
                       
                         dt 
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   21 
                   ) 
                 
               
             
           
         
       
     
     In Equations 20 and 21, t 0  is a time before pressures P ur  and P dr  in both the upstream and downstream reservoirs  110  and  120  start to merge. Equations 19, 20 and 21 can be used for determining the matrix permeability k with measurement observations. An algorithm can be developed to evaluate G(t) and f(t) and then to determine the matrix permeability k. 
     Example Measurement of Matrix Permeability 
     Measurement of matrix permeability of a source rock (e.g., carbonate rock) sample can be performed using a measurement system, e.g., the measurement system  100  of  FIG.  1   . The source rock sample can have fractures. The measurement system can include a processor configured to process measurement data to determine the matrix permeability. 
     First, a core sample, e.g., the sample  150  of  FIGS.  1 ,  2 A and  2 B , is collected and prepared. The core sample can have a cylindrical shape. A dimension of the core sample can be measured and recorded. The dimension can include a length L of the core sample and a cross-sectional area A of the sample. The core sample can be sealed, e.g., with one or more layers of sleeves. For example, the core sample can be wrapped in layers of plastic and metal sheet and put between two end pieces with metal lines to connect the sample and upstream and downstream reservoirs, e.g., the reservoirs  110 ,  120  of  FIGS.  1  and  2 B . 
     Second, measurement on the core sample is performed. The core sample is loaded into a core holder or a measurement cell, e.g., the measurement cell  130  of  FIG.  1   . A confining pump can be used to introduce distilled water as a confining fluid, e.g., the fluid  140  of  FIG.  1   , into the measurement cell. Once the measurement cell is full, the confining pump continues to flow the distilled water into the measurement cell to ensure that there is no gas remaining inside the measurement cell. After the confining of the measurement cell is locked in, an effective stress pressure (e.g., a pressure difference between a confining pressure and a pore pressure) is immediately applied to the core sample to reduce the possibility of water leakage into the sleeve of the core sample. The effective stress pressure can be about 500 psi. After this step, the pore pressure and the confining pressure can be adjusted to the desired values. 
     Both upstream and downstream pumps that are connected to the core sample through upstream and downstream reservoirs, e.g., the reservoirs  110  and  120  of  FIG.  1   , are used to slowly impose a pore pressure to the core sample. The gas is allowed to flow through the sample from the upstream gas reservoir to the downstream gas reservoir. The pore pressure can be about 2,500 psi. The use of the pore pressure is to minimize the impact of diffusion on permeability measurements. As the gas from the upstream and downstream reservoirs is introduced into the core sample, the confining pressure in the measurement cell is increased simultaneously to maintain a pressure difference between the confining pressure and the pore pressure. The pressure difference can be the effective stress pressure, e.g., 500 psi. Once the confining pump and the upstream and downstream pumps are stable at a desired pressure, the core sample is isolated by closing the valves  102  and  104 . A pressure pulse can be applied in the upstream reservoir  110 . By opening the valve  102 , a pressure pulse decay permeability (PDP) measurement can be carried out. 
     During an early stage of the pressure pulse decay measurement on the core sample with fractures, the gas mainly flows within the fractures until the pressures in the upstream and downstream gas reservoirs become practically the same value. Once upstream and downstream reservoir pressures become the same, the pressure within the system including upstream reservoir, downstream reservoir and fracture reach equilibrium. After this time (referred to as a pressure merging point), flow in fractures essentially stops because of zero pressure gradient, but flow still occurs from gas reservoirs into the rock matrix, resulting in further pressure decline. 
       FIG.  3 A  is a diagram  300  showing changes of pressures P ur  and P dr  of upstream and downstream reservoirs with time during an example measurement. Plot  310  shows the change of pressure P ur  with time, and plot  320  shows the change of pressure P ur  with time. As shown in  FIG.  3 A , starting from a pressure declining point t d  (about 15 second), both P ur  and P dr  start to decline, indicating that the transport from the upstream and downstream gas reservoirs to an interior of the core sample is the dominant process while gas is still being transported from the upstream gas reservoir to the downstream reservoir. At a pressure merging point t m    334  (e.g., about 21 second), P ur  and P dr  merge to have the same pressure. P ur  and P dr  continue to decline until the measurement concludes, e.g., at about 55 second, when there is no pressure drop, i.e., P ur =P dr =P ∞ . 
     A time period can be selected within a time range from the pressure declining point t d  to the pressure merging point t m , and measurement data in the time period can be used to determine matrix permeability of the core sample. For example, as illustrated in  FIG.  3   , measurement data  330  within a time period from a starting point t s    332  to the pressure merging point t m    334  can be extracted to determine the matrix permeability of the core sample. 
     Third, matrix permeability of the sample is determined based on the measurement data extracted before the upstream and downstream reservoir pressures merge. Pressure P ms , at an axial center of the core sample can be calculated according to Equation 12 using the data P ur , P dr , and 
                 P   ∞     :           P     m   ⁢   s         =         -     (           V   ur       V   p       ⁢     3   2       +     1   4       )       ⁢     P   ur       -       (           V     d   ⁢   r         V   p       ⁢     3   2       +     1   4       )     ⁢     P   dr       +       3   2     ⁢     (         V   p     +     V     u   ⁢   r       +     V     d   ⁢   r           V   p       )     ⁢       P   ∞     .               
Volume ratios AA and BB are also calculated as follows: AA=V ur /V p  and BB=V dr /V p .
 
     An initial time point to is selected from the extracted time period [t s , t m ), where p ur &gt;p ms  and P dr &gt;p ms . The initial time point to can be the starting point t s , but not the merging point t m . Two arrays G(t) and f(t) that change with time t can be constructed according to Equations 20 and 21: 
                 G   ⁡   (   t   )     =         (       5   ⁢   AA     +     1   2       )     ⁢     (         P   ur     (   t   )     -       P   ur     (     t   0     )       )       +       (       5   ⁢   BB     +     1   2       )     ⁢     (         P   dr     (   t   )     -       P   dr     (     t   0     )       )           ,         
and
 
     
       
         
           
             
               f 
               ⁡ 
               ( 
               t 
               ) 
             
             = 
             
               
                 ∫ 
                 
                   t 
                   0 
                 
                 t 
               
               
                 
                   1 
                   
                     μ 
                     ⁢ 
                     Φ 
                     ⁢ 
                     
                       C 
                       g 
                     
                     ⁢ 
                     
                       L 
                       2 
                     
                   
                 
                 ⁢ 
                 4 
                 ⁢ 
                 8 
                 ⁢ 
                 
                   ( 
                   
                     
                       
                         ( 
                         
                           AA 
                           + 
                           
                             1 
                             2 
                           
                         
                         ) 
                       
                       ⁢ 
                       
                         P 
                         ur 
                       
                     
                     + 
                     
                       
                         ( 
                         
                           BB 
                           + 
                           
                             1 
                             2 
                           
                         
                         ) 
                       
                       ⁢ 
                       
                         P 
                         dr 
                       
                     
                     - 
                     
                       
                         ( 
                         
                           1 
                           + 
                           
                             A 
                             ⁢ 
                             A 
                           
                           + 
                           
                             B 
                             ⁢ 
                             B 
                           
                         
                         ) 
                       
                       ⁢ 
                       
                         P 
                         ∞ 
                       
                     
                   
                   ) 
                 
                 ⁢ 
                 d 
                 ⁢ 
                 
                   t 
                   . 
                 
               
             
           
         
       
     
     The matrix permeability k is determined based on the two arrays G(t) and f(t) according to Equation 19: G(t)=−k f(t)+constant. In some examples, linear regression is performed on the two arrays G(t) and f(t) to get the slope −k, that is the negative of the matrix permeability. In some examples, f(t) and G(t) are plotted. The matrix permeability is determined from the slope −k.  FIG.  3 B  is a diagram  350  showing the plot between f(t) and G(t) based on measurement data  360  selected from the measurement data  330  of  FIG.  3 A . It is shown that the two arrays f(t) and G(t) have a linear relationship which can be fitted by a fitted line  370 . The slope of the fitted line  370  is −k. Thus, the matrix permeability k of the core sample can be determined from the slope, e.g., 360 nD. 
     The fracture permeability of the core sample can be also determined using the measurement data  330  of  FIG.  3 A . The fracture permeability can be obtained using the formula:
 
 P   ur ( t )− P   dr ( t )= f (ϕ, L,A,V   bulk   ,c,μ,V   p   ,V   ur   ,V   dr   ,k,t )  (22),
 
where P ur (t) is the pressure of the upstream side, P dr (t) is the pressure downstream side, t is time, k the fracture permeability, ϕ the sample&#39;s total porosity, L is the sample length, and A is the cross section area of the sample, V bulk  the bulk volume of the sample (V bulk =LA; and
 
                 ϕ   =       V   p       V   bulk         )     ,         
c is the compressibility of the fluid, μ is the gas viscosity, V p  the pore volume of the sample, V ur  is the upstream reservoir volume, and V dr  is the downstream reservoir volume. Thus, within the same measurement, both matrix permeability and fracture permeability of the core sample can be obtained. The obtained permeability values can be used in reservoir simulation models to predict well productions from the corresponding unconventional reservoirs.
 
     Example Process 
       FIG.  4    is a flowchart of an example process  400  of determining matrix permeability of a subsurface formation. The process  400  can be performed by a system including a measurement system, e.g., the measurement system  100  of  FIG.  1   , and a computing system. In some implementations, the computing system is separate from the measurement system. In some implementations, the measurement system includes the computing system. The computing system can include at least one processor, and at least one non-transitory machine readable storage medium coupled to the at least one processor having machine-executable instructions stored thereon that, when executed by the at least one processor, cause the at least one processor to perform one or more operations. 
     At  402 , a sample of the subsurface formation is positioned in a measurement cell of the measurement system. The subsurface formation can include at least one of shale, limestone, siltstone, or sandstone. The sample of the subsurface formation can be a fractured sample. The sample can be prepared using the subsurface formation. For example, the sample can be prepared to be a core plug sample, e.g., the sample  150  of  FIGS.  1 ,  2 A and  2 B . The sample can be sleeved so that the sample can be sealed. For example, the core sample can be wrapped in layers of plastic and metal sheet and put between two end pieces with metal lines. The core sample can have a cylindrical shape. A dimension of the core sample can be measured and recorded. The dimension can include a length L of the core sample and a cross-sectional area A of the sample. The measurement cell can be the measurement cell  130  of  FIG.  1   . 
     At  404 , an inlet and an outlet of the sample are fluidly connected to an upstream reservoir and a downstream reservoir, respectively. The upstream reservoir and the downstream reservoir can be the upstream reservoir  110  and the downstream reservoir  120  of  FIG.  1   . 
     A confining pump can be used to introduce a confining fluid, e.g., water, into the measurement cell. Once the measurement cell is full, the confining pump continues to flow the confining fluid into the measurement cell to ensure that there is no gas remaining inside the measurement cell. After the confining of the measurement cell is locked in, an effective stress pressure is immediately applied to the sample, e.g., from all directions of the sample, to reduce the possibility of water leakage into the sleeve of the core sample. The effective stress pressure can be about 500 psi. After the effective stress pressure is applied, a pore pressure and a confining pressure can be adjusted to any desired values. 
     At  406 , a fluid is flowed through the sample from the upstream reservoir to the downstream reservoir. The fluid can include one or more gases from a group of gases including methane, argon, nitrogen, carbon dioxide, helium, ethane, and propane. Both upstream and downstream pumps that are connected to the sample through the upstream and downstream reservoirs can be used to slowly impose a pore pressure to the core sample. The pore pressure can be smaller than the confining pressure The pore pressure can be about 2,500 psi. The use of the pore pressure is to minimize the impact of diffusion on permeability measurements. As the gas from the upstream and downstream reservoirs is introduced into the core sample, the confining pressure in the measurement cell is increased simultaneously to maintain a pressure difference between the confining pressure and the pore pressure to be the effective stress pressure. Once the confining pump and the upstream and downstream pumps are stable at a desired pressure, a pressure pulse decay permeability (PDP) measurement can be carried out. 
     At  408 , changes of an upstream pressure and a downstream pressure in a measurement time period are measured. The upstream pressure can be associated with the upstream reservoir. In some examples, the upstream pressure can be a pressure in the upstream reservoir, which can be measured by a pressure gauge, e.g., the pressure gauge  160 , that is inserted into the upstream reservoir. In some examples, the upstream pressure can be a pressure before an inlet of the sample measured by a pressure gauge, e.g., the pressure gauge  162  between the valve  102  and the inlet of the sample of  FIG.  1   . The downstream pressure can be associated with the downstream reservoir. In some examples, the downstream pressure can be a pressure after the outlet of the sample measured by a pressure gauge, e.g., the pressure gauge  164  after the outlet of the sample before the valve  104  of  FIG.  1   . 
     The upstream pressure and the downstream pressure can be measured by carrying out a pressure pulse decay permeability (PDP) measurement. The measurement data shows changes of the upstream pressure and the downstream pressure in the measurement time period, as illustrated in  FIG.  3 A . 
     At  410 , the matrix permeability of the subsurface formation is determined based on measurement data before the upstream pressure and the downstream pressure merge at a merging time point. The computing system can process the measurement data to determine the matrix permeability of the subsurface formation. 
     The computing system can determine an initial time point when the upstream pressure and the downstream pressure are larger than a central pressure at an axial center of the sample. The initial time point can be a time point no earlier than a declining time point when both the upstream pressure and the downstream pressure start to decline, e.g., the declining time point td of  FIG.  3 A . The computing system can determine the central pressure at the axial center of the sample using the upstream pressure and the downstream pressure according to a formula: 
                 P     m   ⁢   s       =         -     (           V     u   ⁢   r         V   p       ⁢     3   2       +     1   4       )       ⁢     P   ur       -       (           V     d   ⁢   r         V   p       ⁢     3   2       +     1   4       )     ⁢     P   dr       +       3   2     ⁢     (         V   p     +     V     u   ⁢   r       +     V     d   ⁢   r           V   p       )     ⁢     P   ∞           ,         
where P ms , P ur , and P dr  represent the central pressure, the upstream pressure, the downstream pressure, respectively, V ur  and V dr  represent volumes of the upstream reservoir and the downstream reservoir, respectively, V p  represents a total volume of pores in the sample, and P ∞  represents an equilibrium pressure at an end of the measurement time period when P ur  and P dr  are stable.
 
     The computing system can extract measurement data within a particular time range from the initial time point to an ending time point no later than the merging time point, e.g., the measurement data  330  of  FIG.  3 A , and determine the matrix permeability of the subsurface formation using the extracted measurement data within the particular time range. In some examples, the computing system can determine determining the matrix permeability based on data arrays G(t) and f(t) obtained from the extracted measurement data within the particular time range, where: 
                 G   ⁡   (   t   )     =         (       5   ⁢   AA     +     1   2       )     ⁢     (         P   ur     (   t   )     -       P   ur     (     t   0     )       )       +       (       5   ⁢   BB     +     1   2       )     ⁢     (         P   dr     (   t   )     -       P   dr     (     t   0     )       )           ,         
and
 
                 f   ⁡   (   t   )     =       ∫     t   0     t         1     μ   ⁢   Φ   ⁢     C   g     ⁢     L   2         ⁢   4   ⁢   8   ⁢     (         (     AA   +     1   2       )     ⁢     P   ur       +       (     BB   +     1   2       )     ⁢     P   dr       -       (     1   +     A   ⁢   A     +     B   ⁢   B       )     ⁢     P   ∞         )     ⁢   d   ⁢   t         ,         
where AA=V ur /V p , BB=V dr /V p  V ur  and V dr  represent volumes of the upstream reservoir and the downstream reservoir, respectively, V p  represents a total volume of pores in the sample, t 0  represents the starting time point, t represents a time point within the particular time range, μ represents a viscosity of the fluid, Φ represents a porosity of the sample, L represents a length of the sample, and C g  represents a compressibility of the fluid, P ur  and P dr  represent the upstream pressure and the downstream pressure, respectively, and P ∞  represents an equilibrium pressure at an end of the measurement time period when P ur  and P dr  are stable after P ur  and P dr  merge.
 
     The computing system can determine the matrix permeability based on the data arrays permeability according to a formula:
 
 G ( t )=− kf ( t )+ C,  
 
where k represents the matrix permeability, and C represents a constant. The formula can be derived based on mass balance and Darcy law. Then the computing system can perform linear regression or linear fitting on the data arrays to get a slope and determine the matrix permeability to be a negative of the slope.
 
     In some cases, the computing system can further determine a fracture permeability of the subsurface formation based on the measurement data before the upstream pressure and the downstream pressure merge. Thus, within the same measurement, both matrix permeability and fracture permeability of the core sample can be obtained. The obtained permeability values can be used in reservoir simulation models to predict well productions from the corresponding unconventional reservoirs. 
     Implementations of the subject matter and the functional operations described in this specification can be implemented in digital electronic circuitry, in tangibly embodied computer software or firmware, in computer hardware, including the structures disclosed in this specification and their structural equivalents, or in combinations of one or more of them. Implementations of the subject matter described in this specification can be implemented as one or more computer programs, such as, one or more modules of computer program instructions encoded on a tangible, non-transitory computer-storage medium for execution by, or to control the operation of, data processing apparatus. Alternatively or in addition, the program instructions can be encoded on an artificially generated propagated signal, such as, a machine-generated electrical, optical, or electromagnetic signal that is generated to encode information for transmission to suitable receiver apparatus for execution by a data processing apparatus. The computer-storage medium can be a machine-readable storage device, a machine-readable storage substrate, a random or serial access memory device, or a combination of one or more of them. 
     The terms “data processing apparatus,” “computer,” or “electronic computer device” (or equivalent as understood by one of ordinary skill in the art) refer to data processing hardware and encompass all kinds of apparatus, devices, and machines for processing data, including by way of example, a programmable processor, a computer, or multiple processors or computers. The apparatus can also be or further include special purpose logic circuitry, for example, a central processing unit (CPU), an FPGA (field programmable gate array), or an ASIC (application-specific integrated circuit). In some implementations, the data processing apparatus and special purpose logic circuitry may be hardware-based and software-based. The apparatus can optionally include code that creates an execution environment for computer programs, for example, code that constitutes processor firmware, a protocol stack, a database management system, an operating system, or a combination of one or more of them. The present specification contemplates the use of data processing apparatuses with or without conventional operating systems. 
     A computer program, which may also be referred to or described as a program, software, a software application, a module, a software module, a script, or code, can be written in any form of programming language, including compiled or interpreted languages, or declarative or procedural languages, and it can be deployed in any form, including as a stand-alone program or as a module, component, subroutine, or other unit suitable for use in a computing environment. A computer program may, but need not, correspond to a file in a file system. A program can be stored in a portion of a file that holds other programs or data, for example, one or more scripts stored in a markup language document, in a single file dedicated to the program in question, or in multiple coordinated files, for example, files that store one or more modules, sub-programs, or portions of code. A computer program can be deployed to be executed on one computer or on multiple computers that are located at one site or distributed across multiple sites and interconnected by a communication network. While portions of the programs illustrated in the various figures are shown as individual modules that implement the various features and functionality through various objects, methods, or other processes, the programs may instead include a number of sub-modules, third-party services, components, libraries, and such, as appropriate. Conversely, the features and functionality of various components can be combined into single components as appropriate. 
     The processes and logic flows described in this specification can be performed by one or more programmable computers executing one or more computer programs to perform functions by operating on input data and generating output. The processes and logic flows can also be performed by, and apparatus can also be implemented as, special purpose logic circuitry, such as, a CPU, an FPGA, or an ASIC. 
     Computers suitable for the execution of a computer program can be based on general or special purpose microprocessors, both, or any other kind of CPU. Generally, a CPU will receive instructions and data from a read-only memory (ROM) or a random access memory (RAM) or both. The essential elements of a computer are a CPU for performing or executing instructions and one or more memory devices for storing instructions and data. Generally, a computer will also include, or be operatively coupled to, receive data from or transfer data to, or both, one or more mass storage devices for storing data, for example, magnetic, magneto-optical disks, or optical disks. However, a computer need not have such devices. Moreover, a computer can be embedded in another device, for example, a mobile telephone, a personal digital assistant (PDA), a mobile audio or video player, a game console, a global positioning system (GPS) receiver, or a portable storage device, for example, a universal serial bus (USB) flash drive, to name just a few. 
     Computer-readable media (transitory or non-transitory, as appropriate) suitable for storing computer program instructions and data include all forms of non-volatile memory, media and memory devices, including by way of example semiconductor memory devices, for example, erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), and flash memory devices; magnetic disks, for example, internal hard disks or removable disks; magneto-optical disks; and CD-ROM, DVD-R, DVD-RAM, and DVD-ROM disks. The memory may store various objects or data, including caches, classes, frameworks, applications, backup data, jobs, web pages, web page templates, database tables, repositories storing business and dynamic information, and any other appropriate information including any parameters, variables, algorithms, instructions, rules, constraints, or references thereto. Additionally, the memory may include any other appropriate data, such as logs, policies, security or access data, reporting files, as well as others. The processor and the memory can be supplemented by, or incorporated in, special purpose logic circuitry. 
     To provide for interaction with a user, implementations of the subject matter described in this specification can be implemented on a computer having a display device, for example, a cathode ray tube (CRT), liquid crystal display (LCD), light emitting diode (LED), or plasma monitor, for displaying information to the user and a keyboard and a pointing device, for example, a mouse, trackball, or trackpad by which the user can provide input to the computer. Input may also be provided to the computer using a touchscreen, such as a tablet computer surface with pressure sensitivity, a multi-touch screen using capacitive or electric sensing, or other type of touchscreen. Other kinds of devices can be used to provide for interaction with a user as well; for example, feedback provided to the user can be any form of sensory feedback, for example, visual feedback, auditory feedback, or tactile feedback; and input from the user can be received in any form, including acoustic, speech, or tactile input. In addition, a computer can interact with a user by sending documents to and receiving documents from a device that is used by the user; for example, by sending web pages to a web browser on a user&#39;s client device in response to requests received from the web browser. 
     The term “graphical user interface,” or “GUI,” may be used in the singular or the plural to describe one or more graphical user interfaces and each of the displays of a particular graphical user interface. Therefore, a GUI may represent any graphical user interface, including but not limited to, a web browser, a touch screen, or a command line interface (CLI) that processes information and efficiently presents the information results to the user. In general, a GUI may include multiple user interface (UI) elements, some or all associated with a web browser, such as interactive fields, pull-down lists, and buttons operable by the business suite user. These and other UI elements may be related to or represent the functions of the web browser. 
     Implementations of the subject matter described in this specification can be implemented in a computing system that includes a back-end component, for example, as a data server, or that includes a middleware component, for example, an application server, or that includes a front-end component, for example, a client computer having a graphical user interface or a web browser through which a user can interact with an implementation of the subject matter described in this specification, or any combination of one or more such back-end, middleware, or front-end components. The components of the system can be interconnected by any form or medium of wireline or wireless digital data communication, for example, a communication network. Examples of communication networks include a local area network (LAN), a radio access network (RAN), a metropolitan area network (MAN), a wide area network (WAN), worldwide interoperability for microwave access (WIMAX), a wireless local area network (WLAN) using, for example, 902.11 a/b/g/n and 902.20, all or a portion of the Internet, and any other communication system or systems at one or more locations. The network may communicate with, for example, internet protocol (IP) packets, frame relay frames, asynchronous transfer mode (ATM) cells, voice, video, data, or other suitable information between network addresses. 
     The computing system can include clients and servers. A client and server are generally remote from each other and typically interact through a communication network. The relationship of client and server arises by virtue of computer programs running on the respective computers and having a client-server relationship to each other. 
     In some implementations, any or all of the components of the computing system, both hardware and software, may interface with each other or the interface using an application programming interface (API) or a service layer. The API may include specifications for routines, data structures, and object classes. The API may be either computer language-independent or -dependent and refer to a complete interface, a single function, or even a set of APIs. The service layer provides software services to the computing system. The functionality of the various components of the computing system may be accessible for all service consumers via this service layer. Software services provide reusable, defined business functionalities through a defined interface. For example, the interface may be software written in any suitable language providing data in any suitable format. The API and service layer may be an integral or a stand-alone component in relation to other components of the computing system. Moreover, any or all parts of the service layer may be implemented as child or sub-modules of another software module, enterprise application, or hardware module without departing from the scope of this specification. 
     While this specification contains many specific implementation details, these should not be construed as limitations on the scope of any invention or on the scope of what may be claimed, but rather as descriptions of features that may be specific to particular implementations of particular inventions. Certain features that are described in this specification in the context of separate implementations can also be implemented in combination in a single implementation. Conversely, various features that are described in the context of a single implementation can also be implemented in multiple implementations separately or in any suitable sub-combination. Moreover, although features may be described as acting in certain combinations and even initially claimed as such, one or more features from a claimed combination can in some cases be excised from the combination, and the claimed combination may be directed to a sub-combination or variation of a sub-combination. 
     Particular implementations of the subject matter have been described. Other implementations, alterations, and permutations of the described implementations are within the scope of the following claims as will be apparent to those skilled in the art. While operations are depicted in the drawings or claims in a particular order, this should not be understood as requiring that such operations be performed in the particular order shown or in sequential order, or that all illustrated operations be performed (some operations may be considered optional), to achieve desirable results. In certain circumstances, multitasking or parallel processing may be advantageous and performed as deemed appropriate. 
     Moreover, the separation or integration of various system modules and components in the implementations described earlier should not be understood as requiring such separation or integration in all implementations, and it should be understood that the described program components and systems can generally be integrated together in a single software product or packaged into multiple software products. 
     Accordingly, the earlier provided description of example implementations does not define or constrain this specification. Other changes, substitutions, and alterations are also possible without departing from the spirit and scope of this specification.