Analytical Modeling and Contradictions in Limestone Reservoirs : Breccias , Vugs , and Fractures

Modeling of limestone reservoirs is traditionally developed applying tectonic fractures concepts or planar discontinuities and has been simulated dynamically without considering nonplanar discontinuities as sedimentary breccias, vugs, fault breccias, and impact breccias, assuming that all these nonplanar discontinuities are tectonic fractures, causing confusion and contradictions in reservoirs characterization. The differences in geometry and connectivity in each discontinuity affect fluid flow, generating the challenge to develop specific analytical models that describe quantitatively hydrodynamic behavior in breccias, vugs, and fractures, focusing on oil flow in limestone reservoirs. This paper demonstrates the differences between types of discontinuities that affect limestone reservoirs and recommends that all discontinuities should be included in simulation and static-dynamic characterization, because they impact fluid flow. To demonstrate these differences, different analytic models are developed. Findings of this work are based on observations of cores, outcrops, and tomography and are validated with field data. The explanations and mathematical modeling developed here could be used as diagnostic tools to predict fluid velocity and fluid flow in limestone reservoirs, improving the complex reservoirs static-dynamic characterization.


Introduction
A significant portion of the world oil reservoirs is found in carbonates [1,2]; early identification of the types of geological discontinuities that are dominant in limestone rock is one of the principal problems for advanced development and efficient exploitation of the reservoirs, because of impact in the reservoir oil flow.These complex limestone reservoirs are associated with geological events (tectonic fractures, sedimentary breccias, impact breccias, vugs, and fault breccias) that require to be understood dynamically and geologically to properly derive an analytical model that would describe flow behavior.
The phenomenology of fluid flow in limestone carbonate reservoirs depends on type of discontinuities; connected discontinuities have a dependence on, at least, diagenesis, structural history, and lithology.A key concept is that planar or nonplanar discontinuities can lead to different flow response during oil production.Moreover, limestone carbonate reservoirs are modeled, confused, and associated with planar discontinuities named tectonic fracture networks causing contradictions.Processes that generate nonplanar discontinuities are usually excluded, and models do not represent the reservoir reality, when based on an equivalent flow medium.A case well-known is Cantarell field, which is the eighth largest oil field in the world, with fractures, breccias, and vugs [3][4][5] located in the Campeche Sound, Yucatán Platform in Mexico.This reservoir is modeled using tectonic fracture networks without considering vugs, fault and impact breccias, and caverns [6][7][8]; in effect, these tectonic fractures are used to represent all types of discontinuities.
In 1972, Neale and Nader [9] described the flow dynamics in a vuggy porous medium using the creeping Navier-Stokes and the Darcy equations in the surrounding homogeneous and isotropic medium.Koenraad and Bakker in 1981 [10] proposed a theoretical study about fracture/vug, collapse breccias, and brecciated karst based on geological information.Wu et al. in 2011 [11] proposed numerical models multiphase and single-phase flow for vuggy reservoirs using balance mass equations; Darcy and Forchheimer equations; and Hagen Poiseuille tube flow.In 1999, McKeown et al. [12] developed a numerical groundwater flow model for Sellafield, in which they modeled hydrogeological Brockram fault breccia using Darcy's equation applied to a limestone with porosity between 1 and 10 percent, obtaining a range of hydraulic conductivity.Gudmundsson [13] studied fluid flow behavior using Darcy equation applied to breccias and faults.
Sedimentary breccias are related to debris flow.The physics of debris flow published by Iverson in 1997 [14] analyzed flows of dry, granular solids and solid-fluid mixtures and described the grains behavior through momentum, mass, and energy balances.Enos [15] described debris reservoir and sedimentary breccias in Tamabra Limestone of the Poza Rica field in Mexico, with field ranges of porosity of 3.7-9.7%and permeability of 0.01-700 md with depth ranging between 1980 and 2700 m and cumulative production to July 1983 of 1.98 billion oil barrels.There are many impact breccias, but, associated with oil reservoir, the most known is KT limit in Cantarell located in Campeche Bay, Yucatan Platform in Mexico, related to the Chicxulub impact [16].From the point of view of fluid flow, Mayr et al. [17] in 2008 estimated rocks hydraulic permeability associated with the impact crater Chicxulub using Nuclear Magnetic Resonance, Kozeny-Carman equation, and the fractal PaRis-model obtaining low hydraulic permeability.
The purpose of this paper is to demonstrate how qualitative and quantitative differences between tectonic fractures, fault breccias, sedimentary breccias, vugs, and impact breccias affect the exploitation strategy and advanced development of carbonate reservoirs.Our hypothesis considers that all the discontinuities in carbonate reservoirs are distinct and cannot be called or analyzed as tectonic fractures without knowing their genesis and discontinuity impact in the reservoir, or only focusing in fluid flow.To demonstrate this hypothesis, we understood the rocks nature and used subsurface geological analogs (outcrops), cores, computerized tomography, and analytical modeling, making use of fluid flow kinematics.

First Contradiction: Tectonic Fracture or
Nonplanar Discontinuity?
Unfortunately, whatever the discontinuity (planar or nonplanar), it is considered as a natural fracture.Usually, oil reservoirs with discontinuities are "naturally fractured, " causing confusion in reservoirs simulation.Our discussion begins with the following: the word fracture derived from the Latin fractus means "broken" and its definitions are as follows.
"A natural fracture is a macroscopic planar discontinuity, which results from stresses that exceed the rupture strength of the rock" [18].
"A reservoir fracture is a naturally occurring macroscopic planar discontinuity in rock due to deformation or physical diagenesis" [19].
"A fracture is any planar or subplanar discontinuity that is very narrow in one dimension compared to the other two and forms as a result of external (e.g.tectonic) or internal (thermal or residual) stress" [20].
These definitions converge in that fractures result from stresses and deformation on the rock, mechanical diagenesis that generates volume reduction by compaction during burial.Whether natural tectonic fractures result as the effect of stresses and mechanical compaction, why are sedimentary breccia, impact breccia, fault breccia, and vug considered natural fractures regardless of the genesis of these geological events?This is the first contradiction.
The discontinuities observed in limestone reservoirs are planar or nonplanar.The observation of the core of limestone reservoirs in Figure 1 shows a planar natural fracture with inclination; it is a preferential path for oil flow; its aperture has been developed on the plane of weakness (with low cohesion) as consequence of stresses on the rock.Although the core was badly preserved, its natural fracture presents an aperture increase, but it is natural.
However, many reservoirs normally contain nonplanar discontinuities caused by diagenesis.The dissolution processes (chemical diagenesis) as result of expelled water generate vugs, caverns, and clasts dissolution.This type of discontinuities affects the fluid flow and creates an interconnected system through pores.The Figure 2(a) presents Cladocoropsis that could act as barrier because these fragments have not been fragmented or dissolved.Also, Figure 2(b) shows a core with fragmented Cladocoropsis grainstone due to chemical diagenesis, creating a secundary porosity and permeability, connecting pores system (primary and secondary porosity).
Reservoir engineers have assumed discontinuities as a type of pores regardless of their origin or genesis, and tectonic fractures are considered as structural discontinuities impacting fluid flow [21].So the discontinuities or heterogeneities such as vugs, tectonic fractures, and porosity associated with breccias can be interpreted with an equivalent medium [22].
In this case, "all is tectonic fractures" regardless of their genesis creating confusion.When we compare and observe Figure 1 with Figure 2, it can be inferred that in dissimilar media oil flows have different velocities.From a modeling point of view, the discontinuities (planar or nonplanar) in the reservoirs should be related to their genesis and fluid characteristics flow to understand the phenomenology of limestone systems.

