# Patent

## Title
Flow of self-diverting acids in carbonate reservoirs

## Abstract

Two new flow parameters derived from laboratory core-flood experiments are used in building mathematical models to predict the performance of an acid treatment when treatment is made with self diverting fracturing acids. The two new variables are:

    
    
        ΔPr is defined as the value of Δp (in the core flood experiment) when Δp switches from a first to a second linear trend at time tr 
        ⊖r is the number of pore volumes injected when the switch occurs.

## Background
FIELD OF THE INVENTION
The invention relates to acid stimulation of hydrocarbon bearing subsurface formations and reservoirs. In particular, the invention relates to methods of optimizing field treatment of the formations.
BACKGROUND
Matrix acidizing is a process used to increase the production rate of wells in hydrocarbon reservoirs. It includes the step of pumping an acid into an oil- or gas-producing well to increase the permeability of the formation through which hydrocarbon is produced and to remove some of the formation damage caused by the drilling and completion fluids and drill bits during the drilling and completion process.
In order to predict the outcome in the field of the pumping of an acid, or of acid stages, into a reservoir, engineers go through a design process, which can be divided into several steps. In the first step, for example, core flood experiments are carried out, where different acids are injected, for testing, into cylindrical rock cores under various conditions. During such tests, many parameters can be varied, such as an injection rate Q, a temperature T, an acid formula Ac, and a rock type Ro.
In the core flood experiment, as acid flows into the rock, it dissolves part of the rock matrix and increases the overall permeability of the core with time. Depending on the combination of the above parameters, the dissolution pattern inside the rock can vary between face dissolution (also known as compact dissolution), wormholing dissolution and uniform dissolution. Face dissolution corresponds to the regime where acid flows so slowly that it dissolves the rock through the rock face only, located at the interface between the acid and the core. This interface moves slowly in the flow direction as more and more rock gets dissolved with time. Wormholing dissolution happens when acid flows faster than in the face dissolution regime and not all the acid is spend at the rock face. Live acid enters the core and, due to instable dissolution fronts, fingers of live acids propagate into the rock forming structures known as wormholes. If acid is pumped fast enough for the amount of acid spent during the residence time of the fluid into the core is very small, then, the acid concentration is constant within the rock and the matrix is dissolve in a uniform way. These three known dissolution regimes give rise to different acid efficiencies. Acid efficiency is measured as the amount of acid that is required by the rock core to increase its permeability to a pre-set value kw, for instance 100 times larger than the initial permeability k0 of the sample. The smaller this volume of acid is, the higher the efficiency is. The moment at which this target value of permeability increase is reached is called the breakthrough time, t0. The corresponding volume of acid is called the breakthrough volume, V0.
The measure of pore volumes to breakthrough, denoted ⊖0, (i.e. the breakthrough volume divided by the pore volume of the core PV, where PV is the volume of fluid that can be contained in the core, within the pore network), and its use to predict acid performance during a treatment job has been known to the industry for a long time. For example, pore volume to breakthrough has widely been used as a measure of the velocity at which wormholes propagate into the formation, under various conditions such as mean flow-rate Q, temperature T, rock-type Ro, and acid formulation Ac.
In order to measure pore-volume to breakthrough, acid is pumped at a constant rate Q and the pressure drop Δp across the core is monitored. The initial pressure drop when the acid reaches the inlet core face is called Δp0. When non-self diverting acids such as hydrochloric and acetic are used, as acid flows into the core, the pressure drop declines, mostly linearly. When Δp is virtually equal to 0 (i.e., the core permeability has reached a value kw orders of magnitude larger than the initial permeability k0) the pore-volume injected is recorded as the pore-volume to breakthrough ⊖0.
Recently, acid systems have been developed with the goal of achieving maximum zonal coverage in heterogeneous reservoirs. Such fluids are designed to self-divert into lower permeability zones of the reservoir after having penetrated and stimulated higher-permeability zones. When such systems are pumped using the same procedure as the one described above, the pressure drop Δp across the core may evolve in a very different way as for non-self diverting acids: the pressure does not decline linearly with time and might increase significantly over a certain period of time.
SUMMARY
In various aspects, the methods of the invention are related to the discovery of two new key flow parameters that can be derived from laboratory core-flood experiments, and to their use in building mathematical models to predict the performance of an acid treatment when treatment is made with self diverting fracturing acids. In one embodiment, predictions of the performance of acid treatments based on the models are used to enhance or optimize such treatment.
One important difference in self diverting acid treatment is that the pressure drop Δp across the core observed during the core-flood experiment either increases and then decreases with time or decreases with time at two different rates. In particular, it is observed that Δp has a piece-wise linear evolution. First, Δp evolves according to a first linear relationship with time (or equivalently with volume or pore volume injected). Then, at a certain time tr, it switches to a second linear behavior. Associated with this behavior, two new variables are provided:

    
    
        ΔPr is defined as the value of Δp (in the core flood experiment) when Δp switches from the first to the second linear trend at time tr 
        ⊖r is the number of pore volumes injected when the switch occurs.
    
    


In various embodiments, the two variables are utilized and exploited in methods of predicting the performance of self-diverting acids. Where necessary, mathematical models and algorithms are developed.
When we use the term “acid” here we include other formation-dissolving treatment fluid components, such as certain chelating agents. Further areas of applicability will become apparent from the description provided herein. It should be understood that the description and specific examples are intended for purposes of illustration only and are not intended to limit the scope of the present disclosure.



BRIEF DESCRIPTION OF THE DRAWINGS
FIG. 1 shows a typical experimental apparatus for acid injection into a rock core.
FIG. 2 illustrates pressure-drop for non-diverting acid systems such as HCl. Left: schematic, Right: actual data.
FIG. 3 illustrates pressure drop for self-diverting acid systems such as VDATM. Left: schematic, Right: example of actual data.
FIG. 4 shows a multi pressure tap/transducer core-flooding apparatus.
FIG. 5 shows the evolution of the effective viscosity μe with the number of pore volumes injected for a self-diverting acid.
FIG. 6 illustrates a flow pattern in the core when a self diverting acid is pumped.
FIG. 7 illustrates axisymmetric flow around a wellbore.
FIG. 8 shows an experimental setup for radial flow.
FIG. 9 gives a comparison between method and experiment for radial flow.
FIG. 10 shows a treatment design methodology in the field.
FIG. 11 is a diagram of a reservoir description and wellbore trajectory. The wellbore [32] enters the reservoir [34] at the reservoir top [48], and passes through multiple layers in the reservoir.
FIG. 12 shows HCl treatment results. The wellbore trajectory [32] is shown, along with stimulated regions [50] and virgin un-treated matrix [52].
FIG. 13: VDA treatment results. The wellbore trajectory is shown, along with stimulated regions and virgin un-treated matrix.
FIG. 14 shows the wellbore [32], a wormholed region [54], and a low-mobility region [56], in an optimized VDA treatment. A wormhole penetration profile [58] is shown on the left side of the figure and a low fluid mobility front penetration profile [60] is shown on the right side of the figure.



DESCRIPTION
In one embodiment, the invention provides a method for optimizing the flow rate of a self diverting acid into an acid soluble rock formation during an acid fracturing process. The method comprises

    
    
        predicting treatment performance in the self-diverting acid system on the basis of two flow parameters, the parameters derived from core flood experiments with the self-diverting acid, wherein a fluid is injected into a core and a pressure drop Δp is measured against time at a constant flow rate,

