Patent Publication Number: US-2013231867-A1

Title: System and method for modeling shale discontinuity

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
     The present application is based on and claims priority to U.S. Provisional Patent Application No. 61/586,370, filed on Jan. 13, 2012. The entire content of which is incorporated herein by reference. 
    
    
     BACKGROUND 
     1. Field 
     The present invention relates generally to modeling reservoir characteristics and more particularly to detection or prediction of shale discontinuities based on reservoir flow data. 
     2. Background 
     A good geological model is critical to model the fluid flow in a waterflood project. Especially when we have limited water resources, the choices of injecting locations can become an important issue in waterflood management. For a layered reservoir with interlayer cross flow, the earlier depletion of the higher permeability layer is predictable. In this situation, it may be useful to change the injection schedule to enhance the oil recovery in the lower permeability layer. To achieve this goal, identifying the shale discontinuities can be helpful in modeling cross layer fluid flow. 
     In the past few decades, many different methods have been developed to estimate the parameters for layered reservoir. Seismic crosshole tomography (Chapman and Pratt, 1992) can provide very high resolution results if enough sensors are installed. Tracer tests (Vasco et al. 1999) are also widely used to investigate the flow property in the reservoir. 
     However, all of the above methods require additional cost and may interrupt the daily operation. Seismic crosshole testing usually takes a long period of time and has difficulty in characterizing the flow property directly. Tracer tests can map the inter-well flow property conveniently, but usually are not repeatable. 
     SUMMARY 
     An aspect of an embodiment of the present invention a method of locating discontinuities in a reservoir. The method includes identifying a plurality of injection and production sites in the reservoir; injecting fluid into selected injection sites in accordance with an injection schedule; monitoring an output at the production sites; determining a time delay between the selected injection sites and the production sites based on the monitored output; and determining a location of discontinuities based on the determined time delay. 
     An aspect of an embodiment of the present invention includes a system for locating discontinuities in a reservoir. The system includes a processor configured to: (a) determine a time delay between injection sites and production sites in a the reservoir based on monitored output at the production sites from injecting fluid to selected injection sites in accordance with an injection schedule; and (b) determine a location of discontinuities based on the determined time delay. 
     Other aspects of embodiments of the present invention include computer readable media encoded with computer executable instructions for performing any of the foregoing methods and/or for controlling any of the foregoing systems. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       Other features described herein will be more readily apparent to those skilled in the art when reading the following detailed description in connection with the accompanying drawings, wherein: 
         FIG. 1  is a flow chart illustrating steps of a method in accordance with an embodiment of the invention; 
         FIG. 2  is a schematic representation of an injector/producer relationship in a multilayer reservoir having a shale discontinuity, according to an embodiment of the present invention; 
         FIG. 3  is a schematic illustration of a flow in an injector/producer relationship illustrating Fermat&#39;s principle in determining the location of shale discontinuities, according to an embodiment of the present invention; 
         FIGS. 4   a - 4   c  illustrate a testing scenario in a simulated application of a method in accordance with an embodiment of the invention; 
         FIG. 5   a  illustrates injection rates, according various embodiments of the present invention 
         FIG. 5   b  illustrates measured production rates corresponding to the injection rates shown in  FIG. 5   a;    
         FIG. 5   c  illustrates an estimated FIR and a CM to determine a time delay constant for the simulated method, according to an embodiment of the present invention; 
         FIG. 6   a  illustrates an initial structure, according to an embodiment of the present invention; 
         FIG. 6   b  illustrates an output of the algorithm, according to an embodiment of the present invention; 
         FIG. 6   c  illustrates the simulated reservoir structure, according to an embodiment of the present invention; 
         FIG. 7  shows a comparison between the estimated CM to predict the production rate and the measurements, according to an embodiment of the present invention; 
         FIG. 8   a  shows the relative well locations for the injector and the producer wells; 
         FIG. 8   b  shows an estimate probability map for a pilot area, upper side; 
         FIG. 9   a  shows the relative well locations for the injector and producer wells; 
         FIG. 9   b  shows an estimate probability map for a pilot area, lower side; and 
         FIG. 10  depicts the estimated probability map for the whole field with a known fault, according to an embodiment of the present invention. 
     
    
    
     DETAILED DESCRIPTION 
     In accordance with an embodiment of the present invention, injection and production data are used to estimate high permeability channels. This may be extended to the 3D case to identify the shale discontinuities which causes cross layer flow. The present method may be implemented without stoppage in production schedules. 
     In particular, the inventors have determined that it may be difficult to identify shale discontinuities, particular in a layered reservoir with large permeability contrast. In order to overcome this difficulty, the inventors have investigated use of flow rates as an indicator. In particular, because the shale discontinuities provide channels for the fluid to flow between different layers, when an injection rate is changed, and producers in the same layer are shut-in, it should be possible to detect changes in production in different layers. By controlling injection schedule design and including shut-in for producers, well-pair response time in a 3D layered reservoir may be estimated. 
     With the injection/production data, a capacitance model (CM) may be used to characterize the flow property for each well-pair. Specifically, the model can be viewed as a virtual streamline approach to explain the production according to the interaction between well-pairs. In CM, a ‘time delay’ constant is introduced to characterize the time delay effect of the injection signal at the producers. This “time delay” constant is associated with the cross layer travel path, which is determined by the location of shale discontinuities. In this regard, the injector and producer are treated as the transmitter and receiver, respectively, and use the time delay constant as the travel time for each well-pair. One special property of this problem is that shale discontinuities may have very different flow properties. This approach may be referred to as “high contrast travel time tomography.” 
     After determining a time delay constant, the problem of finding the shale discontinuities can be modeled as a cross-hole travel time tomography. Travel time tomography aims to reconstruct an interior velocity model based on measured first-arrival times between transmitters and receivers. The velocity model captures the physical properties of the region where the signal transmission occurs. The transmitted signal can be a seismic wave, an acoustic sound wave, and even a fluid pressure wave. This technique can characterize large scale elastic media, in applications such as seismic geophysical exploration and acoustic tomography, in the atmosphere&#39;s surface or the ocean layers. However, the geophysical problem is different from many other transmission tomographic reconstruction problems (X-Ray CT, Positron Emission Tomography CT, etc.) where the straight-line trajectory assumption is commonly used. Indeed, in the geophysical problem, the acoustic trajectory may bend severely if the local velocity variation is high. Instead of using the straight line trajectory assumption, travel paths can be characterized by Fermat&#39;s principle, where the travel time observed corresponds to the path traversed in the shortest time, where the time can be obtained as a path integral of the slowness (reciprocal of velocity) function integrated along the travel path. 
     In many geophysical applications, it is very common to have heterogeneous structures where the velocity contrast is high, e.g., an area may have twice the velocity than other areas. Example scenarios where this situation can be encountered arise in different applications, e.g., using seismic waves to find permeable fracture zones in ground-water flow, monitoring the water/oil saturation in vegetable oil bio-remediation projects, etc. Different from the case that the velocity distribution only has small variation, the travel path in high velocity contrast media not only bends severely but is almost dominated by these high velocity structures that form high transmissive channels. Because the straight line travel path approximation is not valid, the reconstruction of the velocity model becomes a nonlinear inverse problem. Furthermore, conventional reconstruction methods are based on iterative linearized algorithms to approximate the travel path. Given that in many practical situations, only sparse measured travel time data is available, these methods suffer from problems due to low ray-coverage and severe path bending in the high contrast velocity medium. 
     Water-flooding is used as a technique in enhanced oil recovery (EOR). Water is injected in a controlled manner in order to provide pressure support that can slowly sweep oil into the production wells. During this process, the permeability (measurement of ability to transmit fluid) of open fractures can be orders of magnitude higher than that of surrounding tight rocks, providing fast pathways for the flow. Thus, travel time through a fracture (which could be modeled as a line in 2D or a plane in 3D) is much faster than through surrounding rocks. The flow properties of the reservoir are dominated by these highly ‘transmissive’ fracture structures. If a fracture is close to both an injection well and a production well, most water will directly flow through the fracture and thus fail to displace oil in other areas which is called “water cycling” and thus significantly reduces the efficiency of oil recovery. Therefore, understanding these fractures is useful in flow characterization to enhance oil recovery. 
     In one embodiment, high velocity structures (HVS) in a relatively slow homogeneous background are determined. The velocities are considered discrete which calls for solving a ‘discrete’ travel time tomography problem. For example, a reservoir can be modeled as a layered structure, where each layer having a different hydraulic conductivity. 
     Object-based models can be used to represent the HVS. For example, objects can be pre-defined convex polygons, and the geometrical shape of HVS is represented by a combination of these fundamental objects. This allows an arbitrary shape to be represented with convex objects, and prior information about the HVS can be input into the object-based models. For instance, in the fracture characterization problem, the geometrical shape of fractures can be approximated by lines in 2D models or planes in 3D models. The shape of the convex fundamental object can be defined as a ‘line’. As a result, a small number of objects is sufficient to approximate the shape of HVS, which reduces the dimensionality and uncertainty of the problem. Furthermore, the travel path tracking procedure can be simplified by only considering the shortest time path between objects instead of all grids in the spatial domain. This approach can reduce the number of unknown variables to the object parameters and avoids the ‘lack of path coverage’ problem arising in grid-based models which provides a more stable result in inversion. 
     In order to estimate the velocity model with the measured travel time, a nonlinear inverse problem may be solved. In inverse problems, the estimated velocity model (the solution) may not be unique to the given measurements. One popular approach to handle the non-uniqueness of the solution is to apply a regularization to favor certain properties in the model. The regularization methods can lead to solutions that balance the data-fitting and model-penalty. 
     In another embodiment, an alternative approach is to estimate the probability distribution for the model space according to the data-fitting. This gives a full description of the relative probabilities of the different models so that all likely solutions may be considered. 
     In one embodiment, the Bayesian approach may be selected to estimate the probability density map for the model parameter space. Due to the large dimension of the model parameter space, an accelerated random walk algorithm is provided to explore the model space. In one embodiment, based on the Hamiltonian Monte Carlo method (HMC), an additional friction term can be added in the simulation of dynamic system. The samples are more likely to fall into local minimum and achieve efficient sampling on high probability regions. The HVS property can then be used to approximate the neighboring probability distribution. The results are presented as a probability map showing where the high velocity objects are more likely to appear in the underlying structure of the system. 
     In one embodiment, the HVS properties can be used to simplify the path tracking step, which allows the Hamiltonian Monte Carlo (HMC) sampling process to be more efficient. Furthermore, the monotonicity of the travel time as a function of high velocity object size can be exploited to approximate the probability distribution in the space of model parameters. 
     In one embodiment, the high contrast travel time tomographic reconstruction algorithm described above can be implemented in the geophysical context and applied in petroleum or gas fields. Understanding the heterogeneous structure of oil and gas reservoirs can provide information to sweet spot mapping, drilling, production forecasting and operation efficiency optimization. In enhance oil recovery (EOR) process, fluids such as water, gas, or chemicals are injected to increase the amount of oil that can be extracted from the reservoir. The water-flood process is a process where water is used as the injected fluid. Water injection provides pressure maintenance in the reservoir while displacing part of residual oil to production wells. To maximize operation efficiency of the water-flood, steps may be taken to uniformly “push” the residual oil where the presence of high permeability channels affect the horizontal and vertical sweep efficiency. For example, the water-flood process can be applied to a tight fractured reservoir. In “tight” fractured reservoir, fractures provide conduits for the fluid flow and these conduits provide fast pathways for the injected fluid. Hence, these pathways dominate the sweep efficiency in water-flood. 
     For example, if a high permeability channel is in close proximity to one injection-production well-pair, most injected water may flow through this channel and may fail to sweep the oil in other regions. This phenomenon is called “water cycling”, which decrease the sweep efficiency and increase the water cut (the ratio of water produced compared to the volume of total liquids produced). To optimize the water-flood efficiency, it is desirable to identify the locations of these high permeability channels in the reservoir by appropriate changes in injection rates in different wells. 
     The seismic cross-hole tomography method provides very detailed geophysical structures. However, in this method, it can be difficult to use the seismic velocity to infer the flow property directly. Other methods use tracer tests which measure the diffusion process between the injector-producer pairs and provide the flow permeability distribution in a direct way. A common aspect of these methods is that they all require additional equipments and may interrupt daily operations. Seismic cross-hole testing requires inducing seismic waves and deploying sensors to monitor the reflected waves. In tracer tests, chemicals or radioactive materials are injected and concentration response is monitored to map inter-well flow properties. In addition, repeating the testing in the same area is complicated and one needs to wait a long time for the tracer concentration to return back to background levels. 
     In a water-flood field, the injection/production data is the most abundant data source and there is freedom to control the injection rates. Thus, the injection rate can be changed while monitoring the change in production wells. The estimated response time between well pairs enables a user to infer the structure of a reservoir using travel time tomography. In one embodiment, to estimate the travel time, the water-flood reservoir can be modeled as a linear system by considering water injection rates as inputs and total fluid production rates as outputs. For example, in pulse testing, sensitive differential pressure gauges can be used to monitor the resulting pressure response at adjacent wells. 
     In another embodiment, the reservoir can be treated as a multi-input and multi-output system and estimate all the inter-well response by the production rates. By using a capacitance model (CM), the relation between the injection-production response and the “time delay” constant can be built. The relation between the injection-production response and the time delay is approximately proportional to the pressure wave propagation time. Therefore, the time delay can be estimated from the injection/production data. The obtained time delay can be used as the travel time between well pairs. This procedure has the benefit of being able to be applied without additional cost or significantly affecting daily operation. 
     In tight fractured reservoirs, fractures can be viewed as high permeability channels where the permeability contrast is relatively high, about 10 5  times permeable than the host rock formation. In this case, the travel time for the pressure wave to propagate through the high permeability channels can be considered almost negligible. The travel time data is restricted to the injection-production well-pair locations sparsely located in the field. The reconstructed result is related to the density of travel ray-path. Thus, the reconstructed reservoir image resolution is fundamentally limited by the spatial distribution of well locations. 
     Based on these conditions, the high permeability channels in a reservoir can be characterized by solving a high contrast travel time tomography problem with sparse data. If the grid model is used and the conventional iterative least square methods are applied, the reconstruction results would be poor. To identify these high permeability channels an object-based model which employs “lines” as fundamental objects can be used to represent fractures in 2D. 
     To estimate the travel time from historical data, a physical model is used to describe the injection/production response. The capacitance model can be used to characterize the system response, where the time delay constant in CM is related to the control volume. When the flow passes through a homogeneous region, the control volume can be assumed to be proportional to the travel path. By definition, the travel time is the travel path divided by the local velocity. The time delay constant is proportional to the control volume divided by the productivity index. In a tight fractured reservoir, a simple reservoir model with embedded high permeability channels in a homogeneous background may be used. Thus, if the productivity index in high permeability channels is higher than the background, the time delay constant inside the high permeability channels can be ignored. The time delay constant is proportional to the travel time in the high contrast reservoir model. 
     The capacitance model provides a powerful tool to evaluate the water-flood performance by treating the injection/production as the system input/output and uses signal processing techniques to estimate the unknown system parameters, which relates to reservoir physical properties. The mathematical formulation for one injection-production well pair in CM can be expressed by the following equation (1). 
     
       
         
           
             
               
                 
                   
                     
                       τ 
                        
                       
                           
                       
                        
                       
                         
                            
                           q 
                         
                         
                            
                           t 
                         
                       
                     
                     + 
                     
                       q 
                        
                       
                         ( 
                         t 
                         ) 
                       
                     
                   
                   = 
                   
                     
                       i 
                        
                       
                         ( 
                         t 
                         ) 
                       
                     
                     - 
                     
                       τ 
                       * 
                       J 
                        
                       
                           
                       
                        
                       
                         
                            
                           
                             p 
                             wf 
                           
                         
                         
                            
                           t 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   1 
                   ) 
                 
               
             
           
         
       
     
     Where i(t) is the injection rate and q(t) is the total production rate. P wf  represents the flowing bottom hole pressure (BHP) and J stands for the productivity index. In the case where the reservoir is mature and the bottom-hole pressure is fixed, then it is possible to approximate the j-th well production rates based on only the contribution from injectors by using the following equation (2). 
     
       
         
           
             
               
                 
                   
                     
                       q 
                       j 
                     
                      
                     
                       ( 
                       t 
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       i 
                     
                      
                     
                       
                         ∫ 
                         
                           t 
                           0 
                         
                         t 
                       
                        
                       
                         
                           
                              
                             
                               
                                 - 
                                 
                                   ( 
                                   
                                     t 
                                     - 
                                     ζ 
                                   
                                   ) 
                                 
                               
                               / 
                               
                                 τ 
                                 ij 
                               
                             
                           
                           
                             τ 
                             ij 
                           
                         
                          
                         
                           
                             I 
                             i 
                           
                            
                           
                             ( 
                             ς 
                             ) 
                           
                         
                          
                         
                            
                           ς 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   2 
                   ) 
                 
               
             
           
         
       
     
     The discrete form of equation (2) can be written as equation (3). 
     
       
         
           
             
               
                 
                   
                     q 
                      
                     
                       ( 
                       t 
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       k 
                     
                      
                     
                       ( 
                       
                         
                           1 
                           τ 
                         
                          
                         
                            
                           
                             - 
                             
                               ( 
                               
                                 k 
                                 / 
                                 τ 
                               
                               ) 
                             
                           
                         
                          
                         
                           i 
                            
                           
                             [ 
                             
                               t 
                               - 
                               k 
                             
                             ] 
                           
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     Where t is equal to 1, 2, . . . , N. For multiple injection-production well pairs, additional parameters l ij  for inter-well connectivity, which represent flow distribution from injector i to different producers j. The j-th well production rate is the sum of flow from different injector wells i. This can be expressed by the following equation (4). 
     
       
         
           
             
               
                 
                   
                     
                       q 
                       j 
                     
                      
                     
                       ( 
                       t 
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       i 
                     
                      
                     
                       
                         ∑ 
                         k 
                       
                        
                       
                         
                           l 
                           ij 
                         
                          
                         
                           ( 
                           
                             
                               1 
                               τ 
                             
                              
                             
                                
                               
                                 - 
                                 
                                   ( 
                                   
                                     k 
                                     / 
                                     τ 
                                   
                                   ) 
                                 
                               
                             
                              
                             
                               i 
                                
                               
                                 [ 
                                 
                                   t 
                                   - 
                                   k 
                                 
                                 ] 
                               
                             
                           
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
           
         
       
     
     By using equation (4), a linear multi-input multi-output (MIMO) system can be used to model the relation between injection/production rate. The impulse response is controlled by the “time delay constant” τ ij  between each injector-producer pair which is defined by the total compressibility c t , the productivity index J, and the pore volume Vp which are associated with the control volume between injector i and producer j pair. The time delay constant τ ij  between each injector—producer pair can be expressed by the following equation (5). 
     
       
         
           
             
               
                 
                   
                     τ 
                     ij 
                   
                   = 
                   
                     
                       ( 
                       
                         
                           
                             c 
                             t 
                           
                            
                           
                             V 
                             p 
                           
                         
                         J 
                       
                       ) 
                     
                     ij 
                   
                 
               
               
                 
                   ( 
                   5 
                   ) 
                 
               
             
           
         
       
     
     From above, it is clear that if there is an area with higher productivity index J, the time delay constant will be much smaller. This implies that if the time delay constant for a specific well-pair is very small, it is probable that a majority of the flow path is through a high permeability layer. 
     For a homogeneous reservoir, the control volume is roughly proportional to the distance between the injector and producer. When a cross-layer flow in a layered reservoir is considered, the travel path will be the path with the least time cost in accordance with Fermat&#39;s principle. This implies that the flow will take advantage of the higher permeability layer. 
     A total control volume can be modeled as the cascade of several small control volumes, e.g., for example two control volumes which can be two different layers. Each control volume is proportional to a line segment which represents the flow path through the control volume. Therefore, the time-delay constant between a well pair will be the summation of two parts: the time delay constant of the first layer (e.g., upper layer) and the time delay constant of the second layer (e.g., lower layer). Each is given by the ratio of flow path in different layers, which is determined by the location of the shale discontinuities because they provide the entries for the fluid to travel through different layers. 
     In order to get a reliable estimate of the parameters for CM, a Pseudo-Noise (PN) sequence may be used as the injection schedule. By using PN sequence as the input, it may be possible to achieve the lowest covariance for the parameter estimation. In addition, by using PN sequences as inputs, lowest average cross-correlation between inputs may also be achieved. 
     When the injection rate in the lower permeability layer is changed, usually it is possible to detect the change of production even for the producers located in the higher permeability layer. However, in some instances, it may be difficult to detect this cross-layer flow if the injection rate in the high permeability layers is changed. This is because the high permeability layer provides a fast pathway for the fluid so only a small amount will flow through the shale discontinuity to the low permeability layer. Therefore, when designing the injection schedule it may be useful to shut-in the producers in the same layer to force the fluid to flow through different layers. In this situation, the cross layer flow can be detected by monitoring the change of production. 
     The algorithm to identify the shale discontinuities with cross-layer time delay constant can be referred as the “high contrast travel time tomographic reconstruction”. One of the issues with this approach is that the sensor locations may be limited in the injectors/producers, and the flow velocity contrast of the shale discontinuities may be very high. If the conventional grid based inversion method is used, the result may result in a blurred image and it may be difficult to identify the shale discontinuities. 
     As such, instead of using the grid based method, an object-based representation of the shale discontinuities may be chosen. The reconstruction algorithm will modify the object parameters for the shale discontinuities to match the measured delay time constant. It can be described as a ‘Forward-Backward’ iterative procedure. 
       FIG. 1  is a flow chart illustrating steps of a method for determining the time delay constant, according to an embodiment of the present invention. In the forward step, a time delay is calculated based on the current or initial locations of shale discontinuities based on a current or initial reservoir model, at S 10 . At S 12 , a test is performed to determine whether or not a difference between the time delay and a measured time delay is small enough, i.e., smaller than a set threshold (e.g. smaller than 1% to 5% of the measured value). If the calculated result fits the measurement (i.e., the difference between the measurement and the current estimate is small or smaller than the set threshold), then the selected current reservoir model can be considered as a solution and the algorithm is stopped, at S 14 . If not, the ‘backward’ step is performed to update the locations of shale discontinuities (i.e., the initial reservoir model is updated) based on the difference of the current predicted response time with the measurements, at S 16 . Therefore, if the difference is greater than the set threshold, the method includes updating the initial reservoir model based on the difference, calculating the time delay based on the updated reservoir model; and determining whether or not a difference between the time delay based on the updated reservoir model and the measured time delay is smaller than the threshold. This loop procedure can be repeated until convergence between the calculated time delay and the measured time delay is achieved. That is, the loop procedure is repeated until the difference between the calculated time delay and the measured time delay is smaller than the set threshold. 
     The detail of the “forward” step can be described as follows. With the current location of shale discontinuities, the flow travel path can be calculated by Fermat&#39;s principle. The shale discontinuities act as entry points to the different layers, and the time delay constant can be approximated by the cascade of different layers which is proportional to the travel distance inside. In this case, the time delay constant is well defined by the locations of shale discontinuities. 
       FIG. 2  is a schematic representation of an injector/producer relationship in a multilayer reservoir having a shale discontinuity, according to an embodiment of the present invention. Injector well  10  and producer well  12  provided within rock formation  13  are separated by distance L. As illustrated in  FIG. 2 , there two shale layers  14 A and  14 B. Shale layer  14 A has a lower permeability and shale layer  14 B has a higher permeability. As shown in  FIG. 2 , there is a shale discontinuity  16  between shale layer  14 A and shale layer  14 B, the discontinuity corresponding to a fracture in the rock formation  13  where fluid flow (indicated schematically by the arrows) is privileged. For example, if it is assumed that the flow velocity of the layer (shale layer  14 A) having lower permeability is equal to 1 (VL=1) and the flow velocity of the layer (shale layer  14 B) having a higher permeability is 10 (VH=10), the flow path can be separated as two segments belonging to the two different layers  14 A and  14 B and the delay time can be calculated using the following equation (6). 
       Travel time= L   1 /1 +L   2 /10, with  L   1   +L   2   =L   (6)
 
     where L represents the geometrical distance between the injector-producer well  10  and  12 , and the ratio of two segments is controlled by the location of shale discontinuities. The distance L 1  corresponds to the distance traversed by the injected fluid through layer  14 A. The distance L 2  corresponds to the distance traversed by the injected fluid through layer  14 B. The travel time within the discontinuity or fracture  16  is considered small relative to the travel times of fluid flow within the layers  14 A and  14 B. Therefore, the travel time within fracture  16  is not taken into account when calculating the total travel time using equation (6). 
       FIG. 3  is a schematic illustration of a flow in an injector well and producer well relationship illustrating Fermat&#39;s principle in determining the location of shale discontinuities, according to an embodiment of the present invention. Injector well  20  and producer well  22  provided within rock formation  23  are separated by distance L. As illustrated in  FIG. 3 , there are two types of shales  24 A and  24 B. Shale  24 A has a lower permeability and shale  24 B has a higher permeability. As shown in  FIG. 3 , there is a shale discontinuity  26  between shale  24 A and shale  24 B, the discontinuity corresponding to a fracture in the rock formation  23 . In this model of the earth, the relationship between the various shales (i.e., shales with various permeabilities) and the position of the discontinuity is not as straight forward as in the example shown in  FIG. 2 . Therefore, in order to determine an appropriate model that provides an estimate of the travel time and thus the location of the discontinuity, the backward step in the loop S 16  is performed. In the ‘backward’ step, the locations for the shale discontinuities  26  is updated base on the mismatch of the measured and predicted time delay constant τ. Time delay constant τ is the sum of time delay τ(1) within shale  24 A and time delay τ(2) within shale  24 B. Time delay τ(1) is proportional to the traversed length L 1  within shale  24 A and inversely proportional to the flow velocity within shale  24 A. Similarly, time delay τ(2) is proportional to the traversed length L 2  within shale  24 B and inversely proportional to the flow velocity within shale  24 B. 
     Because the shale discontinuities are represented by objects with parameters, we can formulate the updating as a parameter estimation problem. Current object parameter can be expressed by the following equation (7). 
     
       
         
           
             
               
                 
                   θ 
                   = 
                   
                     
                       Argmin 
                       θ 
                     
                      
                     
                        
                       
                         
                           T 
                            
                           
                             ( 
                             θ 
                             ) 
                           
                         
                         - 
                         t 
                       
                        
                     
                   
                 
               
               
                 
                   ( 
                   7 
                   ) 
                 
               
             
           
         
       
     
     where T(θ) is the calculated response time based on current object parameters θ and t is the measurement. A gradient search is used to find the optimal parameters to represent the location of the shale discontinuities. 
       FIGS. 4   a - 4   c  illustrate a testing scenario in a simulated application of a method in accordance with an embodiment of the invention. A two layer reservoir model with high permeability contrast was built and a commercial simulator (CMG) is used to test the method. The testing case was a line drive with 5 injectors and 5 producers. Three of the injectors/producers are in the low permeability layers, and the two other injectors/producers are in the high permeability layers. A simple shale discontinuity was placed and high permeability assigned to it to simulate the effect. For the reservoir model configuration, a high permeability channel is simulated. In the simulation, the high permeability channel connects injector well  5  with producer well  3 . A pulse testing is applied in the simulation. The pulse testing uses a step or square function or wave as the change of injection rate. The injection rate is changed each time so the response can be easily estimated from the change of production rate. The travel time is defined by the time to reach 90% of the final change of the production rate.  FIG. 4   a  shows the configuration of injector wells I 2  and I 4  relative to production wells P 2  and P 4  in the first layer.  FIG. 4   b  shows the configuration of injector wells I 1 , I 3  and I 5  relative to producer wells P 1 , P 3  and P 5  in the second layer.  FIG. 4   c  illustrates the shale layers and their respective permeabilities (shown on the scale strip) within the rock formation and the relative position of the 5 injector wells I 1  through I 5  and  5  producer wells P 1  through P 5 . 
     A PN sequence was used as the injection schedule and the producers in the same layer were shut-in for the testing period. The changes of production can be measured for the producers located in different layers. A time delay constant was retrieved in CM by injection/production data, and used to identify the location of shale discontinuities as shown in  FIGS. 5   a - 5   c.    
       FIG. 5   a  depicts two plots with two different injection rates.  FIG. 5   b  depicts two plots with two production rates corresponding to the two injection rates shown in  FIG. 5   a .  FIG. 5   c  represents a plot of an estimated FIR and a CM to determine a time delay constant for the simulated method. The plot represents a correlation between injection and production as a function of time (the horizontal or x-axis representing time and the vertical or y-axis representing the normalized response). As it can be noted, the average time delay to obtain a change of production at the production site is in this case about 5 to 6 days. 
     To apply the method in a real field, one of the useful parameters is the sampling frequency. Another useful parameter is accuracy of the data. Usually it is possible to obtain reliable and frequent injection data, but production rates are obtained from the well-test, which can performed periodically, e.g., on a weekly basis. This factor tends to limit the time-resolution when the time delay constant is measured. For example, assuming the time delay constants for the well-pairs are 1 and 4 days, when injection rates are increased, it is possible to measure almost the same changes in production rate if it is measured by well-test. To deal with this problem, it is possible to use an inferred production value provided by a controller pump which provides more frequent production data. 
       FIG. 6   a  illustrates an assumed initial structure of the high permeability channel or discontinuity, according to an embodiment of the present invention. The assumed initial structure is used in the initial model to estimate an initial time delay.  FIG. 6   b  illustrates an output of the algorithm, according to an embodiment of the present invention.  FIG. 6   b  depicts the position of the high permeability channel or shale discontinuity obtained through refining the initial model until the difference between the estimate time delay and the measured time delay is relatively small and below a certain threshold.  FIG. 6   c  illustrates the simulated reservoir structure, according to an embodiment of the present invention. The scale strip next to the graph provides a scale of the permeability of the rock formation. The bar within the graph indicates the position of the discontinuity. A good match is obtained between the measurement and the calculated position or location of the discontinuity (fault) as well as the direction of the discontinuity (fault). 
     A field experiment under water-flood in an oil field is performed. The field is considered to be a tight fractured reservoir. In the experiment, a pilot area with 12 injection wells is considered. The well installation is a line-drive and roughly parallel to the direction of fractures which are estimated from seismic survey. 
     Two types of systems in the oil field are used to provide the well production data. One is the “well-test” data, which accumulates a three phase production rate of one specific well for a short period of time. However, many wells share the same separation facility. Therefore, these wells are tested sequentially. This limits the sampling rate of the production rates, because the sampling period is the sum of the testing period for all production wells. In this case, one measurement point every two weeks is taken which makes our production rate sparsely sampled in the well-test data. 
     The other data source is the pump-of-control (POC) data, which measures the pump load change in every stroke. If the number of strokes counted is multiplied by the load change, the theoretical total daily production can be calculated. To get a reliable data with high sampling rate, the POC data is used and the pumping efficiency is calibrated with the well-test data. After calibrating with the pumping efficiency, the POC data would be able to provide a daily production rate. However, we do observe unusual low production rate, which might be outliers due to the power outage or loss of transmission. Thus, in the linear model estimation we choose l 1  instead using l 2  norm as the error metric, which is more balance between fitting good data points and outliers. 
     To estimate the CM parameters, the direct approach is chosen in this case. The reason is because the input data length is limited (about 90), and thus, the multi-stage approach is not suitable because the unknown parameters in FIR model will be much larger than the measured data. From a previous discussion in the above paragraphs, it is established that the direct estimation uses fewer unknown variables but needs to solve a nonlinear optimization. Therefore, the estimated CM might be a solution in a local minimum and not represent the system well. To verify the estimated CM, cross-validation is used which divides the measured data into two parts. The first part is the training data, which is used to estimate the CM parameters. The second part is the testing data, where estimated CM is used to predict the production rate and compare with the measurements. If the error is too high, the estimated CM is determined to be not accurate and the non-linear optimization is re-run.  FIG. 7  shows a comparison between the estimated CM to predict the production rate and the measurements, according to an embodiment of the present invention. 
     Following the same principle in numerical simulation, the travel time is only used when the inter-well connectivity is greater than 10%. The pilot area is separated into upper and lower side and the travel time between well pairs inside each region is estimated. The reason why the estimation for the well pair is not performed for the whole field is that that solving for the whole field needs a greater number of variables. In the line drive water-flood, it can be assumed the flow only comes from the nearby two rows of injectors. Hence, solving two smaller fields is selected instead. 
       FIG. 8   a  shows the relative well locations for the injector and the producer wells.  FIG. 8   b  shows an estimate probability map for pilot area  1 , upper side. It can be noted that the estimated high permeability is roughly parallel to the installation of the wells which is consistent with the prior information of fracture direction.  FIG. 9   a  shows the relative well locations for the injector and producer wells.  FIG. 9   b  shows an estimate probability map for pilot area  1 , lower side. Similar to  FIG. 8   b , it can also be noted that the estimated high permeability is roughly parallel to the installation of the wells which is consistent with the prior information of fracture direction. 
     The estimate travel time results for the upper part of pilot area is shown in  FIG. 8   b , and the lower part is shown in  FIG. 9   b . In the upper part, two reliable production rate measurements are obtained, and the travel time used are t_(1,1), t_(3, 2), t_(4, 2), t_(5,1), t (5, 2), t_(6,1), t_(7,1), t_(7,2), t_(8,1), t_(8,2), a total of 8 well pairs. The result shows in the upper side, the estimated high permeability channels are close to production well  2  and roughly parallel to the installation of wells. This is consistent with the fact that the well production is a 24 hours run, high production well, which implies that there are high permeability channels nearby. In the lower part, one production well is present and the estimated travel time are t_(7,3), t_(9,3), t_(10,3), t_(11,3), and t_(12,3) where there are fewer measured travel times to estimate the high permeability map. The result in two areas are then combined and compared with the known fault location. The result shows the appearance probability of high permeability channels in the fault zones is very low, and roughly parallel to the well installation which agrees with the seismic survey.  FIG. 10  depicts the estimated probability map for the whole field with a known fault, according to an embodiment of the present invention. The position of the known fault is indicated on the probability map. The results show the permeability channel has a very low appearance probability in the fault zone consistent with the prior survey. 
     The “time delay constant” can be used in a CM method to detect high permeability channels. The method uses the injection-production. The method uses the change of injection as the active probing and detect the field changes in real-time without altering the average daily production. In order to apply this to real field data, some practical issues are noted. First, the data sampling period and quality can play an important role. Usually, a reliable daily injection rate data is obtained. However, production rates are often obtained from bi-weekly well-test data. This may reduce the time resolution of estimated travel times. The POC data calibrated with well-test data can be used to get the daily production rate. In addition, in practical applications, the distribution of well locations is related to spatial resolution. For example, in an extreme case where one injector-producer pair is in the horizontal direction, any high permeability channel exactly in the vertical direction will not affect the lag time. Therefore, the high permeability channels may not be visible in this situation. In order to detect the high permeability channel in any arbitrary direction, the wells can be uniformly located in the field and to cover all angles. 
     As will be appreciated, the method as described herein may be performed using a computing system having machine executable instructions stored on a tangible medium. The instructions are executable to perform each portion of the method, either autonomously, or with the assistance of input from an operator. In an embodiment, the system includes structures for allowing input and output of data, and a display that is configured and arranged to display the intermediate and/or final products of the process steps. A method in accordance with an embodiment may include an automated selection of a location for exploitation and/or exploratory drilling for hydrocarbon resources. Where the term processor is used, it should be understood to be applicable to multi-processor systems and/or distributed computing systems. 
     Those skilled in the art will appreciate that the disclosed embodiments described herein are by way of example only, and that numerous variations will exist. The invention is limited only by the claims, which encompass the embodiments described herein as well as variants apparent to those skilled in the art. In addition, it should be appreciated that structural features or method steps shown or described in any one embodiment herein can be used in other embodiments as well.