Second Contradiction: How Do Different Discontinuities Impact?
Despite the fact that there exist different kinds of discontinuities, such as tectonic fractures, vugs, caves, collapse breccias, sedimentary breccias, fault breccias, and impact breccias, in transient pressure analysis and reservoir simulation all of them are considered as fractures or fissures.In consequence, they are represented as planes in the reservoir, with orientation, intensity, or fracturing frequency, with specific geometry, and with fluid dynamical properties.This study is based on open and connected discontinuities.The discontinuities clearly affect fluid flow in the reservoir.However, each discontinuity creates a particular hydrocarbons flow behavior, specifically, in limestone reservoirs.It is essential to understand them to properly study its influence on the production behavior of the reservoirs.
The phenomenology of discontinuities begins with their geometry.In principle, a planar geometry describes a parabolic flow profile for conditions of laminar flow for tectonic fractures; the fluid velocity is associated with the aperture and roughness of the discontinuity, which is linked to the permeability.The pressure profile changes are due to friction, tortuosity, the complex interconnected system, and cross flow.This phenomenon is present in tectonic fractures.
A nonplanar discontinuity shows nonparabolic flow profile because of high oil velocity and heterogeneous media.The fluid high velocity is directly linked with geometry, shape, aperture, roughness, diameter, discontinuities distribution, irregular surfaces, and pressure.In the case of vugs and caverns, their diameters or apertures are large and generate zones of high permeability, or superhighway production.
Discontinuities (planar or nonplanar) can be described through visible geological features related to geometry; specifically, cross section flow area, which has a direct impact in fluid flow.This can be understood using the flow concept  = .A nonplanar discontinuity area is greater than a planar discontinuity area; in effect, flow increases because pores may be connected and their size are increased due to dissolution.In addition, when the flow area is increased, oil velocity  increases too.This will be demonstrated using analytical models in this study.Moreover, it can be observed when vugs and caverns are compared with tectonic fractures that they have a greater permeability.In summary, fluid flow in a planar discontinuity (fracture) is different to flow in a nonplanar discontinuity (vugs, breccias, and caverns).This is the second contradiction, considering that hydrodynamics is intrinsic to geometry or flow area because there is a mass and momentum transfer process, oil flow in each discontinuity is unique and cannot be traditionally considered as fractures.The understanding and comparison of these dynamic and geological differences in limestone reservoirs allow the accuracy of fluid flow characteristics and approach their optimum exploitation strategy.
Although discontinuities generate highly heterogeneous and anisotropic systems, a complex reservoir could evidence different types of geological characteristics; production and pressure behavior might identify dominant discontinuities, and, to achieve a reasonable modeling, it should be considered that they might in some cases perform as seal or flow barriers.
However, before carrying out pressure transient test and reservoir dynamic simulations, it is necessary to characterize these discontinuities and establish convergence between geology and fluid flow features, using cores, outcrops, seismic, well logs, tomography, pressures profiles, production, and mathematical modeling that describe the flow kinematic in reservoir fluid flow.
It is becoming increasingly difficult to ignore the differences between types of discontinuities, in which all of them are called fractures by simplicity.Therefore, there is a need to be explicit about the phenomenology of discontinuities.

Analogs Reservoirs and Outcrops for a Flow Analytical Model
Reservoir analog based on outcrop studies provides definitions of geometry and heterogeneities at interwell scales and characterizes reservoirs [24].A subsurface reservoir model using outcrop analog [25] is employed for an understanding of geological parameters associated with reservoir and at outcrop rock (porosity, permeability, and heterogeneities) with the objective of characterizing it.Outcrops were at geological conditions of subsurface; fortunately they are exposed, and we may describe a reservoir with similar features.In this paper, we used outcrops to understand physics of each discontinuity (vugs, breccias, and fractures) and tomography to know how is the fluid flow through the rock.As a result of understanding that rock, we developed an analytical model.

Geological and Tomography Features of Tectonic Fractures
Limestone reservoirs constitute a portion of the basin.Tectonic fractures are planar discontinuities in reservoirs.Generally, the genesis of fracturing is tectonism and it is related to local folds, faults, or a regional system [26].There are different fractures types: shear fractures or slip surface, extension fractures such as joints, fissures, and veins, and contraction or closing fractures as stylolites [20].Open fractures can behave as hydraulic conductors and impact the productivity of the reservoir.Tectonic fractures have been observed in carbonate cores, but it is difficult to obtain complete cores due to its brittle nature (Figure 1).Tectonic fractures can be open, deformed, and mineral-filled.In spite of deformation, it is evident that tectonic fractures are planar discontinuities that could store and produce hydrocarbons.In outcrops, these fractures are studied because they explain paleostresses (Figure 3).Also, these surface systems are used as an analogical model in reservoir characterization.Figure 3 shows different planar fractures systems as a result of paleostresses in the rock, and it is observed that each fracture has its aperture that may be interconnected, inclined, parallel, orthogonal, open, or closed.This outcrop was buried at subsurface conditions during a geological age; nowadays it is exposed.In Figure 3, diverse types of natural fractures can be observed.
X-ray computerized tomography (CT) is a nondestructive technique that allows visualization of the internal structure of rocks, determined mainly by variations in density and atomic composition [27].One well-known advantage of tomography is fracture morphology description inside the rock.This implies that open fractures are observed as planar discontinuities with a specific physical behavior.In consequence, the dominant parameter in tectonic fractures is aperture.Aperture is associated with permeability, parabolic flow profile, maximum and mean velocity, and the flow rate.
A study illustrated how width or aperture impacts fracture permeability using an equation flow [28].It was shown that a single fracture of 0.25 mm aperture has the equivalent permeability of 188 m of unfractured rock, with a uniform permeability of 10 md, and a 1.27 mm aperture fracture is equivalent to 173 m of rock with a permeability of 1,000 md.
Tomography was used to show internal fractures and to verify their morphology.Useful, calcareous cores without dissolutions or visible discontinuities are considered as matrix rock without fractures.In many cases, narrow subplanar fractures can be observed using tomography or would be inferred if the density in the core internal structure shows changes.
A sample core in a limestone reservoir (15.9 cm in length and 10 cm in diameter) was obtained.The scan of a limestone core (Figure 4) at 3 mm spacing showed that rock contains narrow subplanar fractures.They are open, deformed, and/or mineral-filled.In images 16 to 21 an open fracture with great aperture can be seen, and in images 30 to 42 other fractures can be observed, but with small aperture.Moreover, doing a visual inspection of core does not observe any fractures that may be connected (see slides 30 to 53 in Figure 4).Approximately, these limited length fractures have 1 to 2 mm apertures and would allow hydrocarbons flow; in effect they also provide effective porosity because fracture connects spaces in this rock.