wherein the flow parameters are

        ΔPr, the pressure where a plot of Δp vs. time switches from a first linear trend to a second linear trend; and
        tr is the time at which the switch occurs.
    
    


In another embodiment, the invention provides a method of modeling the pressure in a wellbore during acid treatment with a self diverting acid delivered at a velocity Q, the pressure being determined at a depth z, a distance r from the center of the well, and a time t, the method involving use of functions derived from core flooding experiments wherein a self diverting acid is injected into a core and the pressure along the core is measured as a function of time, the modeling method comprising:
calculating at least one of an effective viscosity μr, a mobility Mr, and a permeability kr, wherein



  
    
      μ
      r
    
    =
    
      
        μ
        d
      
      ⁢
      
        
          Δ
          ⁢
          
            
          
          ⁢
          
            p
            r
          
        
        
          Δ
          ⁢
          
            
          
          ⁢
          
            p
            o
          
        
      
      ⁢
      
        
          Θ
          0
        
        
          
            Θ
            o
          
          -
          
            Θ
            r
          
        
      
    
  




  
    
      M
      r
    
    =
    
      
        k
        0
      
      
        
          
            Δ
            ⁢
            
              
            
            ⁢
            
              p
              r
            
          
          
            Δ
            ⁢
            
              
            
            ⁢
            
              p
              0
            
          
        
        ⁢
        
          
            Θ
            0
          
          
            
              Θ
              0
            
            -
            
              Θ
              r
            
          
        
      
    
  




  and




  
    
      k
      r
    
    =
    
      
        k
        0
      
      
        
          
            μ
            d
          
          μ
        
        ⁢
        
          
            Δ
            ⁢
            
              
            
            ⁢
            
              p
              r
            
          
          
            Δ
            ⁢
            
              
            
            ⁢
            
              p
              0
            
          
        
        ⁢
        
          
            Θ
            0
          
          
            
              Θ
              0
            
            -
            
              Θ
              r
            
          
        
      
    
  



wherein:

    
    
        k0 is the initial absolute permeability of the core, before acid is injected;
        μd is the viscosity of the displaced fluid originally saturating the core before acid is injected;
        μ is the viscosity of the acid;
        Δpr is the pressure drop derived from the core flooding experiments and is the pressure drop at the time tr that the pressure drop changes from a first linear trend to a second linear trend;
        ⊖r is the number of pore volumes delivered to the core at the time tr;
        Δpo is the pressure drop at t=o of the core flood experiment; and
        ⊖o is the pore volume to breakthrough measured in the core flood experiment; and

calculating pressures within the formation on the basis of the effective viscosity μr, the mobility Mr, and/or the permeability kr.

    
    


In another embodiment, the invention provides a method of optimizing acid treatment of a hydrocarbon containing carbonate reservoir with a self-diverting acid. The method involves:

    
    
        carrying out linear core flood experiments varying one or more parameters selected form the group consisting of acid formulation, rock type, flow rate, and temperature;
        deriving the following functions from the experiments, as a function of the parameters:
        
            ⊖o—the pore volume to wormhole/dissolution front breakthrough;
            ⊖r—the pore volume to resistance zone breakthrough; and
            Δpr—the pressure drop at resistance zone breakthrough;
        
        
        writing equations of a flow model based on the functions;
        solving the equations in an arbitrary flow field in a simulator;
        using the simulator in an optimization loop together with known and/or estimated reservoir parameters; and
        calculating at least one of the following from the simulator optimization loop:
        
            stage and rate volumes of the acid treatment;
            fluid selection for the acid treatment;
            wormhole invasion profile; and
            skin profile.
        
        
    
    


