Abstract:
The invention is a method of determining all the components of an absolute permeability tensor of a porous medium sample from permeability measurements, obtained for example by placing the sample in a permeameter which is useful for fast determination of permeability anisotropies of rocks and detention of internal heterogeneities. A pressure difference ΔP is applied between the inlet and outlet faces of a laboratory rock sample, with zero flow conditions on the edges parallel to the mean flow obtained by confining it in a sheath under pressure. Starting from conventional permeability measurements in three directions and from the measurements of the two components of the viscous forces transverse to the sample, permeability tensor k can be “inverted” by numerical solution of the corresponding boundary-value problem. Since these quantities can be alternatively obtained from numerical solutions of the same flow in a heterogeneous medium, the method can also be used as an upscaling tool in a reservoir simulator.

Description:
[0001]     A method of determining a permeability tensor of a rock sample from measurements allowing estimation of, with sufficient accuracy, the mean pressures around a sample placed in a permeameter suited to the measurements thereof is described for example in the following publication: 
        Renard P. “ Laboratory Determination of the Full Permeability Tensor ” JGR 106, B11 2001 pp 26 443-26 451.        
 
         [0003]     Another known method of determining a permeability tensor wherein flows which are made deliberately tortuous by changes in the boundary conditions of the injection faces are created in the sample and wherein an inverse problem is solved is for example described in the following publication: 
        Bernabé Y.: “ On the Measurement of Permeability in Anisotropic Rocks ” in “Fault Mechanics and Transport Properties of Rocks” edited by B. Evans and T F Wong pp 147-167, Academic San Diego 1992.        
 
         [0005]     In order to be implemented, the known methods require substantial changes in the flow conditions, which makes them costly and difficult to be applied in practice.  
         [0006]     Within the context of upscaling, well-known to reservoir engineers, laboratory experiments are replaced by the results of a “fine” numerical simulation on the heterogeneous medium, in order to replace it by an equivalent homogeneous medium whose permeability tensor is for example denoted by K eff . In the case of an anisotropic medium, the fine reference simulation is generally performed by considering periodic boundary conditions. Most authors wrongly consider that these are the only conditions allowing obtaining all the elements of tensor K eff .  
         [0007]     As shown in the description hereafter, a suitable interpretation of a numerical simulation of the same permeameter can provide this information. The boundary conditions applied in a permeameter can be more realistic in the case where the large-scale flow is constrained by clay barriers. Furthermore, the comparison between the two resulting tensors can provide information on the scale of the Representative Elementary Volume (REV).  
       SUMMARY OF THE INVENTION  
       [0008]     A method according to the invention allows determination, from known permeability measurements, all the components of an absolute permeability tensor (K eq, ) of a sample of an anisotropic porous medium. The method of the invention essentially comprises the followingsteps: 
        measuring or calculating, by applying two liquid flows successively in two (or three) orthogonal directions (Ox, Oy) (and Oz) through the sample and under an imposed pressure gradient, quantities the liquid flow rates and the transverse forces generated in the sample by passage of the liquid, the quantities being brought back to unitary pressure gradients; and     using a simulator to determine step by step, from an a priori numerical solution of a confined flow in a homogeneous medium, components of the equivalent tensor (k) by means of any inverse method so that corresponding quantities substantially adjust to the measured or calculated quantities.        
 
         [0011]     According to a first implementation mode, the quantities are measured by means of a permeameter and by means for measuring transverse forces applied to the sample in response to flows imposed through the sample.  
         [0012]     The transverse forces applied to the sample in response to flows imposed therethrough can be measured for example by measuring weight variations of the sample in the permeameter when the direction of circulation of the liquid flows imposed through the sample is alternated. The other component can be measured by repeating this procedure after a 90° rotation of the sample about an axis parallel to the flow.  
         [0013]     According to a second implementation mode, the same quantities can be determined by simulation of flows in a medium, from known or estimated permeability data such as permeability maps provided, for example, by geologists in geophysical applications.  
         [0014]     Implementation of the method according to the invention is advantageous in that it allows to simplify and to accelerate determination of the permeability tensor of samples. 
     
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0015]     Other features and advantages of the method according to the invention will be clear from reading the description hereafter of a non limitative example, with reference to the accompanying drawings wherein  
         [0016]      FIG. 1  shows a flowchart allowing iterative determination of the equivalent tensor K eq (n) of a rock sample, from various flow rate and force measurements performed by confining it in a permeameter;  
         [0017]     FIGS.  2  to  4  illustrate various simulation examples showing the relevance of the method, which will be discussed in the course of the description; and  
         [0018]      FIG. 5  shows a diagrammatic example of a permeameter suited for measurements on a rock sample. 
     
    
     DETAILED DESCRIPTION OF THE INVENTION  
       [0019]     The permeability measurements useful for implementation of the method described hereafter can be obtained by means of a permeameter of a well-known type. Such a permeameter comprises for example ( FIG. 5 ) a containment cell  1  for the sample. The lateral wall of the cell is a deformable sheath. At a first end, cell  1  is connected to a circuit  2  for injecting a fluid at an injection pressure P 1 . At the opposite end thereof, the cell is connected to a circuit  3  for discharge of the fluid coming from the sample maintained at a pressure P 2 . Cell  1  is placed in an outer enclosure  4  connected to a device  5  for applying a confining pressure Pc to the sample, so that the flow of liquid on laterally surrounding the sample is zero (impermeable edges). The opposite terminal parts of cell  1  are such that pressures P 1  and P 2  respectively exerted on the upstream and downstream faces of the sample are uniform. The quantity ΔP x,x equalsthequantityP 2 −P 1 .  
         [0020]     In order to simplify the presentation, a 2D porous medium is considered which is, square, of side L. The method transposed to 3D cases with any proportions is described hereafter. In this part, the porous rock sample is considered to be homogeneous but with a permeability tensor k whose proper axes do not necessarily coincide with the axes of the permeameter where it is placed ( FIG. 1 ).  
         [0021]     With reference to  FIG. 1 , the liquid flow is first applied to the sample oriented in a first direction A in relation to the axis of the cell and a series of measurements described below is performed. The experiment is then repeated after orienting the sample in a direction B orthogonal to the first direction, and the same measurements are performed according to the same protocol.  
         [0022]     On the faces of the sample parallel to the imposed flow, the standard Neumann zero flow conditions are present, that is 
 
 n.k.∇p ( r )=0   (1) 
 
 where r is the spatial position vector of components x, y (or x 1 , x 2 ), n is the unitary vector orthogonal to the boundary at the point considered, and ∇p(r) is the pressure gradient. Component by component, a formula is written as follows:  
             ∑     i   ,     j   =   1     ,   D       ⁢       n   i     ·     k   ij     ·       ∂     p   ⁡     (   r   )           ∂     x   j             =   0     ,       
 
 or more concisely, by adopting from now on the Einstein summation convention on the repeated indices:  
                 n   i     ⁢     k   ij     ⁢       ∂     p   ⁡     (   r   )           ∂     x   j           =   0           (   2   )             
 
         [0023]     Pressure equation: the pressure equation is given conventionally, the point (·) corresponding to tensor contractions  
               ∇     ·     (       k   μ     ·     ∇     p   ⁡     (   r   )           )         =   0           (   3   )             
 
 i.e., in extenso, still with the Einstein convention:  
                 ∂     ∂     x   i         ⁡     [       k   ij     ⁢       ∂     p   ⁡     (   r   )           ∂     x   j           ]       =   0           (   4   )             
 
         [0024]     The pressure equation and the boundary conditions define a problem denoted by problem (P), well set and having a single solution. In the particular case where direction x coincides with a proper axis of k, the linear solution  
         p   ⁡     (     x   ,   y   ,   z     )       =       P   1     +     Δ   ⁢           ⁢       P     x   ,   x       ·     x   L               
 
 is the single solution to the problem. We then obtain, in this case, the conventional relation  
         Q     x   ,   x       =       -     k     x   ,   x         ⁢       Δ   ⁢           ⁢     P     x   ,   x         μ           
 
 from which k x,x  is inferred if Q x,x  and ΔP x,x  have been measured simultaneously. 
 
         [0025]     In the general anisotropic case, solution p(x, y, z) has no simple structure and there is no known analytical solution. Therefore the assumption of that the situation presented above is applicable and in practice estimations of the proper values of tensor k are provided.  
         [0026]     It is now shown that, provided that some additional quantities are measured, it is possible to determine all the elements k ij  of tensor k. Therefore observe that equation (4) can be rewritten in the form as follows:  
                 ∫   V     ⁢         ⅆ   D     ⁢   r     ⁢       ∂     ∂     x   i         ⁡     [       k   ij     ⁢       ∂     p   ⁡     (   r   )           ∂     x   j           ]       ⁢     f   ⁡     (   r   )           =     0   ⁢     ∀     f   ⁡     (   r   )                   (   5   )             
 
 where r is the D-component position vector and d D r the integration element. By selecting then function f(r) suitably (f(r)=x, y then z) in 3D, by integrating by parts the previous equality, the following set of equalities in 2D is obtained:  
           μ   ⁢           ⁢     q   xx       =         k   xx     ⁢     1     L   2       ⁢       ∫   0   L     ⁢     ⅆ     y   ⁡     (       p   ⁡     (       x   =   L     ,   y     )       -     p   ⁡     (       x   =   0     ,   y     )         )             +     
     ⁢           ⁢       k   xy     ⁢     1     L   2       ⁢       ∫   0   L     ⁢     ⅆ     x   ⁡     (       p   ⁡     (     x   ,   y     )       =   L     )             -     p   ⁡     (     x   ,     y   =   0       )           )       
 
 i.e., with more condensed notations:  
               μ   ⁢           ⁢     q   xx       =         k   xx     ⁢       {       ∂   p       ∂   x       }     x       +       k   xy     ⁢       {       ∂   p       ∂   y       }     x                 (   6   )                 μ   ⁢           ⁢   δ   ⁢           ⁢     y   x     ⁢     q   xx       =         k   yx     ⁢       {       ∂   p       ∂   x       }     x       +       k   yy     ⁢       {       ∂   p       ∂   y       }     x                 (   7   )             
 
 where q xx  corresponds to the flow rate normalized by the volume of sample q xx =Q x,x  /V. The {. . . } x  designate a spatial mean and index x represents here the direction of the pressure gradient imposed. Quantities δy x  and δx y  are defined by:  
         δ   ⁢           ⁢     y   x       =           ∫   0   L     ⁢         yk   xx     ⁡     [           ∂   p       ∂   x       ⁢     (     L   ,   y     )       -         ∂   p       ∂   x       ⁢     (     0   ,   y     )         ]       ⁢     ⅆ   y             ∫   0   L     ⁢       k   xx     ⁢       ∂   p       ∂   x       ⁢     (     0   ,   y     )     ⁢     ⅆ   y           ⁢     
     ⁢           =         ∫   0   L     ⁢         yk   xx     ⁡     [           ∂   p       ∂   x       ⁢     (     L   ,   y     )       -         ∂   p       ∂   x       ⁢     (     0   ,   y     )         ]       ⁢     ⅆ   y           Vq   xx             
 
         [0027]     This quantity corresponds to the difference of coordinates y of the barycenters of the flow rates in local x on the downstream and upstream faces of the sample during a flow imposed in direction x. Similarly, the following relationship is applicable:  
         δ   ⁢           ⁢     x   y       =           ∫   0   L     ⁢         xk   yy     ⁡     [           ∂   p       ∂   y       ⁢     (     x   ,   L     )       -         ∂   p       ∂   y       ⁢     (     x   ,   0     )         ]       ⁢     ⅆ   x             ∫   0   L     ⁢       k   yy     ⁢       ∂   p       ∂   y       ⁢     (     x   ,   0     )     ⁢     ⅆ   x           ⁢     
     ⁢           =         ∫   0   L     ⁢         xk   yy     ⁡     [           ∂   p       ∂   y       ⁢     (     x   ,   L     )       -         ∂   p       ∂   x       ⁢     (     x   ,   0     )         ]       ⁢     ⅆ   x           Vq   yy             
 
         [0028]     Similar formulas can be shown by means of the same 3D technique, introducing the three flow rates q xx , q yy , q zz , the six quantities (δy x  δz x  δx y  δz y  δx z  and δy z ), and the six additional other quantities  
         (         {       ∂   p       ∂   y       }     x     ,       {       ∂   p       ∂   z       }     x     ,       {       ∂   p       ∂   x       }     y     ,       {       ∂   p       ∂   z       }     y     ,       {       ∂   p       ∂   y       }     z     ,       {       ∂   p       ∂   y       }     z       )     .       
 
         [0029]     The Dirichiet conditions on the upstream and downstream faces of course provide the relation  
           {       ∂   p       ∂   x       }     x     =         Δ   ⁢           ⁢     P   xx       L     .         
 
 Considering an experiment carried out in another orthogonal direction, the following relationships are obtained:  
               μ   ⁢           ⁢     q   yy       =         k   yx     ⁢       {       ∂   p       ∂   x       }     y       +       k   yy     ⁢       {       ∂   p       ∂   y       }     y                 (   8   )                 μ   ⁢           ⁢   δ   ⁢           ⁢     x   y     ⁢     q   yy       =         k   xx     ⁢       {       ∂   p       ∂   x       }     y       +       k   xy     ⁢       {       ∂   p       ∂   y       }     y                 (   9   )             
 
         [0030]     In order to simplify the notations, unitary imposed pressure gradients  
           {       ∂   p       ∂   x       }       x   ⁢           ⁢   ou   ⁢           ⁢   y       =         Δ   ⁢           ⁢     P     xx   ⁢           ⁢   ou   ⁢           ⁢   yy         L     =   1         
 
 are utilized. By inverting these equalities, in this case the following relations are obtained:  
               k   xx     =         q   xx     -       q   yy     ⁢   δ   ⁢           ⁢     x   y     ⁢       {       ∂   p       ∂   y       }     x           1   -         {       ∂   p       ∂   x       }     y     ⁢       {       ∂   p       ∂   y       }     x                   (   10   )                 k   yy     =         q   yy     -       q   xx     ⁢   δ   ⁢           ⁢     y   x     ⁢       {       ∂   p       ∂   x       }     y           1   -         {       ∂   p       ∂   y       }     x     ⁢       {       ∂   p       ∂   x       }     y                   (   11   )                 k   xy     =         -       {       ∂   p       ∂   x       }     y       ⁢     k   xx       +     δ   ⁢           ⁢     x   y     ⁢     q   yy                 (   12   )                 k   yx     =         -       {       ∂   p       ∂   y       }     x       ⁢     k   yy       +     δ   ⁢           ⁢     y   x     ⁢     q   xx                 (   14   )             
 
         [0031]     These relations can be grouped together in a symbolic form  
             k   =       F   ⁡     (       q   xx     ,     q   yy     ,       {       ∂   p       ∂   y       }     x     ,       {       ∂   p       ∂   x       }     y     ,     δ   ⁢           ⁢     y   x       ,     δ   ⁢           ⁢     x   y         )       .             (   15   )             
 
         [0032]     To sum up, in the 2D homogeneous case, by measuring quantities  
         q   xx     ,     q   yy     ,       {       ∂   p       ∂   x       }     y     ,       {       ∂   p       ∂   y       }     x     ,     δ   ⁢           ⁢     y   x       ,     δ   ⁢           ⁢     x   y       ,       
 
 tensor k can be found directly. Since the tensor is defined by three 2D numbers and six 3D numbers (but, in this case, fifteen additional quantities are determined), this means that there are implicit relations between quantities  
         q   xx     ,     q   yy     ,       {       ∂   p       ∂   y       }     x     ,       {       ∂   p       ∂   x       }     y     ,     δ   ⁢           ⁢     y   x       ,     δ   ⁢           ⁢       x   y     .           
 
 In theory, if at least three values are maintained, it should still be possible to find k. 
 
         [0033]     The following subset is maintained:  
               q   xx     ,     q   yy     ,       {       ∂   p       ∂   y       }     x     ,       {       ∂   p       ∂   x       }     y             (   16   )             
 
         [0034]     The main interest of these data is to have an immediate physical sense so that they can be directly measured: q xx , q yy  correspond to the usual flow rates (brought back to the volume of the experiment). Quantities  
           {       ∂   p       ∂   y       }     x     ,       {       ∂   p       ∂   x       }     y         
 
 are measurable because:  
           {       ∂   p       ∂   y       }     x     =       1     L   2       ⁢       ∫   0   L     ⁢     ⅆ     x   ⁡     (       p   ⁡     (     x   ,     y   =   L       )       -     p   ⁡     (     x   ,     y   =   0       )         )                 
 
 is merely a resultant of the transverse force per volume unit applied to the sample. Furthermore, by using the general relation:  
             1   V     ⁢       ∫   V     ⁢         ⅆ   D     ⁢   r     ⁢     ∇     p   ⁡     (   r   )               =       1   V     ⁢       ∫     ∂   V       ⁢         ⅆ     D   -   1       ⁢     rp   ⁡     (   r   )         ⁢   n           ,       
 
 it is seen that this interpretation of the mean pressure gradient in terms of surface force exerted on the sample remains. 
 
         [0035]     All these quantities or forces can be measured directly in the laboratory. To measure the transverse forces, it is possible to use, for example, dynamometers or strain gages applied locally against the lateral walls of the sample in its sheath (See  FIG. 5 ).  
         [0036]     To measure the resultant of the transverse forces more conveniently, it is also possible to measure the weight variations of the sample according to whether the transverse forces add to or are subtracted from its static weight. All things being equal, a flow is created through the sample successively in a first direction and in the opposite direction. In the first case, the transverse force adds for example to the weight of the cell and, in the other case, it is deducted therefrom. By subtracting the measurements, the measurement of the transverse force can be directly obtained. For implementation, enclosure  4  (see  FIG. 5 ) is placed horizontally on a weighing device (not shown). The weight variation is measured for two positions oriented 90° apart, corresponding to directions A and B shown by arrows in  FIG. 1 .  
         [0037]     2D Example  
         [0038]     It is assumed that the quadruplet  
         q   xx     ,     q   yy     ,       {       ∂   p       ∂   y       }     x     ,       {       ∂   p       ∂   x       }     y         
 
 of measurements of the flow rates, and transverse forces applied to the sample in the successive positions imposed thereon in the permeameter (four global 2D force and flow rate measurements, nine 3D measurements) is obtained. From now on, the porous media considered can be heterogeneous, that is characterized by a fine map of permeability tensors depending on the position k(r) that can be provided by the geologist or any other technique. 
 
         [0039]     In order to go back to tensor k-K eq  is denoted to identify that it is an effective tensor—from knowledge of  
         q   xx     ,     q   yy     ,       {       ∂   p       ∂   y       }     x     ,       {       ∂   p       ∂   x       }     y     ,       
 
 it is necessary to have a numerical solution to the anisotropic Laplace equation. Any known inversion method can be used. The user can select, for example, an “optimization” method (to minimize the difference between the measured values of  
         q   xx     ,     q   yy     ,       {       ∂   p       ∂   y       }     x     ,       {       ∂   p       ∂   x       }     y         
 
 and the predictions of the numerical model), a trial-and-error method, a dichotomy method, etc. An example is given hereafter. 
 
         [0040]     In order to find the tensor K eq , associated with the equivalent homogeneous medium, a numerical simulator of a well-known type available on the market is used, operating by finite elements or volumes, to solve the 2D (or 3D) boundary-value problem. The procedure, for example, constructs a sequence of tensors K eq (1), . . . , K eq (n), etc., defined by its first term representing a first estimation of K eq  obtained using relation (15):  
             K   eq     ⁡     (   1   )       =     F   ⁡     (       q   xx     ,     q   yy     ,       {       ∂   p       ∂   y       }     x     ,       {       ∂   p       ∂   x       }     y     ,   0   ,   0     )         ,       
 
         [0041]     Since the moments δy x  and δx y  are unknown, because they have not been measured, this can only be an approximation since, as discussed above, there is an implicit relation between the moments and  
         q   xx     ,     q   yy     ,       {       ∂   p       ∂   y       }     x     ,         {       ∂   p       ∂   x       }     y     .         
 
 This first estimation does not take it into account and must therefore be corrected. The estimation of K eq (n) is modified for example by means of the following recurrence relation:  
           K   eq     ⁡     (     n   +   1     )       =       F   ⁡     (       q   xx     ,     q   yy     ,       {       ∂   p       ∂   y       }     x     ,       {       ∂   p       ∂   x       }     y     ,     δ   ⁢           ⁢     x   y   n       ,     δ   ⁢           ⁢     y   x   n         )       .         
 
         [0042]     Here, quantities δx y   n , δy x   n  are calculated by means of the numerical simulator: the boundary-value problem (P) is solved by taking as the permeability tensor the previous estimation K eq (n). The iterations are stopped once the quantities  
       {       q   xx   ′     ,     q   yy   ′     ,       {       ∂   p       ∂   y       }     x   ′     ,       {       ∂   p       ∂   x       }     y   ′       }       
 
 calculated by the simulator with the corresponding estimation K eq (n) are ε from the desired value, ε being a stop criterion set by the user in relation to a norm selected beforehand. The flowchart is given in  FIG. 1 . 
 
         [0043]     The measurement quadruplet used to completely determine tensor k may also have been calculated by means of a “fine simulation” of the experiment on a heterogeneous medium whose equivalent permeability tensor is to be defined by means of this method. The permeability of this medium is known entirely or partly, for example on the basis of a permeability map provided by the geologist. During this simulation, the Darcy&#39;s laws are applied to the flows in the medium and the applicable quadruplet of values is deduced therefrom.  
         [0044]     Validation  
         [0045]     By simulation, the existence of transverse forces are checked, measurable and quite significant in some cases, which appear when a stream flows through a confined medium, due to the presence of internal heterogeneities or to an anisotropy of the medium, and the possibility of identifying all the components of the equivalent permeability tensor by means of simple measurements of forces on samples obtained by means of a permeameter has been validated. Furthermore, this study allows establishing the error made on the behavior of an equivalent medium (defined by one method or another) in relation to the behavior of the real medium to be homogenized. It is on this basis of a better adjustment between the respective behaviors that the method according to the invention was developed.  
         [0046]      FIG. 2B  shows the evolution of the pressure field for an anisotropic homogeneous medium with an anisotropy ratio 10 and where the principal axes of the diagonal tensor ( FIG. 2A ) are oriented 45° to the horizontal (θ=45°).  FIG. 2C  shows the evolution of the ratio of the cumulative pressure difference on the impermeable edges (edges perpendicular to the imposed pressure fields) to the pressure difference applied to the inlet-outlet of the medium (edges parallel to the imposed pressures) for θ (0°≦θ≦180°).  
         [0047]      FIG. 3A  shows a lognormal permeability map (geometric mean Kg=121, σ 2 =2) generated with the known stochastic generation code FFTAM. Locally, the permeability is isotropic (scalar), on the other hand its correlation function is spatially anisotropic, oriented 45°.  FIG. 3B  shows the evolution of the pressure on the upper and lower edges (impermeable edges for a horizontal confined flow) for the real medium and its anisotropic homogeneous equivalent determined by means of the permeameter method.  FIGS. 3C and 3D  respectively show the mappings of the pressure field for the real medium (left) and the equivalent medium (right). After inversion, the components of the permeability tensor of the equivalent medium are estimated: K xx =122.05; k yy =119.86; k xy =k yx =58.08. The proper values of this tensor are: K′ xx =179.05; K′ yy =62.86. The system of principal axes is indeed oriented 44.46° in relation to the imposed direction of flow. This corresponds to the direction of the initial anisotropy.  
         [0048]      FIG. 4A  shows a network of fractures of a fracture volume fraction of 10%, the permeability of the matrix being of the order of 1 Darcy whereas the permeability of the fractures is of the order of 100 Darcy (left).  FIG. 4B  shows the evolution of the pressure on the upper and lower edges (impermeable edges for a horizontal confined flow) for the real medium and its anisotropic homogeneous equivalent determined by means of the permeameter method.  FIGS. 4C and 4D  respectively show the mappings of the pressure field for the real medium and its equivalent. The components of the permeability tensor of the equivalent medium are: K xx =5.79 D; k yy =4.37 D; k xy =k yx =0.26 D. The proper values of this tensor are: K′ xx =5.83 D; K′ yy =4.32 D. The system of principal axes is rotated 9.83° in relation to the system of axis made up of the horizontal and the vertical direction.