Analytical Model for the Profile Velocity of Tectonic Fractures
Fluid flow is defined by a velocity vector field and a pressure scalar field.In order to understand how permeability influences fluid flow in fractures, we analyze flow and fluid kinematics equations.Kinematic equations are applied to an isothermal, low viscosity, irrotational, and single-phase fluid flow and are referred to as potential flows and stream function.Applicability of kinematic equations is limited to Reynolds numbers, Re, more than unity because the inviscid fluid flow theory corresponds to high flow values and small fluid viscosity.Calcareous reservoirs heterogeneities may produce high velocity flow.At larger Reynolds numbers, kinematic equations are suggested for Newtonian liquids of minimal viscosity and gases.
Fluid kinematics describes fluid motion using streamlines and potential function, with respect to a plane, and a flow equation describes a close relationship between fracture aperture, pressure, and velocity and flow rate, using classical methods in fluid mechanics.Kinematics of fluid flow would show a particular behavior of flow lines for each discontinuity, in this case for tectonic fractures (Figure 5(a)).Classical methods in fluid mechanics are illustrated by an analytical model for velocity profile in parallel plates [29].Two symmetrically inclined parallel surfaces with respect to the -axis are shown in Figure 5(b).Fluid motion is in the  direction and flow velocity distribution in the  direction, for steady state of an incompressible fluid of constant density  is indicated in this figure; pressure variation is a function of  coordinate, and the surface longitudes are larger than the aperture.
Figure 5 shows top and lateral views in an open fracture simulated with two inclined plates that may be fixed.Figure 5(a) presents straight and parallel flow lines and Figure 5(b) shows a parabolic velocity profile.Our contribution is the application of the Couette flow exact solution for natural fracture to obtain specific solution that describes flow profile for natural fissure.When a plate is moving; namely, system may be interpreted as a stress-sensitive system (use Couette flow).Fixed plates may be understood nonstress sensitive (use Poiseuille flow).So Couette equation is a general case with respect to Poiseuille equation.
Fundamentally, Newtonian fluid flow of a single phase through a fractured rock is governed by the Navier-Stokes equations [30,31].To simplify the discussion, we will consider a "one-dimensional" fracture and use the Navier-Stokes equations exact solution known as Couette flow, and this discussion about laminar flow between plates may be detailed in [29].Using Navier-Stokes equations, Considering the previously stated fracture flow problem regarding Figure 5(a), (1) can be simplified as follows: In Figure 5(b), sen  = / and applying in (2) Boundary conditions are  = 0 for  = 0.In the limit, when  nearly reaches the fracture aperture  and fluid velocity () is close to the terminal upper plate velocity (), then  =  for  = , and  2 / 2 = , where  = constant.
To obtain the solution of (3), it is necessary to integrate and apply boundary conditions.This solution is () = (/2) 2 +  + , which is a parabola.The constants , , and  are obtained by integration and finally result in the following equation: Equation ( 4) describes the Couette flow because there is a linear plate movement and can be used for stress-sensitive tectonic fractures.However, (4) can be written as Equation ( 5) corresponds to the steady flow pressure distribution through two inclined parallel surfaces when  = 0 (Poiseuille flow); it can be used to derive an expression for the rate flow of nonstress-sensitive tectonic fractures using  = ∫ , where  =  ⋅ 1 gives Equation ( 6) is known as The Cubic Law which has been used in fractured rocks [30], and the mean velocity () is obtained using  = /; substituting (6) gives Deriving with respect to  equation ( 5), substituting  = /2, and applying the maximum and minimum criterion, the maximum velocity can be obtained: Equations described based on the classical methods in fluid mechanics have been developed to calculate flow rate in tectonic fractures.The negative sign in ( 6), (7), and ( 8) is related to flow in the direction of the negative pressure gradient.Although the cubic equation was derived for tectonic fractures, it is used for diverse kinds of discontinuities (breccias, vugs, and channels) yet.
For natural fractures, our contribution is the application of the Couette flow general solution that describes flow profile in a natural fissure.But it will be compared with flow kinematics equations to know the range of velocity or when maximum and mean flow velocity can be used.

Kinematic Analytical Modeling for Tectonic Fractures
For flow visualization, three types of flow lines are used.Three types of fluid element trajectories are defined: streamlines, pathlines, and streaklines.They are all equivalent for steady flow but differ conceptually for unsteady flows.
Streamlines are lines tangents to the velocity vector and perpendicular at lines of constant potential, called equipotential lines, a pathline is a line traced out in time by a given fluid particle as it flows, and a streakline is a line traced out by a neutrally buoyant marker fluid which is continuously injected into a flow field at a fixed point in space [32].
To demonstrate the difference between the types of discontinuities, streamlines are used.These are associated with stream analytical function (Ψ) and potential analytical function (Φ).In tectonic fractures, their streamlines are parallel and would tend to be uniform (Figure 6).
The theory of complex variable guarantees that Laplace's equation for the velocity potential ∇ 2 Φ = 0 and for the stream function ∇ 2 Ψ = 0 is satisfied and can be solved.Then, Equations ( 9) are recognized as the Cauchy-Riemann equations for Φ(, ) and Ψ(, ), and the analytical expressions used for uniform flow in Cartesian and polar coordinates are the following: Equations ( 10) and ( 12) are Laplace's equation solutions.Thus, they are satisfied to ∇ 2 Φ = 0 and ∇ 2 Ψ = 0.For this demonstration in one-dimension is used the stream function Ψ(, ), as follows: Following a similar procedure for velocity potential Φ(, ), Thus, Laplace's equations for the velocity potential ∇ 2 Φ = 0 and for the stream function ∇ 2 Ψ = 0 have been satisfied.

Third Contradiction: Darcy's Flow or Couette General Flow for Planar Discontinuity (Tectonic Fracture)
For stress-sensitive system we recommend Couette flow, and for nonstress-sensitive system use Poiseuille flow.Initially, it has been referenced the use of Darcy's flow to describe behavior fluid flow in vugs [9,11], fault [12], fault breccias [12,13], and fractures [33].
The application of Darcy's flow is recommended when the range of value of Reynolds number (Re) is between zero and unity.Last information was reported by Muskat [34] and is applied for low laminar velocity, namely, in porous media, although it is also applied for tectonic fractures or tectonic fractured media without considering value of Reynolds number.Our third contradiction is based on calculating Re, and if this number is greater that unity, then The Cubic Law for planar discontinuity (tectonic fracture) is necessary, in which

Geological and Tomography Features of Fault Breccias
Fault breccias consist of planar discontinuities associated with faulting.These breccias are incohesive, characterized by predominant angular to subangular fragments and internal fractures of cataclastic rock present in a fault zone.Fragments are slickensided with variable slip directions, diverse sizes, and greater than 30 percent visibility.Cataclastic rocks have internal planar structure, are incohesive when faulting occurs above depths of 1 to 4 km or upper crustal fault zones, where brittle deformation increases breccia formation processes because of stresses and tectonic activity.Fault breccia zones enhance permeability [35].In effect, fault breccias are a superhighway production in limestone reservoirs.High permeability is linked with unconsolidated faulted rock that can be a channel for the hydrocarbons flow.
A fault breccia is presented in Figure 7, which shows an outcrop of limestone rock with subangular fragments, embedded in the matrix, with absence of primary cohesion and erosion; its fault breccia zone is approximately 7 meters, which implies huge movement rock mass and stresses (Figure 7).This outcrop is a point of comparison for fault breccia present in limestone reservoir, because it was clearly buried rock in the past geological time.In effect, fault breccias in limestone reservoirs are production paths with approximate 7-10 meters (width).
Limestone reservoirs show fault breccias associated with faulting; these channels have been described as horizontal, inclined, and vertical flux surfaces, which can be regarded as planar discontinuities.A sample core (10 cm in diameter) was obtained in a limestone reservoir of the Campeche Sound. Figure 8 shows a limestone core with fault breccia; its geometry corresponds to successive flux planes between barriers.These planes are connected networks in the whole medium and generate interconductivity associated with permeability and effective porosity.
On the other hand, computerized tomography (CT) was used to study the internal morphology of fault breccias.Limestone cores with fault breccias present fragments embedded in the matrix.We used information of the Integrated Ocean Drilling Program [36].Tomography images of limestone rock show the contrast between matrix and fragments that can be observed based on the different CT numbers.
Normally, fault breccia fragments are denser than the matrix, which indicates higher bulk density and low porosity, and are identified with high CT numbers.Moreover, fragments origin should be studied considering their age, lithology, CT numbers, total minerals, composition, and physical properties, to explain their density and porosity.Voids have minor CT numbers with respect to matrix and fragments.
Fault and drilling-induced breccias show high contrast between fragments and matrix because of voids that CT has identified.Fault breccias present low contrast between fragments and matrix with voids, too.These empty spaces allow the oil flow through the rock; in effect, fragments act as barriers and fluid movement is linked to voids between matrix and fragments.This can be observed in Figure 9.
Observing Figures 7, 8, and 9, the following observations can be listed: (1) Embedded clasts in a fault breccia are subangular and angular, creating a huge flow area.
(2) Fault breccias are consequence of stresses in the limestone rock, with planar and connected discontinuities that contain embedded clasts.
(4) In cores, the proportion of brecciated rock is greater than matrix.
(5) In an outcrop, the brecciated rock has dimensions of meters.
(6) Fault breccias can be described with multiple connected planes and embedded clasts.
These observations are useful for the development of an analytical model that will follow next.