In order to predict the outcome of the pumping of an acid, or of acid stages, into a reservoir, engineers go through a design process, which can be divided into several steps. In the first step, different acids are injected, for testing, into cylindrical rock cores, under various conditions. FIG. 1 is an illustration of a typical experimental setup used for injecting acid into a core. A pump [2] pumps a fluid, for example an acid, through an accumulator [4] into a core [6] held in a core holder [8]. During such tests, the following parameters will normally be varied:
Injection rate: Q
Temperature: T
Acid formulation: Ac
Rock type: Ro
As acid flows into the rock, it dissolves part of the rock matrix and increases the overall permeability of the core with time. Depending on the combination of the above parameters, the dissolution pattern inside the rock can vary between face dissolution (also known as compact dissolution), wormholing dissolution, and uniform dissolution. These three dissolution regimes give rise to different acid efficiencies. Acid efficiency is measured as the amount of acid that is required by the rock core to increase its permeability to a pre-set value kw, for instance 100 times larger than the initial permeability k0 of the sample. The smaller this volume of acid is, the higher the efficiency is. The moment at which this target value of permeability increase is reached is called the breakthrough time, t0. The corresponding volume of acid is called the breakthrough volume, Vol0.
The measure of pore volumes to breakthrough, denoted ⊖0, (i.e. the breakthrough volume divided by the pore volumes of the core PV (the volume of fluid that can be contained in the core), and its use to predict acid performance during a treatment job has been known to the industry for a long time. If we define Vol as being the geometrical volume of the core and φ0 the initial porosity of the core (i.e. the fraction of the core volume that can be occupied by a fluid through the pore space network), these parameters are linked to each other as follows:




  
    
      
        
          Θ
          0
        
        =
        
          
            
              Vol
              0
            
            PV
          
          =
          
            
              
                
                  Qt
                  0
                
                PV
              
              ⁢
              
                
              
              ⁢
              where
              ⁢
              
                
              
              ⁢
              PV
            
            =
            
              
                φ
                0
              
              ×
              Vol
            
          
        
      
    
    
      
        (
        1
        )
      
    
  




Pore volume to breakthrough has widely been used as a measure of the velocity at which wormholes propagate into the formation, under various conditions such as mean flow-rate Q, temperature T, rock-type Ro, and acid formulation Ac.
Typically, multiple pressure taps are installed down the length of the core holder; FIG. 1 shows an inlet pressure tap [10], that has an inlet pressure pi, and a second pressure tap [12], that has a pressure away from the inlet pL, at a distance [14], denoted L, from the inlet. The cross sectional area of the core, A, for example at the core face, is shown at [16]. In order to measure pore-volume to breakthrough for a non-self diverting acid, acid is pumped at a constant rate Q and the pressure drop Δp across the core is monitored. The initial pressure drop when the acid reaches the inlet core face is called Δp0. Then, as acid flows into the core, the pressure drop declines mostly linearly as illustrated in FIG. 2A, in which the breakthrough time, to, is shown at [18], and in FIG. 2B in which the pore-volume to breakthrough, ⊖0, is shown at [20]. When Δp is virtually equal to 0 (i.e., the core permeability has reached a value kw orders of magnitude larger than the initial permeability k0) the pore-volume injected is recorded as the pore-volume to breakthrough ⊖0.
More recently, new acid systems, also known as self-diverting acids such as Viscoelastic Diverting Acid (VDATM), have been used to improve the performance of more classical acid systems such as HCI or organic acids. When such systems are pumped using the same procedure as the one described above, very different Δp behavior can be observed, as is illustrated in FIG. 3. FIG. 3a illustrates the development of Δp with time of pumping (or equivalently, with volume pumped) at a constant rate for two arbitrary systems designated A and B. Results with one self-diverting acid 1, in rock R1, at temperature T1, and rate Q1, are shown by the solid line; results with another self-diverting acid 2, in rock R2, at temperature T2, and rate Q2, are shown by the dotted line.
One important difference is that Δp may increase and then decrease with time or decrease in two regimes at different rates. In particular, it is observed that Δp has a piece-wise linear evolution. First, Δp evolves according to a first linear relationship with time (or equivalently with volume or pore volume injected) in the regions marked as A1 and A2 for two illustrative fluids. Then, at a certain time tr, (or volume Volr) it switches to a second linear behavior, as depicted by B1 and B2 in FIG. 3a. Associated with this behavior, we define two new parameters ΔPr (see FIG. 3a) and the number of pore-volumes to reach Δpr, denoted ⊖r. Δpr is defined as the value of Δp when Δp switches from the first to the second linear trend at time tr. The parameter ⊖r is given by:




  
    
      
        
          Θ
          r
        
        =
        
          
            
              Vol
              r
            
            PV
          
          =
          
            
              Qt
              r
            
            PV
          
        
      
    
    
      
        (
        2
        )
      
    
  




where PV is the pore volume of the core, measured by known methods to determine the volume of liquid held in the core at saturation.

These two parameters constitute a means of predicting the performance of self-diverting acids when used in mathematical models and algorithms as will be explained below. Real data are shown in FIG. 3b. 
Additional experiments have shown that the pressure drop evolution described in FIG. 3, and obtained for self-diverting acid, is due to the existence of a region of low fluid mobility propagating ahead of the wormholes, or ahead of the dissolutions fronts in general. For illustration, a setup as in FIG. 1 is fitted with multiple pressure taps and transducers to measure the pressure along the core during the acid injection experiments, local pressure drops Δpe along the core can be measured. Such a new experimental setup is represented in FIG. 4, in which the inlet pressure tap and transducer is shown at [22] and additional pressure taps and transducers at distances down the core holder are shown at [24].
For a given pair of successive transducers (taps), Le is the distance between the two taps, ke is the permeability of the core and μe is the fluid viscosity between the two taps. According to Darcy's law regarding fluid flow, the measured parameters are interrelated:




  
    
      
        Q
        =
        
          
            
              Ak
              e
            
            
              μ
              e
            
          
          ⁢
          
            
              Δ
              ⁢
              
                
              
              ⁢
              
                p
                e
              
            
            
              L
              e
            
          
        
      
    
    
      
        (
        3
        )
      
    
  




where A is the cross sectional area of the core and Q is the rate of fluid flow. The fluid mobility Me is defined as:





  
    
      
        
          M
          e
        
        =
        
          
            k
            e
          
          
            μ
            e
          
        
      
    
    
      
        (
        4
        )
      
    
  




With the apparatus in FIG. 4, one can:

    
    
        measure Δpe for every pair of transducers, against time,
        and use equations (3) and (4) to determine the fluid mobility Me between every pair of transducers, against time
    
    


From the knowledge of Me at any time, either an effective viscosity or an effective permeability can also be determined:

    
    
        assuming the core permeability k0 is unchanged, equation (4) gives
    
    






  
    
      
        
          μ
          e
        
        =
        
          
            k
            0
          
          
            M
            e
          
        
      
    
    
      
        (
        5
        )
      
    
  




    
    
        assuming the acid viscosity μe is known, equation (4) gives:

ke=μMe  (6)

    
    


The effective viscosity μe of the fluid flowing between pairs of transducers can be monitored against time, or equivalently, against the number of pore volumes injected. The results of one example of such monitoring are illustrated in FIG. 5. The five curves labeled 1, 2, 3, 4, and 5 in FIG. 5 are the values of μe calculated from equations (3), (4), and (5) at the five locations Le in FIG. 4.
Line number 1 (see FIG. 5) corresponds to the zone between the core inlet and the first pressure tap on the core. Line number 2 corresponds to the zone between the first and second pressure taps on the core. The other lines represent the remaining successive pairs in order.
From FIG. 5, it can be seen that, as the self-diverting acid flows into the core, a first zone of finite effective viscosity μe propagates along the core (observed from the viscosity peaks) followed be a zone of virtually zero effective viscosity, or equivalently (using equation (4)), a zone of very large effective permeability ke. The flow pattern in the core when acid is being pumped (from left to right as shown in the figure) can therefore be represented as in FIG. 6.
The zone of high fluid mobility [26] can be parameterized by an effective fluid mobility Me=Mw and a propagation velocity Vw. Equivalently, the zone can also be characterized by an effective fluid viscosity μw or an effective permeability kw, derived according to equation (4).
Similarly, the zone of resistance or low fluid mobility [28] can be parameterized by an effective fluid mobility Me=Mr (and therefore according to Equation 4 an effective fluid viscosity μe=μr or an effective permeability ke=kr), as well as a propagation velocity Vr. Finally, there is a zone of displaced fluid [30] that was originally in the core prior to injection.
The velocities can be determined as follows




  
    
      
        {
        
          
            
              
                
                  
                    V
                    w
                  
                  ⁡
                  
                    (
                    
                      
                        (
                        
                          Q
                          /
                          A
                        
                        )
                      
                      ,
                      T
                      ,
                      Ro
                      ,
                      Ac
                    
                    )
                  
                
                =
                
                  
                    (
                    
                      Q
                      A
                    
                    )
                  
                  ⁢
                  
                    1
                    
                      
                        θ
                        0
                      
                      ⁡
                      
                        (
                        
                          
                            (
                            
                              Q
                              /
                              A
                            
                            )
                          
                          ,
                          T
                          ,
                          Ro
                          ,
                          Ac
                        
                        )
                      
                    
                  
                
              
            
          
          
            
              
                
                  
                    V
                    r
                  
                  ⁡
                  
                    (
                    
                      
                        (
                        
                          Q
                          /
                          A
                        
                        )
                      
                      ,
                      T
                      ,
                      Ro
                      ,
                      Ac
                    
                    )
                  
                
                =
                
                  
                    (
                    
                      Q
                      A
                    
                    )
                  
                  ⁢
                  
                    1
                    
                      
                        θ
                        r
                      
                      ⁡
                      
                        (
                        
                          
                            (
                            
                              Q
                              /
                              A
                            
                            )
                          
                          ,
                          T
                          ,
                          Ro
                          ,
                          Ac
                        
                        )
                      
                    
                  
                
              
            
          
        
      
    
    
      
        (
        7
        )
      
    
  




The parentheses indicate that the velocities and pore volumes to breakthrough are themselves functions of fluid velocity Q/A, temperature T, rock formation Ro, and acid formulation Ac. The functions ⊖0 and ⊖r are determined experimentally from the core flood experiments.
Using effective viscosities to express the effective mobilities, and rearranging the formulae, the effective viscosity μr is given by (8), and the derivation of (8) is given below.




  
    
      
        
          μ
          r
        
        =
        
          
            μ
            d
          
          ⁢
          
            
              Δ
              ⁢
              
                
              
              ⁢
              
                p
                r
              
            
            
              Δ
              ⁢
              
                
              
              ⁢
              
                p
                0
              
            
          
          ⁢
          
            
              Θ
              0
            
            
              
                Θ
                0
              
              -
              
                Θ
                r
              
            
          
        
      
    
    
      
        (
        8
        )
      
    
  




Where μd is the viscosity of the displaced fluid, originally saturating the core before acid is injected; Δp0 is the value of the pressure drop across the core when only the displaced fluid is pumped at the same conditions (typically brine). (8) is derived as follows. Let Lw be the distance traveled by the wormholes, measured from the core inlet, during the core-flood experiment, where the fluid mobility is Mw (see FIG. 6). Let Lr be the distance traveled by the front of low fluid mobility, where the fluid mobility is Mr (see FIG. 6). At the moment when Lr=L, L being the length of the core, Δpr is measured and using Darcy's law, we find that,





  
    
      
        
          
            
              
                Δ
                ⁢
                
                  
                
                ⁢
                
                  p
                  r
                
              
              =
              
                Q
                ⁢
                
                  
                    μ
                    r
                  
                  
                    Ak
                    0
                  
                
                ⁢
                
                  (
                  
                    L
                    -
                    
                      L
                      w
                    
                  
                  )
                
              
            
          
        
        
          
            
              
                =
                
                  Q
                  ⁢
                  
                    
                      μ
                      r
                    
                    
                      Ak
                      0
                    
                  
                  ⁢
                  
                    L
                    ⁡
                    
                      (
                      
                        1
                        -
                        
                          
                            Θ
                            r
                          
                          
                            Θ
                            0
                          
                        
                      
                      )
                    
                  
                
              
              ,
            
          
        
      
    
    
      
        (
        9
        )
      
    
  




and since, by definition,





  
    
      
        
          
            Δ
            ⁢
            
              
            
            ⁢
            
              p
              0
            
          
          =
          
            Q
            ⁢
            
              
                μ
                d
              
              
                Ak
                0
              
            
            ⁢
            L
          
        
        ,
      
    
    
      
        (
        10
        )
      
    
  




we then find (8) by simple algebra.

For the zone of high fluid mobility, we find that the effective fluid viscosity μe=μw in this region can be expressed as:




  
    
      
        
          μ
          w
        
        =
        
          
            μ
            d
          
          ⁢
          
            
              Δ
              ⁢
              
                
              
              ⁢
              
                p
                bt
              
            
            
              Δ
              ⁢
              
                
              
              ⁢
              
                p
                0
              
            
          
        
      
    
    
      
        (
        11
        )
      
    
  




where ΔPbt is the value of μp when the wormholes have broken through the outlet face of the core (this is the final value of Δp). (11) is derived as follows. When, Lw=L, L being the length of the core, Δpbt is measured. Using Darcy's law, we then find that,





  
    
      
        
          Δ
          ⁢
          
            
          
          ⁢
          
            p
            bt
          
        
        =
        
          Q
          ⁢
          
            
              μ
              w
            
            
              Ak
              0
            
          
          ⁢
          L
        
      
    
    
      
        (
        12
        )
      
    
  




then, using (10) and (12), we find (11) by simple algebra.

Equivalently, (8) and (11) can be used to define an effective mobility or an effective permeability in each zone, using Equation (4). This leads to equation (13).




  
    
      
        
          
            
              {
              
                
                  
                    
                      
                      ⁢
                      
                        
                          M
                          r
                        
                        =
                        
                          
                            k
                            0
                          
                          
                            
                              μ
                              d
                            
                            ⁢
                            
                              
                                Δ
                                ⁢
                                
                                  
                                
                                ⁢
                                
                                  p
                                  r
                                
                              
                              
                                Δ
                                ⁢
                                
                                  
                                
                                ⁢
                                
                                  p
                                  0
                                
                              
                            
                            ⁢
                            
                              
                                Θ
                                0
                              
                              
                                
                                  Θ
                                  0
                                
                                -
                                
                                  Θ
                                  r
                                
                              
                            
                          
                        
                      
                    
                  
                
                
                  
                    
                      
                      ⁢
                      
                        
                          M
                          w
                        
                        =
                        
                          
                            k
                            0
                          
                          
                            
                              μ
                              d
                            
                            ⁢
                            
                              
                                Δ
                                ⁢
                                
                                  
                                
                                ⁢
                                
                                  p
                                  bt
                                
                              
                              
                                Δ
                                ⁢
                                
                                  
                                
                                ⁢
                                
                                  p
                                  0
                                
                              
                            
                          
                        
                      
                    
                  
                
              
            
          
          
            
              {
              
                
                  
                    
                      
                      ⁢
                      
                        
                          k
                          r
                        
                        =
                        
                          
                            k
                            0
                          
                          
                            
                              
                                μ
                                d
                              
                              μ
                            
                            ⁢
                            
                              
                                Δ
                                ⁢
                                
                                  
                                
                                ⁢
                                
                                  p
                                  r
                                
                              
                              
                                Δ
                                ⁢
                                
                                  
                                
                                ⁢
                                
                                  p
                                  0
                                
                              
                            
                            ⁢
                            
                              
                                Θ
                                0
                              
                              
                                
                                  Θ
                                  0
                                
                                -
                                
                                  Θ
                                  r
                                
                              
                            
                          
                        
                      
                    
                  
                
                
                  
                    
                      
                      ⁢
                      
                        
                          k
                          w
                        
                        =
                        
                          
                            k
                            0
                          
                          
                            
                              
                                μ
                                d
                              
                              μ
                            
                            ⁢
                            
                              
                                Δ
                                ⁢
                                
                                  
                                
                                ⁢
                                
                                  p
                                  bt
                                
                              
                              
                                Δ
                                ⁢
                                
                                  
                                
                                ⁢
                                
                                  p
                                  0
                                
                              
                            
                          
                        
                      
                    
                  
                
              
            
          
        
      
    
    
      
        (
        13
        )
      
    
  




The use of Equations (8) and (11) in the case of axisymmetric radial flow around the wellbore in the reservoir as illustrated in FIGS. 7A and 7B. A wellbore [32] passes through a reservoir [34] and connects first to a wormholed or dissolved zone [36], bounded by a wormhole tip or dissolution front [38], and then to a resistance zone [40], bounded by a resistance zone front [42].
In FIGS. 7A and 7B, q(z,t) is the flow-rate per unit height into the reservoir at a time t, at a distance z along the well-bore. Let rw(z,t) be the radius of the wormhole-tip front or dissolution front and let rr(z,t) be the radius of the front of the resistance zone, both at the same time t and depth z. The evolution with time of both radii is then determined by solving the following set of equations.




  
    
      
        {
        
          
            
              
                
                ⁢
                
                  
                    
                      ∂
                      
                        ∂
                        t
                      
                    
                    ⁢
                    
                      (
                      
                        
                          r
                          w
                        
                        ⁡
                        
                          (
                          
                            z
                            ,
                            t
                          
                          )
                        
                      
                      )
                    
                  
                  =
                  
                    
                      
                        V
                        w
                      
                      ⁡
                      
                        (
                        
                          
                            V
                            ⁡
                            
                              (
                              
                                z
                                ,
                                
                                  r
                                  w
                                
                              
                              )
                            
                          
                          ,
                          
                            T
                            ⁡
                            
                              (
                              
                                z
                                ,
                                
                                  r
                                  w
                                
                              
                              )
                            
                          
                          ,
                          
                            Ro
                            ⁡
                            
                              (
                              
                                z
                                ,
                                
                                  r
                                  w
                                
                              
                              )
                            
                          
                          ,
                          
                            Ac
                            ⁡
                            
                              (
                              
                                z
                                ,
                                
                                  r
                                  w
                                
                              
                              )
                            
                          
                        
                        )
                      
                    
                    
                      
                        Φ
                        0
                      
                      ⁡
                      
                        (
                        
                          z
                          ,
                          
                            r
                            w
                          
                        
                        )
                      
                    
                  
                
              
            
          
          
            
              
                
                ⁢
                
                  
                    V
                    ⁡
                    
                      (
                      
                        z
                        ,
                        
                          r
                          w
                        
                      
                      )
                    
                  
                  =
                  
                    
                      q
                      ⁡
                      
                        (
                        
                          z
                          ,
                          t
                        
                        )
                      
                    
                    
                      2
                      ⁢
                      
                        
                      
                      ⁢
                      π
                      ⁢
                      
                        
                      
                      ⁢
                      
                        
                          r
                          w
                        
                        ⁡
                        
                          (
                          
                            z
                            ,
                            t
                          
                          )
                        
                      
                    
                  
                
              
            
          
        
      
    
    
      
        (
        14
        )
      
    
  
  
    
      
        {
        
          
            
              
                
                ⁢
                
                  
                    
                      ∂
                      
                        ∂
                        t
                      
                    
                    ⁢
                    
                      (
                      
                        
                          r
                          r
                        
                        ⁡
                        
                          (
                          
                            z
                            ,
                            t
                          
                          )
                        
                      
                      )
                    
                  
                  =
                  
                    
                      
                        V
                        r
                      
                      ⁡
                      
                        (
                        
                          
                            V
                            ⁡
                            
                              (
                              
                                z
                                ,
                                
                                  r
                                  r
                                
                              
                              )
                            
                          
                          ,
                          
                            T
                            ⁡
                            
                              (
                              
                                z
                                ,
                                
                                  r
                                  r
                                
                              
                              )
                            
                          
                          ,
                          
                            Ro
                            ⁡
                            
                              (
                              
                                z
                                ,
                                
                                  r
                                  r
                                
                              
                              )
                            
                          
                          ,
                          
                            Ac
                            ⁡
                            
                              (
                              
                                z
                                ,
                                
                                  r
                                  r
                                
                              
                              )
                            
                          
                        
                        )
                      
                    
                    
                      
                        Φ
                        0
                      
                      ⁡
                      
                        (
                        
                          z
                          ,
                          
                            r
                            r
                          
                        
                        )
                      
                    
                  
                
              
            
          
          
            
              
                
                ⁢
                
                  
                    V
                    ⁡
                    
                      (
                      
                        z
                        ,
                        
                          r
                          r
                        
                      
                      )
                    
                  
                  =
                  
                    
                      q
                      ⁡
                      
                        (
                        
                          z
                          ,
                          t
                        
                        )
                      
                    
                    
                      2
                      ⁢
                      
                        
                      
                      ⁢
                      π
                      ⁢
                      
                        
                      
                      ⁢
                      
                        
                          r
                          r
                        
                        ⁡
                        
                          (
                          
                            z
                            ,
                            t
                          
                          )
                        
                      
                    
                  
                
              
            
          
        
      
    
    
      
        (
        15
        )
      
    
  




Equations (14) and (15) are integrated by numerical means. Solving (14) and (15) allows the tracking of the wormhole tip and low-mobility front, respectively. In order to compute the pressure profile in the treated zone, i.e. at any z and for r between rwb and rr, (rwb is the wellbore radius at the depth z and therefore the pressure in the wellbore during the treatment, we make use of μr as follows:




  
    
      
        {
        
          
            
              
                
                  V
                  ⁡
                  
                    (
                    
                      z
                      ,
                      r
                      ,
                      t
                    
                    )
                  
                
                =
                
                  
                    
                      q
                      ⁡
                      
                        (
                        
                          z
                          ,
                          t
                        
                        )
                      
                    
                    
                      2
                      ⁢
                      
                        
                      
                      ⁢
                      π
                      ⁢
                      
                        
                      
                      ⁢
                      r
                    
                  
                  =
                  
                    
                      -
                      
                        
                          
                            k
                            e
                          
                          ⁡
                          
                            (
                            
                              z
                              ,
                              r
                              ,
                              t
                            
                            )
                          
                        
                        
                          
                            μ
                            e
                          
                          ⁡
                          
                            (
                            
                              z
                              ,
                              r
                              ,
                              t
                            
                            )
                          
                        
                      
                    
                    ⁢
                    
                      ∂
                      
                        ∂
                        r
                      
                    
                    ⁢
                    
                      p
                      ⁡
                      
                        (
                        
                          z
                          ,
                          r
                          ,
                          t
                        
                        )
                      
                    
                  
                
              
            
          
          
            
              
                
                  
                    μ
                    e
                  
                  ⁡
                  
                    (
                    
                      z
                      ,
                      r
                      ,
                      t
                    
                    )
                  
                
                =
                
                  {
                  
                    
                      
                        μ
                      
                      
                        
                          
                            if
                            ⁢
                            
                              
                            
                            ⁢
                            
                              
                                r
                                wb
                              
                              ⁡
                              
                                (
                                z
                                )
                              
                            
                          
                          <
                          r
                          <
                          
                            
                              r
                              w
                            
                            ⁡
                            
                              (
                              
                                z
                                ,
                                t
                              
                              )
                            
                          
                        
                      
                    
                    
                      
                        
                          μ
                          r
                        
                      
                      
                        
                          
                            if
                            ⁢
                            
                              
                            
                            ⁢
                            
                              
                                r
                                w
                              
                              ⁡
                              
                                (
                                
                                  z
                                  ,
                                  t
                                
                                )
                              
                            
                          
                          <
                          r
                          <
                          
                            
                              r
                              r
                            
                            ⁡
                            
                              (
                              
                                z
                                ,
                                t
                              
                              )
                            
                          
                        
                      
                    
                  
                
              
            
          
          
            
              
                
                  
                    k
                    e
                  
                  ⁡
                  
                    (
                    
                      z
                      ,
                      r
                      ,
                      t
                    
                    )
                  
                
                =
                
                  {
                  
                    
                      
                        
                          k
                          w
                        
                      
                      
                        
                          
                          ⁢
                          
                            
                              if
                              ⁢
                              
                                
                              
                              ⁢
                              
                                
                                  r
                                  wb
                                
                                ⁡
                                
                                  (
                                  z
                                  )
                                
                              
                            
                            <
                            r
                            <
                            
                              
                                r
                                w
                              
                              ⁡
                              
                                (
                                
                                  z
                                  ,
                                  t
                                
                                )
                              
                            
                          
                        
                      
                    
                    
                      
                        
                          k
                          0
                        
                      
                      
                        
                          
                          ⁢
                          
                            
                              if
                              ⁢
                              
                                
                              
                              ⁢
                              
                                
                                  r
                                  w
                                
                                ⁡
                                
                                  (
                                  
                                    z
                                    ,
                                    t
                                  
                                  )
                                
                              
                            
                            <
                            r
                          
                        
                      
                    
                  
                
              
            
          
        
      
    
    
      
        (
        16
        )
      
    
  




Equations (14)-(16) are integrated by analytical or numerical means and allow calculation of the pressure drop between the wellbore and rr, anywhere along the wellbore. The pressure at the wellbore p(z,rwb,t) can be determined from the pressure p(z,rr,t) at the resistance front using the following formula.




  
    
      
        {
        
          
            
              
                
                  p
                  ⁡
                  
                    (
                    
                      z
                      ,
                      
                        r
                        wb
                      
                      ,
                      t
                    
                    )
                  
                
                =
                
                  
                    p
                    ⁡
                    
                      (
                      
                        z
                        ,
                        
                          r
                          w
                        
                        ,
                        t
                      
                      )
                    
                  
                  +
                  
                    
                      ln
                      ⁡
                      
                        (
                        
                          
                            r
                            w
                          
                          
                            r
                            wb
                          
                        
                        )
                      
                    
                    ⁢
                    
                      
                        
                          q
                          ⁡
                          
                            (
                            
                              z
                              ,
                              t
                            
                            )
                          
                        
                        ⁢
                        μ
                      
                      
                        2
                        ⁢
                        
                          
                        
                        ⁢
                        π
                        ⁢
                        
                          
                        
                        ⁢
                        
                          k
                          w
                        
                      
                    
                  
                
              
            
          
          
            
              
                
                  p
                  ⁡
                  
                    (
                    
                      z
                      ,
                      
                        r
                        w
                      
                      ,
                      t
                    
                    )
                  
                
                =
                
                  
                    p
                    ⁡
                    
                      (
                      
                        z
                        ,
                        
                          r
                          r
                        
                        ,
                        t
                      
                      )
                    
                  
                  +
                  
                    
                      ln
                      ⁡
                      
                        (
                        
                          
                            r
                            r
                          
                          
                            r
                            w
                          
                        
                        )
                      
                    
                    ⁢
                    
                      
                        
                          q
                          ⁡
                          
                            (
                            
                              z
                              ,
                              t
                            
                            )
                          
                        
                        ⁢
                        
                          μ
                          r
                        
                      
                      
                        2
                        ⁢
                        
                          
                        
                        ⁢
                        π
                        ⁢
                        
                          
                        
                        ⁢
                        
                          k
                          0
                        
                      
                    
                  
                
              
            
          
        
      
    
    
      
        (
        17
        )
      
    
  




In (16) and (17), it is possible to substitute the effective viscosity μr and the effective permeability kw with other combinations giving rise to the same fluid mobility, for instance, (16) is equivalent to (18) and (17) to (19).





  
    
      
        {
        
          
            
              
                
                  V
                  ⁡
                  
                    (
                    
                      z
                      ,
                      r
                      ,
                      t
                    
                    )
                  
                
                =
                
                  
                    
                      q
                      ⁡
                      
                        (
                        
                          z
                          ,
                          t
                        
                        )
                      
                    
                    
                      2
                      ⁢
                      
                        
                      
                      ⁢
                      π
                      ⁢
                      
                        
                      
                      ⁢
                      r
                    
                  
                  =
                  
                    
                      -
                      
                        M
                        ⁡
                        
                          (
                          
                            z
                            ,
                            r
                            ,
                            t
                          
                          )
                        
                      
                    
                    ⁢
                    
                      ∂
                      
                        ∂
                        r
                      
                    
                    ⁢
                    
                      p
                      ⁡
                      
                        (
                        
                          z
                          ,
                          r
                          ,
                          t
                        
                        )
                      
                    
                  
                
              
            
          
          
            
              
                
                  M
                  ⁡
                  
                    (
                    
                      z
                      ,
                      r
                      ,
                      t
                    
                    )
                  
                
                =
                
                  {
                  
                    
                      
                        
                          M
                          w
                        
                      
                      
                        
                          
                            if
                            ⁢
                            
                              
                            
                            ⁢
                            
                              
                                r
                                wb
                              
                              ⁡
                              
                                (
                                z
                                )
                              
                            
                          
                          <
                          r
                          <
                          
                            
                              r
                              w
                            
                            ⁡
                            
                              (
                              
                                z
                                ,
                                t
                              
                              )
                            
                          
                        
                      
                    
                    
                      
                        
                          M
                          r
                        
                      
                      
                        
                          
                            if
                            ⁢
                            
                              
                            
                            ⁢
                            
                              
                                r
                                w
                              
                              ⁡
                              
                                (
                                
                                  z
                                  ,
                                  t
                                
                                )
                              
                            
                          
                          <
                          r
                          <
                          
                            
                              r
                              r
                            
                            ⁡
                            
                              (
                              
                                z
                                ,
                                t
                              
                              )
                            
                          
                        
                      
                    
                  
                
              
            
          
        
      
    
    
      
        (
        18
        )
      
    
  
  
    
      
        {
        
          
            
              
                
                  p
                  ⁡
                  
                    (
                    
                      z
                      ,
                      
                        r
                        wb
                      
                      ,
                      t
                    
                    )
                  
                
                =
                
                  
                    p
                    ⁡
                    
                      (
                      
                        z
                        ,
                        
                          r
                          w
                        
                        ,
                        t
                      
                      )
                    
                  
                  +
                  
                    
                      ln
                      ⁡
                      
                        (
                        
                          
                            r
                            w
                          
                          
                            r
                            wb
                          
                        
                        )
                      
                    
                    ⁢
                    
                      
                        q
                        ⁡
                        
                          (
                          
                            z
                            ,
                            t
                          
                          )
                        
                      
                      
                        2
                        ⁢
                        
                          
                        
                        ⁢
                        π
                        ⁢
                        
                          
                        
                        ⁢
                        
                          M
                          w
                        
                      
                    
                  
                
              
            
          
          
            
              
                
                  p
                  ⁡
                  
                    (
                    
                      z
                      ,
                      
                        r
                        w
                      
                      ,
                      t
                    
                    )
                  
                
                =
                
                  
                    p
                    ⁡
                    
                      (
                      
                        z
                        ,
                        
                          r
                          r
                        
                        ,
                        t
                      
                      )
                    
                  
                  +
                  
                    
                      ln
                      ⁡
                      
                        (
                        
                          
                            r
                            r
                          
                          
                            r
                            w
                          
                        
                        )
                      
                    
                    ⁢
                    
                      
                        q
                        ⁡
                        
                          (
                          
                            z
                            ,
                            t
                          
                          )
                        
                      
                      
                        2
                        ⁢
                        
                          
                        
                        ⁢
                        π
                        ⁢
                        
                          
                        
                        ⁢
                        
                          M
                          r
                        
                      
                    
                  
                
              
            
          
        
      
    
    
      
        (
        19
        )
      
    
  




FIGS. 8 and 9 illustrate in a physical way the process described above. To illustrate, an experiment is conducted whereby acid (e.g. 15% HCl) is pumped from the top into a cylindrical core [6] held between two seals [44] as shown in FIG. 8. During the acid injection, performed at a constant flow-rate, the pressure difference between the wellbore [32] and the periphery of the core [46] is logged. The pressure drop is a direct indication of the distance traveled by the wormholes during this experiment. The agreement between the result predicted by the method and the experimental one is very good.
The procedural techniques for pumping stimulation fluids down a wellbore to acidize a subterranean formation are well known. The person who designs such matrix acidizing treatments has available many useful tools to help design and implement the treatments, one of which is a computer program commonly referred to as an acid placement simulation model (a.k.a., matrix acidizing simulator, wormhole model). Most if not all commercial service companies that provide matrix acidizing services to the oilfield have one or more simulation models that their treatment designers use. One commercial matrix acidizing simulation model that is widely used by several service companies is known as StimCADETM. This commercial computer program is a matrix acidizing design, prediction, and treatment-monitoring program that was designed by Schlumberger Technology Corporation. All of the various simulation models use information available to the treatment designer concerning the formation to be treated and the various treatment fluids (and additives) in the calculations, and the program output is a pumping schedule that is used to pump the stimulation fluids into the wellbore. The text “Reservoir Stimulation,” Third Edition, Edited by Michael J. Economides and Kenneth G. Nolte, Published by John Wiley & Sons, (2000), is an excellent reference book for matrix acidizing and other well treatments.
As previously mentioned, because the ultimate goal of matrix acidizing is to alter fluid flow in a reservoir, reservoir engineering must provide the goals for a design. In addition, reservoir variables may impact the treatment performance.
In various embodiments, the overall procedure is implemented into an acid placement simulator to predict the fate of a given design in the field.
A global methodology used by field engineers is described in FIG. 10:
The optimization in FIG. 10 makes use of the above methodology to predict a given acid treatment performance. It is possible to improve a design by

    
    
        Changing operational parameters such as:
        
            Pumping rate
            Acid volume
            Acid formulation
            Number of acid stages
        
        
        Understanding important parameters controlling the treatment outcome such as:
        
            Operational parameters
            Reservoir parameters
            Wellbore parameters
            Conveyance parameters
        
        
    
    


EXAMPLES 
A computer program has been developed to simulate the injection of acid into a carbonate reservoir. The simulator inputs include all the relevant reservoir parameters, schedule and fluid parameters.

    
    
        The simulator predicts the flow of the pumped fluids down the wellbore: location, concentration of acid along the wellbore vs. time and pressure distribution along the wellbore. This is done by mass conservation principle and by using hydrostatic and friction pressure models.
        The wellbore is connected to the reservoir and fluid from the wellbore will flow into the various reservoir layers if the pressure in the wellbore exceeds the pore-pressure in the reservoir. The initial pore pressure is a user input.
        Once the stimulation fluids enter the reservoir at any given depth z along the wellbore, the dissolution fronts (also referred here as the high-mobility front or wormhole-tip front) at this depth, as well as the front of the zone of low mobility, if a self-diverting acid is being pumped) are tracked using equations (14) and (15).
        The effect of acid flowing into the reservoir is to change the fluid mobility distribution and, therefore, the pressure in the reservoir changes. The pore pressure is updated using equations (16) and (17).
        For the two above calculations to be possible, the flow-rate q must be known at the depth z under consideration. The flow-rate q can be estimated using equations (16) and (17) or any equivalent formulations before updating the fluid mobility distribution in the reservoir.
        Then, the location of the dissolution fronts are updated over a certain time-step and the calculations are repeated until the full treatment schedule is complete.
    
    


An example is given in FIG. 11: A well, partly deviated, is to be stimulated. The reservoir from which the well is producing is a limestone reservoir with three producing layers of 100, 20 and 5 mD as depicted in FIG. 12. The dimensions of the layers as well as their petrophysical properties are input into the simulator. These include

    
    
        Permeability, porosity
        Layer fluid saturations and fluid properties
        Layer dimensions, temperatures and pore pressures
        Drilling damage characteristics: skin and depth for each layer
    
    


The well trajectory and dimensions are also input into the simulator. Finally, the type of completion used for this well is also input, in this case the wellbore is open-hole (no casing). The engineer's task is to design the best possible treatment. In other words, the engineer task is to ensure that he delivers the treatment the will provide the best stimulation given some economical and operational constraints.
First, acid core flood experiments, as described above, are performed using core samples from the layers of interest. These are used to calibrate the correlations for θr and μr. θ0 is also determined. These tests are performed at the reservoir temperature, for various rates, and for the candidate stimulation fluids, in this case, 15% HCl and 15% VDATM. The parameters θr, μr and θ0 are tabulated versus flux (V=q/A) and input into the simulator for the various flow-rates tested during the experiment. These tables, or correlations if correlations have been derived, are used in connection with equations (14)-(17) in order to predict the position of the front of the zone of high fluid mobility (where wormholes have increased the virgin permeability) and that of the zone of low fluid mobility.
The task now consists of optimizing acid volumes and rates in order to achieve an optimum treatment. Treatment efficiency is measured by comparing the wellbore skin before and after treatment. The further the wormholes extend into the layers, the lower the wellbore skin and the higher the production rate after treatment.
For such wells, a typical treatment consists of bullheading 15% HCl from the well-head at a constant rate. Given some operational constraints, the rate has to be between 0.5 bbl/min and 5 bbl/min in this example. For economical reasons, only 75 gal/ft of acid will be pumped. The first optimization step consists of running the simulator with different injection rates and choose that one providing the best treatment, with 15% HCl, the most economical acid system. The results are represented in FIG. 12A-12D. It is possible to see that the wormholes extended deeper into the top most-permeable layer of the reservoir that into the middle layer. The lower-permeability zone at the bottom does not get any stimulation. The best treatment with HCl is when the later is pumped at 5 bbl/min. The second step is to do the same exercise with 15% VDA. The results are represented in FIG. 13A-13D. Though wormholes do not propagate as far as with HCI in the top layer, the use of VDA pumped at 5 bbl/min shows that the zonal coverage is better and all layers show similar treatment depth. Because of the good zonal coverage and deep enough wormhole penetration (beyond the damage depth), the preferred treatment consists of pumping 15% VDA at 5 bbl/min. FIG. 14 also illustrates the position of the fronts of the zones of low fluid mobility, where M=Mr, responsible for the diversion.

## Claims

1. A method of modeling the pressure in a wellbore during acid treatment with a self diverting acid delivered at a flowrate Q, the pressure being determined at a depth z, a distance r from the center of the well, and a time t, the method involving use of functions derived from core flooding experiments wherein a self diverting acid is injected into a core and the pressure along the core is measured as a function of time, the modeling method comprising:
calculating at least one of an effective viscosity μr, a mobility Mr, and a permeability kr, wherein




  
    
      μ
      r
    
    =
    
      
        μ
        d
      
      ⁢
      
        
          Δ
          ⁢
          
            
          
          ⁢
          
            p
            r
          
        
        
          Δ
          ⁢
          
            
          
          ⁢
          
            p
            o
          
        
      
      ⁢
      
        
          Θ
          0
        
        
          
            Θ
            o
          
          -
          
            Θ
            r
          
        
      
    
  




  
    
      M
      r
    
    =
    
      
        k
        0
      
      
        
          μ
          d
        
        ⁢
        
          
            Δ
            ⁢
            
              
            
            ⁢
            
              p
              r
            
          
          
            Δ
            ⁢
            
              
            
            ⁢
            
              p
              0
            
          
        
        ⁢
        
          
            Θ
            0
          
          
            
              Θ
              0
            
            -
            
              Θ
              r
            
          
        
      
    
  




  and




  
    
      k
      r
    
    =
    
      
        k
        0
      
      
        
          
            μ
            d
          
          μ
        
        ⁢
        
          
            Δ
            ⁢
            
              
            
            ⁢
            
              p
              r
            
          
          
            Δ
            ⁢
            
              
            
            ⁢
            
              p
              0
            
          
        
        ⁢
        
          
            Θ
            0
          
          
            
              Θ
              0
            
            -
            
              Θ
              r
            
          
        
      
    
  



wherein:
k0 is the initial absolute permeability of the core, before acid is injected;
μd is the viscosity of the displaced fluid originally saturating the core before acid is injected;
μ is the viscosity of the acid;
Δpr is the pressure drop derived from the core flooding experiments and is the pressure drop at the time tr that the pressure drop changes from a first linear trend to a second linear trend;
⊖r is the number of pore volumes delivered to the core at the time tr;
Δpo is the pressure drop at t=o of the core flood experiment; and
⊖o is the pore volume to breakthrough measured in the core flood experiment; and

using a simulator to calculate pressures within the formation on the basis of the effective viscosity μr, the mobility Mr, and/or the permeability kr.

2. The method according to claim 1, comprising deducing the pressure at the wellbore p(z, rwb, t) from the pressure at the resistance front p(z, rr, t) from




  
    p
    ⁡
    
      (
      
        z
        ,
        
          r
          wb
        
        ,
        t
      
      )
    
  
  =
  
    
      p
      ⁡
      
        (
        
          z
          ,
          
            r
            w
          
          ,
          t
        
        )
      
    
    +
    
      
        ln
        ⁡
        
          (
          
            
              r
              
                w
                ⁢
                
                  
                
              
            
            
              r
              wb
            
          
          )
        
      
      ⁢
      
        
          
            q
            ⁡
            
              (
              
                z
                ,
                t
              
              )
            
          
          ⁢
          
            μ
            w
          
        
        
          2
          ⁢
          
            
          
          ⁢
          π
          ⁢
          
            
          
          ⁢
          
            k
            0
          
        
      
    
  




wherein




  
    
      p
      ⁡
      
        (
        
          z
          ,
          
            r
            w
          
          ,
          t
        
        )
      
    
    =
    
      
        p
        ⁡
        
          (
          
            z
            ,
            
              r
              r
            
            ,
            t
          
          )
        
      
      +
      
        
          ln
          ⁡
          
            (
            
              
                r
                
                  r
                  ⁢
                  
                    
                  
                
              
              
                r
                w
              
            
            )
          
        
        ⁢
        
          
            
              q
              ⁡
              
                (
                
                  z
                  ,
                  t
                
                )
              
            
            ⁢
            
              μ
              r
            
          
          
            2
            ⁢
            
              
            
            ⁢
            π
            ⁢
            
              
            
            ⁢
            
              k
              0
            
          
        
      
    
  
  ;








  
    p
    ⁡
    
      (
      
        z
        ,
        
          r
          wb
        
        ,
        t
      
      )
    
  
  =
  
    
      p
      ⁡
      
        (
        
          z
          ,
          
            r
            w
          
          ,
          t
        
        )
      
    
    +
    
      
        ln
        ⁡
        
          (
          
            
              r
              
                w
                ⁢
                
                  
                
              
            
            
              r
              wb
            
          
          )
        
      
      ⁢
      
        
          q
          ⁡
          
            (
            
              z
              ,
              t
            
            )
          
        
        
          2
          ⁢
          
            
          
          ⁢
          π
          ⁢
          
            
          
          ⁢
          
            M
            w
          
        
      
    
  




wherein




  
    
      p
      ⁡
      
        (
        
          z
          ,
          
            r
            w
          
          ,
          t
        
        )
      
    
    =
    
      
        p
        ⁡
        
          (
          
            z
            ,
            
              r
              r
            
            ,
            t
          
          )
        
      
      +
      
        
          ln
          ⁡
          
            (
            
              
                r
                
                  r
                  ⁢
                  
                    
                  
                
              
              
                r
                w
              
            
            )
          
        
        ⁢
        
          
            q
            ⁡
            
              (
              
                z
                ,
                t
              
              )
            
          
          
            2
            ⁢
            
              
            
            ⁢
            π
            ⁢
            
              
            
            ⁢
            
              M
              r
            
          
        
      
    
  
  ;




or




  
    p
    ⁡
    
      (
      
        z
        ,
        
          r
          wb
        
        ,
        t
      
      )
    
  
  =
  
    
      p
      ⁡
      
        (
        
          z
          ,
          
            r
            w
          
          ,
          t
        
        )
      
    
    +
    
      
        ln
        ⁡
        
          (
          
            
              r
              
                w
                ⁢
                
                  
                
              
            
            
              r
              wb
            
          
          )
        
      
      ⁢
      
        
          
            q
            ⁡
            
              (
              
                z
                ,
                t
              
              )
            
          
          ⁢
          μ
        
        
          2
          ⁢
          
            
          
          ⁢
          π
          ⁢
          
            
          
          ⁢
          
            k
            w
          
        
      
    
  




wherein




  
    p
    ⁡
    
      (
      
        z
        ,
        
          r
          w
        
        ,
        t
      
      )
    
  
  =
  
    
      p
      ⁡
      
        (
        
          z
          ,
          
            r
            r
          
          ,
          t
        
        )
      
    
    +
    
      
        ln
        ⁡
        
          (
          
            
              r
              
                r
                ⁢
                
                  
                
              
            
            
              r
              w
            
          
          )
        
      
      ⁢
      
        
          
            q
            ⁡
            
              (
              
                z
                ,
                t
              
              )
            
          
          ⁢
          μ
        
        
          2
          ⁢
          
            
          
          ⁢
          π
          ⁢
          
            
          
          ⁢
          
            k
            r
          
        
      
    
  




wherein
z is the depth in the wellbore;
rwb is the radius of the wellbore at a depth z;
rw is the radius of the dissolution front or the zone of high fluid mobility
rr is the radius of the zone of resistance at a depth z and at a time t;
q is the flow rate of self-diverting acid in the formation at a depth z and at a time t;
μ is the viscosity of the self-diverting acid before the acid is spent; and
kw is the effective permeability in the region of high fluid mobility.

3. A method of optimizing acid treatment of a hydrocarbon containing carbonate reservoir with a self-diverting acid, comprising:
carrying out linear core flood experiments varying one or more parameters selected form the group consisting of acid formulation, rock type, flow rate, and temperature;
deriving the following functions from the experiments, as a function of the parameters:
⊖o—the pore volume to wormhole/dissolution front breakthrough;
⊖r—the pore volume to resistance zone breakthrough; and
Δpr—the pressure drop at resistance zone breakthrough;

writing equations of a flow model based on the functions;
solving the equations in an arbitrary flow field in a simulator;
using the simulator in an optimization loop together with known and/or estimated reservoir parameters; and
calculating at least one of the following from the simulator optimization loop:
stage and rate volumes of the acid treatment;
fluid selection for the acid treatment;
wormhole invasion profile; and
skin profile.

4. A method according to claim 3, comprising deriving fluid mobilities in the resistance zone and in a zone of large mobility on the basis of Darcy's law from measurements of pressure drop along the core in the core flood experiments.

5. A method according to claim 3, comprising optimizing at least one of
the pumping rate;
acid volume;
acid formulations; and
number of acid stages.