Kinematic Analytical Modeling for Fault Breccias
In order to understand fluid flow in fault breccias, their kinematics should be studied.According to geological evidences and computerized tomography, these breccias show angular clasts with similar lithology, which act as barriers; they are poorly sorted and vary in size (1 mm-1 m).A representative scheme could be as shown in Figure 10 in which clasts may be grain-supported or floating and have a matrix or cement to be responsible for close or open discontinuities.
Fractures contained in a fault breccia have a size aperture of millimeters, and zone influences of fault breccias are centimeters or meters.The unfilled spaces between clasts in a brecciated zone are preferential ways, or a complex network of discontinuities for fluid flow.
Flow modeling between unfilled spaces and clasts could be developed using the flow theory near a blunt nose, or Rankine half-body.This premise is based on geometrical similarity between the angular clasts and Rankine halfbody, considering aerofoil shapes, particularly under flow conditions where the viscosity effects could be minimal [37].It is appropriate to use and to adapt the implications and equations of potential flow and streamlines to describe flow through fault breccias, considering the observations previously discussed regarding the physical phenomenon.This problem is solved in cylindrical and spherical coordinates.When cylindrical coordinates are used, physical condition is related to the radial geometry of clast and a radial function is used as a potential function.When spherical coordinates are used, the angular geometry of clast is considered and a potential function is represented using a cosine function to describe this flow problem.Then, using Laplace equation in cylindrical coordinates (, , ), The solution for ( 15) is a function as given by For our physical phenomenon, () is a radial function due to clasts changing radially, and  is an integration constant that depends on the porous media (limestone).The velocities   and   are given as Equation ( 15) is a Partial Differential Equation (PDE) and substituting ( 16) into (15) gives an Ordinary Differential Equation (ODE): Equation ( 18) is the Bessel differential equation of order zero, and its general solution is [38] where  1 and  2 are constants and   and   are the Bessel function of order zero, of the first and second kinds, respectively.A series development for the first kind of Bessel function   () can be written as 10

Journal of Petroleum Engineering
The particular solution of Laplace's equation given by ( 16) is Solution of analytical model has a constant () that is associated with material (limestone) and its roughness.
On the other hand, Laplace's equation can also be solved for a flow described in spherical coordinates (, , 0).If Φ = Φ(, ), then there is a solution Φ =   (), where () is function  = cos  and  is an integer [38].Through the previous transformation, the following Laplace equation can be expressed as follows: The velocities   and   are given as   = Φ/ and   = Φ/.Deriving the solution Φ previously stated with respect to  and substituting it in Laplace's equation given by ( 22), Equation ( 23) is the Legendre differential equation, which is a second-order ODE, and its solution for an integer  is given by the Legendre polynomials   ().Then, the solution for ( 23) is given by Legendre's polynomial of order  (0 and 1) is Equation ( 23) has a term (+1)(), which indicates a linear combination [24].A general solution has an additional term: where   and   are constants.The differential equation given by ( 22) is linear; then solutions can be superposed to produce more complex solutions: Equation ( 27) is a solution for Legendre's equation.According to Legendre's polynomial of order  with   = 1, this potential function can be used in stable steady, irrotational flow and spherical coordinates and can represent some fluid flow problems.This equation describes uniform flow, a source and a sink flow, a flow due to a doublet, a flow around a sphere, a line-distributed source, and a flow near a blunt nose.  and   play an important role in the solution because it is possible to superpose the geological events that influence the fault breccias.In addition, this solution does not depend on porous media (limestone) or experimental constant; namely, it considers fluid flow only.For example, for uniform flow (applied to planar discontinuities, whose conditions are parallel, irrotational flow), the   and   parameters in (27) simplify as Then, the potential, the stream function, and the radial and angular velocities can be expressed as For a source and sink flow (applied to discontinuity with embedded clasts, whose conditions are slow and irrotational flow), Then, the velocity potential, the stream function, and the radial and angular velocities can be expressed as Similarly, the last expressions can be deduced using (27) (the solution of Legendre's equation) and Cauchy equations.
For flow near a blunt nose or Rankine half-body, which can be modeled with the superposition of uniform flow and a source (applied to planar discontinuity with embedded clasts), as illustrated in Figure 11, Equations ( 32), ( 33), (34), and ( 35) can be used to describe the flow kinematics in fault breccias.These analytical expressions apply to a steady, irrotational, and axisymmetric flow.This problem can be considered as an irrotational flow because of the slow laminar movement of a viscous fluid on a planar surface [38].Moreover, two classical methods have been proposed and adapted to describe fluid flow through fault breccias.The first method uses a solution of the Bessel equation with an experimental constant, and second method employs a solution of the Legendre equation.

Geological and Tomography Features of Chicxulub Impact Breccias and Cantarell Reservoir
Impact breccias are a consequence of the impact of asteroids, meteorites, and comets to the Earth.Meteorites are cosmic fragments rich in iron and nickel, which have generated craters on the earth surface.Impact melt breccias and suevite can contain material from the melting of target rocks.In impact breccias, brecciated target rocks, melt fragments, and allochthonous fallback breccia can be observed.Impact breccias have been observed in cores and outcrops, with embedded fragments in the ejected matrix due to the impact.A core sequence was recovered in the Yaxcopoil-1 (Yax-1) borehole located 40 km southwest of Mérida and approximately 60 km from the center of the Chicxulub structure, between 794.65 and 894 m; a 100 m thick impact (suevite) breccia and impactites overlie a 617 m thick sequence of horizontally layered shallow-water lagoonal to subtidal Cretaceous limestone, dolomites, and anhydrites [39,40].In contrast, Grajales-Nishimura et al. [41,42], Rebolledo-Vieyra and Urrutia-Fucugauchi [43], Kring et al. [44], Wittman et al.
[45], Dressler et al. [46], Tuchscherer et al. [47], and Stöffler et al. [48] used outcrops and core sequence from the Yax-1 borehole, but there are differences with respect to lithology, age, and stratigraphic column.Most of the authors coincide closely with the mark of Chicxulub regarding, namely, its suevitic impact breccia.This impact breccia consists primarily of clasts from the underlying shallow-water lithologies, some crystalline rock of continental basement origin, and devitrified glass fragments and spherules [43,49,50].
Representative core samples of Yaxcopoil-1 were analyzed for the estimation of hydraulic permeability.For Tertiary limestones and suevites, their bulk permeabilities range between 10-14 and 10-19 [m 2 ] and their porosities are between 0.08 and 0.35.For Lowermost Suevite Unit, Cretaceous anhydrites, and dolomites (900-1350 m), their permeabilities range between 10 −15 and 10 −23 [m 2 ] and their porosities are between 0.1 and 0.15.The previous permeability and porosity values do not include macroscopic events, such as faults and tectonic fractures [17].
Three decades ago, the impact crater discussion was begun.In 1975, Lopez-Ramos [51,52] and Penfield and Camargo in 1981 [53] reported 210 km diameter circular structure with total magnetic-field data.In 1991, Hildebrand et al. [54] suggested that a buried 180 km diameter circular structure (the Chicxulub impact crater) was located on the Yucatán Peninsula, Mexico, which was formed 65 million years ago [46], proposing it as candidate for the K/Pg boundary impact site.
In the offshore zone of the western margin of Yucatán Platform or Campeche Bay is located the largest oil field in Mexico and has produced more than 18000 million barrels of oil and 10000 billion of cubic feet of gas [55].Outcrops analogs, petrographic analysis, well logs, and cores description indicate that Cantarell Oil Field is genetically related to the Chicxulub meteorite-impact [4].This genetic link is found for the Cretaceous-Tertiary (K/T) in Cantarell oil field that can also been seen in the Guayal (Tabasco Region) and Chiapas outcrops [41].
Impact breccia clasts are impact melt fragments with ripup morphology that are consolidated and embedded in the matrix and can be observed in Figure 12, showing similar geometrical characteristics.Figure 12(a) shows a Campeche Sound core with embedded clasts and rip-up morphology, and Figure 12(b) shows an analog outcrop with rip-up morphology clasts, too.Considering Figure 12, it can be inferred that embedded clasts act as nonflow barriers.Moreover, impact breccias contain ejected material, which generate nonplanar discontinuities, which can be observed in cores and in outcrop samples (Figure 12).
The Cantarell reservoir includes the Upper Cretaceous (Upper Breccia) that is associated with the Chicxulub event, with total average porosity and permeabilities ranging from 8 to 10% and 800 to 5000 md, respectively, with average thickness of 11 to 105 meters, thinning to the south-west, showing dissolution and dolomitization and the Upper Jurassic (Tithonian) with dolomitized limestone.The field oil production is associated with tectonic fractures that provide high permeability [4,8,56,57].Computerized tomography is a tool to observe the differences between matrix and ejected clasts and describe internal texture of rock. Figure 13 indicates that clasts are nonflow barriers in impact breccias and generate nonplanar discontinuities.Figure 13(a) shows different textures in breccia sequence of Yax-1 well; clasts or nonflow barriers are present; it indicates that there is a range of porosities and permeabilities in limestone rocks.Low porosity and permeability are related to fine ejected material.Figure 13(b) shows other impact breccias (Bosumtwi) in three dimensions; the nonflow barriers (high-density clasts) can be observed, with rip-up geometry, generating nonplanar discontinuities, and they are embedded in matrix rock (low-density).
Observing Figures 12 and 13, the following can be inferred: (1) Embedded clasts in an impact breccia have rip-up geometry, creating a nonflow barrier in the porous media (matrix).
(2) Impact breccias are a consequence of the impact of asteroids, meteorites, and comets in the Earth, generating nonplanar discontinuities that contain embedded clasts.
(3) Porosity and permeability in an impact breccia are associated with ejected material (clasts), providing a type of primary porosity in the limestone.
(4) In the core, the proportion of matrix rock is greater than embedded clasts.
(5) In the outcrop and reservoirs, brecciated rock dimensions are related to impact size.
(6) Impact breccias can be described with multiple embedded clasts (high-density) that act as nonflow barriers into connected matrix (low-density).
The previous inferences will be used in the present paper for the development of an analytical model.

Kinematic Analytical Modeling for Impact Breccias
When an impact breccia is observed without the presence of other geological events as fault breccias, tectonic fractures, vugs, sedimentary breccias, dissolution, and dolomitization, its permeability would be ultralow [17] and can have a low to moderate porosity (1-8%) [4].In consequence, impact breccia can act as seal and have fluid storage, but its oil production depends on tectonic fractures, fault breccias, dissolution, and vugs.It is very highly observed in the Cantarell field.Because of its low permeability and porosity, fluid flow velocity in impact breccias is slow and can be considered as irrotational flow in a porous media with primary porosity.In addition to low permeability, rip-up geometry observed in tomography and geological evidences and clasts act as barriers for fluid flow in porous media.Figure 14 illustrates fluid flow with nonflow barriers and rip-up geometry for an impact breccia.
The impact breccia flow may be studied in terms of the streamlines, taking into consideration the axial flow symmetry.To solve this problem, we apply and adapt the flow through a sphere considered by Stokes; then, it is possible to use the Laplace Biharmonic equation to describe this flow [59,60]: For spherical coordinates, (36) can be expressed by The boundary conditions for (37) are given by The boundary conditions given by (38a) and (38b) describe fluid adherence on a clast with rip-up morphology in a porous media.Figure 15 represents this fluid adherence at a clast with similar rip-up morphology.Boundary condition given by (38c) describes the streamlines before-after a breccia clast, which implies a velocity change in fluid flow because the clast is a barrier; namely, for  → ∞, the final clast velocity is obtained.In accordance with this boundary condition is a general solution for Laplace Biharmonic equation expressed as follows: This solution proposes radius and angle change of clast.However, it is necessary to study the Stokes solution and test if it can be adapted to oil flow in an impact breccia.Equation (39) contains (), which is a radial function related to the clast radius.Substituting (39) in ( 37), Considering streamlines with trajectory angle of 90 ∘ , then (40) can be expressed as Equation ( 41) is a fourth-order differential homogeneous linear equation, and its solution is of type () =   , where  can have values (−1, 1, 2, 3) and  includes integration constants (, , , ), such that Deriving (39) with respect to angle , Ψ/ = 2() sin  cos , and substituting it in (38a), Substituting ( 42) in (43), can be expressed as Substituting ( 42) in (39) and deriving the resulting equation, Then, substituting (46) in (45), Applying boundary conditions given by (38a) and (38b) to (44) If  → ∞, clasts should be smaller in layer thickness (two orders of magnitude).Moreover, initial fluid velocity changes when fluid impacts with clast (barrier), but fluid velocity must be different from zero to apply the mass and momentum conservation principles.Thus  = 0 because  ∈ R and  ∈ R. From (50), If   =   = 0 ((38a) and (38b)), then, the previously described procedure regarding velocity   indicates a stagnation point.Following a similar solution for   , the following equation can be derived: Substituting  = /2 and  = 0 and  = −90 ∘ (perpendicular streamline, that describe streamline around clast) in ( 52), For   from (44), Substituting  = /2 and  = 0 and  = 0 ∘ (radius localized at  = 0 ∘ ) in ( 54), Solving the system of ( 53) and ( 55), constants  and  are expressed as follows: Then through the expressions (values) for the parameters , , , and , (44), (47), and (39) can be written as The potential velocity is obtained using   = 1/ (Φ/) and integration.Then, Φ can be derived as As impact breccias clasts do not have spherical geometry and their rip-up morphology shows subrounded sides and other angular sides, in our analytical model we consider symmetry and an angle between 0 ∘ and 90 ∘ .Moreover, we propose that the range of  may be 0 ∘ <  < 180 ∘ and rate change with respect to the angle can be expressed as Equations ( 60), (61), and ( 62) describe the changes in streamline velocities in impact breccias and could be used to study clasts with a variety of angles or subrounded side.In effect, as stated, the Stokes solution was adapted to impact breccias.

Fourth Contradiction: Fluid Flow of the Cantarell Reservoir Modeled without Considerer Chicxulub Impact
Several authors document the influence of the Chicxulub impact in the Campeche Sound and discussed if the Chicxulub impact breccias are seal, production zone [4,41,42], or the impact is a geophysical survey reference as part of an oil exploration program [16,53,58,61].
On the other hand, the Cantarell field is modeled as a tectonic fractures reservoir [6-8, 56, 62].As this paper is focused on oil flow in limestone reservoirs, then there is a question: why is the Cantarell reservoir modeled with tectonic fractures without considering the fluid flow through impact breccias?Or the impact breccia oil does not flow?A contribution of this paper is related to describing fluid flow through an impact breccia.In absence of vugs, fractures, and fault breccias, impact breccia has moderate porosity and low permeability, which indicates low flow velocity, unique in this type of breccia.Clearly, we modeled impact breccia without fractures or other geological events.

Geological and Tomography Features of Sedimentary Breccias
Sedimentary breccias are a type of clastic sedimentary rock with subangular to subrounded clasts.These rocks are deposited by transport and fast moving of a body of sediment particles by water and/or air.These breccias are observed in debris flows, mud flows, and mass flows, such as landslides and talus.
According to type of clasts, sedimentary breccias can be monomict, oligomict, or polymict.Clasts may be extraformational or intraformational, matrix-supported or clastsupported.The morphology of fragments clasts has been reworked by transport.Moreover, rocks with rounded clasts are known as conglomerates, and they are different from breccias.
Sedimentary breccias contain clasts embedded in the matrix.These breccias do not have linear or internal planar structure resulting from lithification and consolidation of rocks, and they have been seen in outcrops (Figure 16). Figure 16 shows reworked clasts embedded in rock matrix with subrounded sides which indicates a short clasts transport.Long transport implies that clasts are rounded and can be modeled as ellipsoids.
In principle, sedimentary breccias affect carbonate reservoirs and may enhance the total porosity.Fluid flow can be significant in the matrix-clast interface due to clast being impermeable and acting as flow barriers.Figure 17 shows a sedimentary breccia core with embedded clasts in a limestone matrix, with subrounded and subangular side.
Brecciation often results in enhanced porosity, but with poor interconnection of pores.A specific example is the Cretaceous Debris Reservoir, Poza Rica Field, VER, Mexico [15].On the other hand, we used information of the Integrated Ocean Drilling Program that had sedimentary breccia tomography images [36].Denser fragments embedded had a contrast with the matrix, indicating high bulk density and low porosity.In many cases, clasts can have high density and ultralow porosity.
In Figure 18, tomography shows dark shading that corresponds to low density and white to high density regions; this figure presents poorly sorted, elongated, and impermeable clasts with low porosity, in addition to Debris flow sediments in different layers.These events indicate that clasts are oligomict or polymict and extraformational that act as flow barriers in limestone reservoirs.
Observing Figures 17 and 18, the following can be inferred: (1) Embedded clasts in a sedimentary breccia have subrounded, reworked, and elongated geometry, creating a nonflow barrier in porous media (matrix).
(3) Sedimentary breccias are consequence of Debris flow, generating nonplanar discontinuities that contain embedded clasts.
(4) Porosity and permeability in a sedimentary breccia are associated with matrix, embedded clast, and interface matrix-clast, producing a type of primary porosity in the limestone rock.(5) In the core, matrix rock portion is greater than embedded clasts.
(6) The presence of Debris flow sediment in different layers indicates that there are diverse processes of sedimentation and that rock was exposed in the surface; namely, it was an outcrop in another geological age at subsurface conditions.(7) In outcrops and reservoirs, sedimentary breccia dimensions are related to Debris flow.
These observations are stated with the objective of the development of an analytical model.

Kinematic Analytical Modeling for Sedimentary Breccia
Fluid kinematics could describe fluid flow in this type of rocks.A representative scheme is shown in Figure 19, with streamlines for sedimentary breccia clasts, where reworked clasts may be grain-supported or matrix-supported.
Modeling between matrix-clast interfaces could be developed using flow around an elliptical body, like the Rankine solids, due to reworked clast.The Rankine methodology is similar to that previously applied to fault breccias to describe fluid kinematics; therefore, the flow around an elliptical body should satisfy Laplace's equation [64].
Considering uniform flow, the Rankine solution is obtained using the superposition of a linear source and a linear sink of equal strength, combined with uniform flow (Figure 19), given by Conceptually, ( 63) and ( 64) express that a sink and a source are equivalent, but geometry shows differences: they have different angles with respect to streamlines (Ψ) and are separated from distance .Distance between origin and source is /2 due to their symmetry.This is observed in geometry indicates angular rate, and oil flows around the impermeable clast, generating a nonplanar discontinuity: where Considering ( 66) and ( 67), the combined flow (Ψ  ) is given by Figure 21 shows two stagnation points associated with a source and sink; the length of the Rankine body is (  ),   =  1 +  2 , and the width of the Rankine body () is related to the maximum value of  on the contour of the body in the -axis direction.
Defining   = Ψ  / and   = −Ψ  / in rectangular coordinates, using (66), horizontal and vertical velocities are given by The total velocity is given by  = √  2 +   2 , and potential velocity is Equation (70a) can also be expressed substituting (67) for  1 and  2 : To determine the width of the Rankine body, the maximum value of  ( max ) in Figure 21 should be estimated at  = 0 (V  max = 0).Geological information could provide the width and length of clasts embedded in cores of a sedimentary breccia, which would be assumed as length and width of the Rankine body.

Geological and Tomography Features of Sedimentary Breccias
Vugs are a type of pores in carbonate rocks.Choquette and Pray (1970) [65] presented a classification based on the pore space genesis.Lucia [66] proposed a classification focused on petrophysical properties and on distribution of pore sizes within the rock.Vuggy pore space is classified into touching-vug pores and separate-vug pores.Touching-vug pores as tectonic fractures, breccias, and caverns are nonfabric-selective in their origin.Separate-vug pores are defined as pore space larger than the particle size and fabric-selective in their origin and are interconnected only through interparticle pore space, such as moldic, intrafossil, intragrain, and shelter pores [66].
According to Lucia [66], vugs compared to separatevug pores are secondary solution pores that are not fabricselective in their origin, with irregular sizes and shapes, which could be interconnected [65].
In absence of tectonic fractures, vugs are nonfabricselective pore spaces or discontinuities from various sizes, with irregular shape, as a result of chemical dissolution, that could be interconnected or separated.Figure 22 shows a core with vugs created by chemical dissolution and irregular shapes.Normally, in a double porosity reservoir, vugs interact with the matrix.
Permeability of vugular porous media depends on its interconnection of the pore space.In this study, vugs are considered as nonplanar discontinuities that may be separated and nonfabric-selective and interact with the matrix of carbonate rocks.
Vuggy porosity can be produced by the interaction of meteoric waters with rock; by grains dissolution, fossils, and matrix pores; by deep-burial fluids; and by exposition of rock surfaces in humid climates.
In limestone outcrops, vugs are interconnected only with matrix; vuggy pore space also can be regarded as interconnected with matrix and other vuggy systems (Figure 23). Figure 23 shows a chemical dissolution process (diagenesis) with multiple size vugs similar to Figure 22; then core and outcrops could be used as analogs.
Tomography images of an oil impregnated core obtained from a vuggy limestone reservoir were taken to determine the internal characteristics, geometry, and connectivity of the pore space.
The obtained one hundred and thirteen images describe the geometry and irregular shapes of the vugs.Figure 24 shows an oil saturated core, with a length of 36 cm, with visible and irregular vugs as result of a chemical diagenesis, and CT slices were taken every 3 mm of distance through the sample (Figure 25).
CT slices for the vuggy core of Figure 24 are presented in Figure 25. Figure 25 shows slices with vugs (black color), matrix (yellow color), filled or recrystallized cavities (green color), and embedded clasts (orange color).Cavities with black color are big pore spaces, or void spaces saturated with oil.When the slices are compared, vugs connectivity changes are observed.If we see a vug with black color in an image, it disappears in subsequent images.In contrast, connected vugs can be observed through different slices.
Considering core axisymmetry, there are permeable zones with isolated porosity and others with interconnected porosity.Isolated zones interchange fluid with matrix, but   connected region presents matrix-vug and vugs system flow.Moreover, dissolution is the dominant process that enhances connection vugs.

Kinematic Analytical Modeling for Vugs
Vugs are irregular spaces, like spheroids and ellipsoids, that interact with the matrix.
To predict the flow field within a vug, we consider the mathematical solution developed by Neale and Nader [9], which was used to describe the pressure distribution within an isotropic and homogeneous porous medium, with a spherical cavity.
The problem and solution described by Neale and Nader [9] are represented as shown in Figures 26 and 27.
To develop the analytical solution for the flow problem described, a composite porous medium (matrix and vugs) was considered, and the creeping Navier-Stokes and the Darcy equations were used.Combining these two equations, the Laplace Biharmonic equation in spherical coordinate can be obtained and solved with its boundary conditions; the solution is given by using the normalized radial coordinate and the normalized radius of the cavity, where Equations ( 71a) and (71b) with their constants (, , , , , , ℵ) predict the streamlines behavior in vugs and porous media.However, the authors Neale and Nadar did not describe angular, radial velocities or potential function.Considering (71a) and (71b) with their constants, we developed the angular, radial velocities and potential function, which are given by Defining potential flow as Φ = ∫ r 0   , then Equations ( 73) and (74) provide a description of the fluid flow in a composite porous media.

Fifth Contradiction: Liquids Retention Paradox for Vugs
well-known in fluidized bed reactors and their design in the chemical engineering area [60].
During the production, oil entrapping is a result of two fluid dynamic processes.Initially, oil can be saturating the cavities, and its flow obeys a pressure gradient and a balance between viscous and inertial forces; thus, this is a dynamic flow process.When the first process concludes, oil flow follows a diffusion process, with slower velocities, generating a second process with static characteristics and fluid entrapping.The described phenomenon is related to the cavity geometry and vuggy pore space; this problem is associated with static-dynamic liquid retention in vugs and should be called liquids retention paradox.
It is a paradox because liquid in the cavities may be trapped in stagnation or dead zones; however, fluid flow will always occur in the vuggy porosity (secondary) producing a change in fluid velocity.In consequence, stagnation zones do not exist because there is always movement of fluid.

Application Examples and Discussion
In this section, examples are presented to illustrate the proposed kinematics models in the analysis of limestone reservoirs with different discontinuities.

Behavior of Fluid Velocities.
To examine the differences between the types of discontinuities, we used analytical kinematics models to calculate and compare fluid velocities.Key data and parameter values employed are presented in Table 1.Note that quantitative models should be implemented with geological characterization for a better prediction.
Additionally, to quantify fluid velocities in this study, we took into consideration a pressure drop, a discontinuity angle with respect to streamline, aperture, clasts angle, and sink and source for sedimentary breccia which are included in Table 2.
Table 3 shows obtained velocities and flow rates values for diverse kinds of discontinuities.The goal is to compare these parameters.
These models and their results would be applicable to oil limestone reservoirs.In this example, the calculated values of Table 3 illustrate that discontinuities related to fluid flow are dissimilar and that flow velocities and volumetric flow can be different.The results suggest that discontinuities in a limestone reservoir may not yield the same level of oil production.This implies that it is necessary to identify and verify the dominant discontinuity using static and dynamic characterization.
Equation ( 12) (The Cubic Law) is used for tectonic fracture.Calculated values show the differences for each type of discontinuity.Horizontal tectonic fracture presents equivalent velocities (radial and angular) because streamlines are  The positive sign in fluid velocities represents its direction from of origin in -axis or -axis direction.
parallel and volumetric flow depends on fracture aperture.Fault breccia has high velocity and flow because it depends on elongated and subangular clast.Also, vugs have superhigh velocity and flow because they are connected cavities without flow barriers.In contrast, impact and sedimentary breccias have low velocity and volumetric flow due to clast geometry (rip-up and ellipsoidal), creating low radial velocity.In addition, clasts are flow barriers and fluid flow is related to interaction matrix-clast and their permeabilities that are normally very low because they behave as primary permeability and porosity.This implies that proportion matrix is greater than the proportion clast.

Carbonate Reservoir Characteristics: Cardenas Field
Application.According to the proposed geological model for the Cardenas Field, which is an anticline with a fault system, controlled by stratigraphic traps rather than structural traps.In accordance with the characteristics of primary porosity (1.5-1.7%),secondary porosity (5-11%), primary permeability (395 × 10 −17 -329 × 10 −17 m 2 ), and secondary permeability (196 × 10 −10 -89 × 10 −10 ), production layers are subdivided into different breccias intervals in the reservoir.Figure 28(a) shows correlation through longitudinal breccias intervals for the Cardenas Field wells; moreover, breccias sequence with a chemical diagenesis (dolomites and vugs) and limestone layers were a criterion for the selection of wells.It is shown in the cross section (Figure 28(a)).Oil production in the Cardenas reservoir comes from lateral interconnected sedimentary breccias and vugs [63].
Interconnected vugs by microfractures provide high permeability and porosity.On the other hand, sedimentary breccias are related to salt tectonics and generate permeability and porosity, which are low.
Application of the flow kinematic models (vugs, sedimentary breccia models) to the Cardenas Field may be used as a diagnosis tool for the fluid velocities prediction and in the discontinuity characterization.In this case, there are wells with a similar pressure behavior (Figure 28(b)) that produce from breccia intervals with vugs and microfractures.
In Figures 28(a) and 28(b), we observe a cross section 1-1  with eight wells and their similar pressure history.In Figure 28(b), 3D seismic section shows wells layers correlation and confirms their lateral continuity.Also, the section correlation shows clearly a connected thickness of Debris flow.As an application, we have chosen three wells: Cardenas-109, Cardenas-129, and Cardenas-308, with geologic heterogeneities related to sedimentary breccias and vugs.The goal is to compare fluid velocities and characteristics flow in sedimentary breccias and vugs and validate them with production data.Wells data and parameters are given in Tables 4, 5, and 6.
Figure 29 shows wells Cardenas-109, Cardenas-129, and Cardenas 308.In accordance with this figure, well Cardenas-109 has a high oil cumulative production (>20 MMBls) due to geological events juxtaposition, of sedimentary breccias and vugs.Vugular porosity was created by dissolution processes associated with pressure and temperature drop, and circulation of high corrosive hydrothermal fluids, and its breccias deposits were exposed to subaerial erosion, after they affected by dissolution and dolomitization.In addition, oil cumulative production for well Cardenas-129 is associated with vugular porosity, although its described production (12-20 MMBls) in Figure 29 is smaller than for Cardenas-109.Well Cardenas-308 geological description is related to sedimentary breccia intervals.Its production (6-12 MMBls) is low with respect to the other two wells (Figure 29), but it can be considered that there is a high oil storage in the breccia intervals; namely, oil storage is localized in its primary porosity.
According to Figures 28(a), 28(b), and 29, diverse volumetric flow and velocities should be calculated because geological events for the Cardenas Field are different and juxtaposed.Juxtaposed geological events imply a high volumetric flow and fluid velocity.
We applied the kinematic equations as a semiquantitative diagnosis tool to determinate fluid velocities and fluid flow for the three study wells.Used equations for sedimentary breccia, vugs, and superposed flow (sedimentary breccia and vugs) are ( 57), (58), and (69) with their geometrical parameters.The fluid velocities can help in the interpretation of the type of geological discontinuity; moreover it is necessary    store oil.This indicates that path for fluid flow will be slow and tortuous, and then its velocity and flow are low.This premise was corroborated in Table 7.
Well Cardenas-129 is studied considering a high oil velocity because vugs are cavities that do not have flow barrier and they are normally connected with limestone matrix.The porosity in this reservoir layer is a preferential way for fluid flow.When there are connected vugs with oil storage in the rock, production conditions are excellent due to oil free path.In contrast, well Cardenas-129 oil production should be greater than well Cardenas-308 oil production.This claim was demonstrated in Table 7. Well Cardenas-109 presents complex geological events juxtaposition.Sedimentary breccias store fluid and vugs store and transport oil through porous media due to their interconnection, in this case, porosity (storage) and secondary permeability (transport).Then, the already stated events juxtaposition is applied to the geological events superposition, which is a characteristic of analytical models developed in this paper because our equations are lineal and they are solutions of the Laplace equation.In consequence, the maximum oil production in this well must be obtained, and this is demonstrated in Table 7.
Table 7 shows the obtained results with input parameters of wells Cardenas-109, Cardenas-129, and Cardenas-308.Chosen wells for models validation have diverse production; in consequence, fluid velocity and flow volumetric are different and they are related to geological event as vugs, sedimentary breccia, and their juxtaposition.
The wells production is associated with phenomenology of complex limestone reservoir.The key point here is that sedimentary breccias contain embedded clasts that are barriers for fluid flow; nevertheless they can store oil.The degree of complexity of clasts interactions with fluid depends on clast diameters and their spacing.In the case of vugs, it depends on their interconnection.When different geological events are juxtaposed and interconnected as in well Cardenas-109, oil production is high.
The Cardenas Field has been chosen because we know, understand, and have geological control of reservoir, which was developed in a doctoral thesis.Data for the field application previously discussed was mainly borrowed from the Ph.D. dissertation of one of the authors of the present paper [63].

Conclusions
This paper has presented a discussion on the effects of planar and nonplanar discontinuities in fluid flow of carbonate reservoirs.The proposed mathematical models have provided a simple method to identify the flow characteristics of fault breccias, vugs, tectonic fractures, impact breccias, and sedimentary breccias.
From the results of the study the conclusions are as follows: (1) It has been demonstrated through the analytical models of this paper, tomography, and observations of outcrops and cores that planar and nonplanar discontinuities in limestone reservoir show distinctive representative flow characteristics.Fault breccias, tectonics fractures, sedimentary breccias, vugs, and impact breccias have flow and geological patterns that affect oil production and generate differences in static and dynamic characteristics that affect flow behavior.
(2) It was demonstrative that discontinuities (vugs, fault breccia, and tectonic fractures) are superhigh production ways and others (impact and sedimentary breccias) are storage zone.In effect, each discontinuity is different and cannot be called or analyzed as tectonic fractures, unknowing geological genesis and flow patterns.
(3) It has been shown that the unknowing of the origin of discontinuities causes confusions and contradictions for static-dynamic characterization, simulation, and exploitation strategy of limestone reservoirs.This is one of the most important conclusions of this study.
(4) Fluid velocity and volumetric flow for vugs (nonplanar discontinuity) are the greatest obtained values between all of discontinuities because they are cavities that do not have barrier flow.In consequence, a limestone reservoir well with connected vugs will have a high oil production.
(5) In the absence of fault breccias, vugs, sedimentary breccias, and tectonic fractures, impact breccias in limestone reservoirs present low fluid velocity and volumetric flow because they have flow barriers (ejected clasts) that act as a type of primary porosity depending on their values of porosity and low permeability.In effect, these impact breccias would not have high oil production.In this case, impact breccias must be superposed with an additional geological event for a high volumetric flow.
(6) Oil production of limestone reservoirs with juxtaposition of geological events (breccias, fractures, and vugs) is high, obeying a dominant geological event in fluid production (vugs, fault breccias, and tectonic fractures) and other dominant geological events in fluid storage (sedimentary breccias and impact breccia).

Figure 4 :
Figure 4: Scan of limestone core, with dark shading associated with low density and white with high density regions with evident macroscopic fractures, Early Cretaceous, Samaria-Luna Region, TAB, Mexico.

Figure 5 :
Figure 5: Top and lateral views of fluid flow in a tectonic fracture.(a) Flow lines in an open fracture (top view).(b) Flow profile between two inclined parallel plates (lateral view) [29].

Figure 6 :
Figure 6: Streamlines in uniform rectilinear flow in tectonic fractures with parallel walls, (a) horizontal direction, (b) inclined direction, and (c) vertical direction.

Figure 9 :
Figure 9: Representative appearance of fault breccia in computerized tomography (CT) images.(a) Slice of fault breccia cut perpendicular to core axis (interval 316-C0004D-28R-2, length: 45.13 cm).(b) Same as (a) with CT numbers of matrix and fragment.Note the small contrast in CT numbers between matrix and fragment (1162 versus 1294) [36].

Figure 13 :
Figure 13: Samples of CT slices through the Bosumtwi and Chicxulub impacts.(a) Breccia sequence in Yaxcopoil-1 borehole.Core images obtained with the Core Scan System for the breccias sections, illustrating the different textures [61].(b) Three-dimensional reconstruction of the clasts population of the Bosumtwi suevite.The image on the left shows the sample with color-coded clasts, with the matrix (groundmass) rendered partly transparent.The center image shows only the high-density clasts, with the low-density clasts (and the groundmass) rendered transparent, and the image on the right shows only the rock edges and the low-density clasts [27].

Figure 14 :Figure 15 :
Figure 14: Streamlines for the flow through in rip-up fragments in an impact breccia.

Figure 21 :
Figure 21: Streamlines for the Rankine solid with sink and source [64].

Figure 24 :
Figure 24: Oil impregnated core of a vuggy limestone reservoir.Late Cretaceous, the Campeche Sound, Marine Region, Mexico.

Figure 26 :Figure 27 :
Figure 26: Lateral view of a composite porous medium with vugs, matrix, and a fluid flow field.

Table 2 :
Additional data for the calculation of fluid velocities.

Table 4 :
Data and parameters for the Cardenas Field, well Cardenas-109.

Table 5 :
Data and parameters for the Cardenas Field, well Cardenas-129.

Table 6 :
Data and parameters for the Cardenas Field, well Cardenas-308.considerer static characterization for estimate values of velocities and rates with reduced uncertainty.Kinematics equations validation is developed considering a low fluid velocity in the sedimentary breccia of well Cardenas-308 because clasts are barriers during fluid flow, but porosity and permeability of the sedimentary breccia may to

Table 7 :
Calculated velocities and rates values for the Cardenas Field, wells C-109, C-129, and C-308.