Impact of Diagenetic Alterations on the Petrophysical and Multiphase Flow Properties of Carbonate Rocks Using a Reactive Pore Network Modeling Approach.

on the and Flow of Résumé — Impact des altérations diagénétiques sur les propriétés pétrophysiques et d’écoulement polyphasique de roches carbonates en utilisant une modélisation par l’approche réseau de pores — Les roches sédimentaires présentent souvent une structure porale hétérogène qui est intrinsèquement liée à la texture de la roche d’origine et aux modifications diagénétiques subies. Ces altérations sont régies par la texture de la roche d’origine, les fluides impliqués (et les interactions rock/fluide), l’histoire de l’écoulement et les conditions physico-chimiques. Tout au long de la paragenèse, l’altération de la dissolution (amélioration de la porosité) et de la précipitation (réduction de la porosité) causée par les changements de conditions chimiques et thermodynamiques peut conduire à des roches hétérogènes à la fois à l’échelle locale et à l’échelle du réservoir. En l’absence d’échantillons de roche de réservoirs permettant de mesurer les propriétés pétrophysiques ( i.e. porosité, perméabilité et facteur de formation) et les propriétés polyphasique ( i.e. pression capillaire, perméabilité relative et indice de résistivité), Abstract — Impact of Diagenetic Alterations on the Petrophysical and Multiphase Flow Properties of Carbonate Rocks Using a Reactive Pore Network Modeling Approach — Sedimentary reservoir rocks generally have complex and heterogeneous pore networks that are related to the original depositional rock texture and subsequent diagenetic alterations. Such alterations are in part controlled by the original mineralogy and sedimentological facies, the compaction history, the involved fluids (and rock/fluid interactions), the flow history and the related physico-chemical conditions. During the diagenetic evolution (paragenesis), cycles of alternating dissolution (porosity enhancement) and precipitation (porosity destruction) caused by changes in chemical and thermodynamic conditions may lead to heterogeneous rock structure at both local and reservoir scale. In the of cored The correlation between these numbers and the dissolved/precipitated layer thickness distribution is investigated. This work contributes to improve the understanding of the impact of dissolution and precipitation on permeability and porosity modification. Using the PNM approach, multiphase flow properties and permeability-porosity relationship have been determined for different reactive flow regimes. These relationships are relevant input data to improve the quality of reservoir simulation predictions.

, the flow history and the related physico-chemical conditions. During the diagenetic evolution (paragenesis), cycles of alternating dissolution (porosity enhancement) and precipitation (porosity destruction) caused by changes in chemical and thermodynamic conditions may lead to heterogeneous rock structure at both local and reservoir scale.
In the absence of cored plugs to measure the petrophysical properties (i.e. porosity, permeability and formation factor) and multiphase flow properties (i.e. capillary pressure, relative permeability and resistivity index), a numerical tool that calculates these properties from pore structure data by predicting its evolution during the diagenetic cycle is of great interest for the petroleum industry and reservoir characterization studies. A Pore Network Model (PNM) provides opportunities to study transport phenomena in fundamental ways because detailed information is available at the pore scale. It has been used over the last decades to understand basic phenomena such as capillarity, multiphase flow or coupled phenomena. In particular, this modeling approach is appropriate to study the rock/fluid interactions since the mass exchange at surfaces can be modeled explicitly. It can provide quantitative information both on the effective transport property modifications due to the reactions and on the structure evolution resulting from dissolution/precipitation mechanisms. In the present paper, this approach is used to study the effect of the diagenetic cycle on the petrophysical properties of carbonate rocks. It involves three discrete steps. The first step consists of replacing the original complex pore structure of real porous media by a conceptual network. The second step consists of resolving the governing equations of the precipitation and dissolution phenomena (i.e. reactive convection diffusion equation) in the conceptual 3D pore network and deducing the local reactive fluxes and the motion of the fluid-solid interface. The third step consists of updating the new pore structure and calculating the new petrophysical properties of the modified porous media. Those steps are repeated in order to mimic a given diagenetic scenario. Finally, the multiphase flow properties of the current porous media are calculated. The impact of one diagenetic cycle of dissolution and precipitation on the pore networks' heterogeneity and consequently on the petrophysical properties (i.e. porosity and permeability) and multiphase flow properties (i.e. relative permeability and capillary pressure) have been investigated. The permeability and porosity evolution during a given diagenetic cycle are calculated and analyzed as a function of the relevant dimensionless numbers (Peclet and Damköhler numbers) that characterize the flow and reaction regime. The correlation between these numbers and the dissolved/precipitated layer thickness distribution is investigated. This work contributes to improve the understanding of the impact of dissolution and precipitation on permeability and porosity modification. Using the PNM approach, multiphase flow properties and permeability-porosity relationship have been determined for different reactive flow regimes. These relationships are relevant input data to improve the quality of reservoir simulation predictions.

INTRODUCTION
Carbonate reservoir rocks generally are characterized by complex and heterogeneous pore networks which are the result of depositional processes and conditions (climate, hydrodynamics, tectonic, etc.) and diagenetic alteration immediately following deposition, through burial and even continuing today (e.g. Choquette and Pray, 1970;Moore, 1989Moore, , 2001Lucia, 1995;Lonoy, 2006). Primary depositional textures (grainstone, packstone, wackestone, etc.) have relatively simple relationships with the hydrodynamic regime and resulting pore network. For instance, a grainstone (Dunham, 1962) is a grain-supported texture lacking a mud matrix (Fig. 1a). The intergranular pores will, therefore, have dimensions that are inherently related to the packing, sorting and size of grains. Grain size has a positive relationship with pore size and permeability (e.g. Lucia, 1995;Lonoy, 2006). On the other hand, packstone, wackestone and mudstone (Dunham, 1962) have increasing mud matrix from grain-supported to mud-supported (Fig. 1b). Carbonate mud is made up of microcrystalline calcite crystals or rhomboids (size ranging from 1-8 μm), which results in a considerable intercrystalline porosity, yet the average size of the pores is relatively in the order of a few microns. This results in very low permeability, although the bulk porosity could be higher than that of coarse grainstone.
Diagenetic processes include dissolution, precipitation, mineral replacement, fracturing and compaction (Moore, 2001). For the sake of simplicity here, only the first two processes are investigated. Dissolution represents a porosity increase due to leaching of grains, mud matrix and/or cement. If the grains are selectively leached, the resulting pores are called molds. If the molds are further dissolved but the original shape of the leached grain is still somehow recognizable, the mold is called solution-enhanced mold (Fig. 1c). With increasing dissolution, the original grain cannot be recognized and the resulting pore is a vug. At any stage during diagenesis, cement may fill the empty pore. Cements can completely occlude the pores or partially line or fill the pores (Fig. 1d).
Dissolution and precipitation processes affect the sediment rocks throughout their geological evolution from deposition to burial and possible uplift. As a consequence, different fluids and extrinsic conditions (pressure, temperature) affect the mineralogically instable rock throughout its history. The porous structure is subjected to a dynamic evolution which has been already discussed (Fig. 2) by Choquette and Pray (1970) for dissolution of grains and precipitation or filling within molds and vugs. Today's petrographic techniques allow to identify the various phases of dissolution and precipitation and to order them chronologically (i.e. propose a sequence of diagenetic phases, a paragenetic evolution). Yet, more work is needed to quantify properly each recognized dissolution and precipitation phase as well as their particular impact on the petrophysical properties of the porous media. When vuggy textures are observed, the first important criterion is related to whether these vugs are inter-connected or not (cf. Lucia, 1995;Lonoy, 2006). Vugs that are isolated and not connected do not improve permeability; inter-connected vugs may result in large and small-scale wormholing and subsequently may increase tremendously the permeability.
Mineral replacement, such as dolomitisation, changes the original mineralogy of the carbonate rocks. The produced crystals induce also a change in the porous space (between the precipitated/replaced crystals). Dolomite rhombs, usually, induce about 13% of increase in porosity with respect to the pre-existing calcite (if mole/mole replacement). However, one must not forget the effects of the crystal-size. If a grainstone (with coarse grains) is dolomitized, the produced dolomite crystals are usually also coarse, but the gain in porosity is minimal (if not negative). However, a dolomitized mudstone will definitely possess higher porosity, especially if the dolomite rhombs are relatively coarse. In all cases, due to the dolomite crystal fabrics and the fact that dolostones get easily fractured compared to original limestones, the resulting permeability is usually enhanced (Lucia, 1995).
Dolomitisation has always been a matter of debate when it comes to reservoir potential enhancement or destruction (compared to the original rock, before dolomitisation). Lucia (1995) considered the types of dolomite crystalline texture and the crystal sizes, as well as the original texture of limestones, in order to present workable functions that tie porosity and permeability of the resulting dolomites. Henceforth, three different classes of dolomites are proposed, according to their porosity/permeability relationships (Fig. 3). Each dolomite class will be defined by a function based on the best-fit line of analyzed data (cored plugs). These equations are also reported here (Lucia, 1995): Class 1: K = (45.35 × 10 8 ) Ø ip 8.537 ; r = 0.71 Class 2: K = (2.040 × 10 6 ) Ø ip 6.38 ; r = 0.80 Class 3: K = (2.884 × 10 3 ) Ø ip 4.275 ; r = 0.81 where, K is the permeability in millidarcys and Ø ip is the fractional interparticle porosity.
In the absence of cored plugs to measure the petrophysical properties and multiphase flow properties, a numerical tool that simulates these data by predicting the pore structure evolution of the original grain deposit is of great interest for the petroleum industry and reservoir characterization studies. It allows providing the properties variation in space and times, having information of the intensity of the flow and the surface reaction in times.
The existing models correlating permeability to porosity (such as the Kozeny-Carman model) are of limited use since they lack physical basis for most of the reservoir rock pore network. Indeed, the petrophysical properties of a porous medium are controlled by the pore structure and the 3D connectivity between pores. They may vary within a wide range of values for a similar porosity. This has been observed in many previous studies and provides a major challenge for rock-typing, which is necessary for reservoir characterization, when aiming at proper field development planning and/or Enhanced Oil Recovery (EOR) strategies (cf. Lonoy, 2006).
Recently, PNM has been proposed to study pores-structure evolution and change in petrophysical properties due to particle capture (Ochi and Vernoux, 1999)  Examples of actual reservoir rock textures emphasizing the complexity of diagenetic modification and challenges for numerical modeling: a) Foraminiferal grainstone, overpacked and featuring intra-particle microporosity (mP); b) Dolomitized mudstone/wackestone with intercrystalline microporosity (mP); c) Moldic "and solution-enhanced moldic" porosity (MP) in a wackestone; d) Cementation reduced moldic porosity (CP) in a wacke-packstone. (Sahimi et al., 2000). Baht and Kovscek (1998) used a PNM to describe deposition and dissolution in diatomite. This modeling approach is flexible, and can account for various phenomena occurring at the pore scale. Algive et al. in 2007 and2009 proposed the PNM approach to study the mineral dissolution and precipitation caused by CO 2 sequestration. They developed a PNM tool, dedicated to reactive flow modeling, which has been validated by experimental data. In the present work, this reactive PNM tool is adapted and used to mimic diagenetic cycles (dissolution/precipitation). The methodology is based on three steps. The first step replaces the original complex pore structure of the actual porous media by a conceptual network which consists of a threedimensional cubic lattice formed by pore-bodies (nodes) Schematic representation of porosity evolution with progressive dissolution (solution) and cementation (infilling). This diagram is brought without modification from Choquette and Pray (1970).
interconnected by pore-throats (bonds) respecting the converging-diverging nature of the pores (Fig. 4). The second step consists of solving the governing equations of the precipitation and dissolution phenomena (i.e. reactive convection diffusion equation) in the conceptual 3D pore network. The concentration field taking into account the local surface reactions is calculated for each pore, as for the pressure when solving the flow. Then, the local reactive fluxes and the motion of the fluid-solid interface which control the growth/reduction of the pore-throat and pore-body diameters are deduced. The third step consists of updating the new pore network and calculating the new petrophysical properties of the modified porous media. Those steps are repeated in order to mimic the given diagenetic cycles. Finally, the petrophysical properties and multiphase flow properties of the current porous media are investigated as function of various isothermal diagenetic alterations. The periods of both dissolution and precipitation are fixed. However, various intensities of the alterations are studied in terms of two dimensionless numbers that govern the flow and the reaction regimes (Peclet and Peclet-Damköhler numbers).

Network Model Construction
To simulate the evolution of petrophysical properties occurring during precipitation and dissolution, a numerical model must Schematic representation of Lucia's (1995) petrophysical and rock-fabric classes of dolomites based on similar capillary properties and interparticle-porosity vs permeability values. This diagram is brought without modification from Lucia (1995).

Figure 4
Schematic depiction of a unit cell of the network model. accurately predict the propagation of the mean chemical disequilibrium causing this change and also correctly take into account the local transport and deformation specificities. The advantage of the pore network approach is to reflect the structure key characteristics (connectivity, tortuosity, etc.) while simplifying the geometry features. The PNM approach was originally developed by Fatt in the 1950s to calculate multiphase flow properties of a porous medium. It is now used as a platform to account for various phenomena occurring at the pore scale, including the effects of the surface reaction and mass transfer. The PNM approach was found appropriate to identify the change in structure caused by the dissolution/precipitation mechanisms (Algive et al., 2009). The pore-throats are modeled as simple geometrical tubes of cylindrical, triangular or elliptical form, whereas the pore-bodies are represented by spheres or cubes. In each element (pore-body or pore-throat; Fig. 4), a macroscopic transport equation can be defined and solved analytically. The main advantage of the PNM approach compared to thr reconstructed model techniques (Adler, 1992;Sales et al., 1993;Bekri et al., 1995), which take into account all the details of the pore space geometry and solve the local equations, is that PNM simulations cab run on much bigger objects with less intensive CPU effort.
The numerical simulations are carried out on a synthetic rock that is relatively homogeneous with a narrow pore-throat size distribution ranging from 5 to 10 μm with an average value of 7.5 μm. The pore-body and its pore-throat are chosen to be correlated with a constant aspect ratio value equal to 3.5. The calculated porosity and permeability of this synthetic porous medium are 30% and 200 mD, respectively.
The constructed pore network consists of a three-dimensional cubic lattice formed by pore-bodies (nodes) interconnected by pore-throats (bonds) respecting thus the convergingdiverging nature of pores. The coordination number (bonds per node), which can be variable, is taken here equal to six. The diameter of the nodes and the bonds are randomly assigned according to the pore size distribution. The throats form 45°angles to the main flow direction in the horizontal plane. Periodic boundary conditions are used perpendicularly to the inlet/outlet direction to minimize finite size effects. Additional details on the model characteristics and construction can be found elsewhere (Laroche and Vizika, 2005).

Macroscopic Transport Equation
The macroscopic reactive transport equation that controls the concentration field at the network scale is usually expressed in the following way (Shapiro and Brenner, 1986): where c -' is the normalized mean chemical disequilibrium of the element (i.e. pore-body or pore-throat of a Pore Network Model) and c -0 is the average concentration over the whole network. Three macroscopic parameters appear in the previous transport equation: -γ -* the volumetric apparent reaction rate coefficient, which takes into account source/sink terms due to surface and possibly volume reactions; -v -* the macroscopic velocity vector of the solute, different from the mean velocity of the fluid because of reactions; -D -* the dispersive tensor, which is different that the one describing the classical Taylor-Aris dispersion (Eq. 4). These parameters must be determined for each element from the local data. The upscaling is done at one time for all the geometries used in the PNM by assuming that they are isolated (periodicity). It has been proven (Shapiro and Brenner, 1988) that the macroscopic parameters can be expressed as functions of the three first spatial moments noted m i (with i equal to 0, 1 or 2): (2) (3) (4) There are at least two methods to determine these moments and thus achieve the upscaling: the moment theory and the random walk. The moment technique was first developed to describe the dispersion in a capillary tube (Taylor, 1953;Aris, 1956;Brenner, 1980;Adler, 1992) and was then completed to account for non conservative phenomena (Sankarasubramanian and Gill, 1973;Shapiro and Brenner, 1986). According to the theory, the moments are calculated by solving an eigenvalue problem issued from the local equation system. This system is composed of the local diffusion-convection equation and of a boundary condition reflecting the surface reaction, here assumed to follow first-order kinetics: where D is the molecular diffusion coefficient, v is the fluid velocity vector and κ is the intrinsic rate constant. n is the unit normal to the wall surface S w pointing to the solid phase and c * the equilibrium concentration.
Analytical solutions are limited in practice to twodimensional geometries, such as cylinder or parallel plates. For more complex cases, such as elliptical or triangular capillaries, it is more convenient and robust to numerically determine the spatial moments from a particle cloud propagation using a random walk technique. We have developed both methods and verified that they give the same results for several cases (Algive et al., 2009).

Dimensionless Numbers
The parameters of the macroscopic transport equation have been determined as a function of two dimensionless numbers: -the Peclet number (Pe) which compares the diffusion and convection characteristic times, and -the Peclet-Damköhler number (PeDa) defined as the ratio of the diffusion and the intrinsic reaction characteristic times.
where vis the mean fluid velocity in the network. Bekri et al. (1995) and Algive et al. (2009) demonstrated that, for a low PeDa, the intrinsic reaction kinetics is slow enough for diffusion to uniformize the concentration. Thus, the solid deformation is uniform whatever the flow rate. However, for a high PeDa, the mass transfer limitations become significant. For a high Pe, the solute can be transported up to the small pore-throats where it reacts preferentially since the mass transfer by diffusion is quicker for the elements that have a small radius. This regime favors important permeability changes. On the contrary, for a low Pe, the solute has reacted before reaching the pore-throats and the permeability is less affected.
To determine the mean concentration field at the origin of a mineral dissolution or precipitation inside the network, we use the macroscopic Equation (1). Its parameters take into account the pore-scale information. The concentration of each node of the network is computed numerically by considering the conservation of mass, i.e. the net fluxes entering the node, the accumulation and the sink terms equal zero. The flux expression is derived analytically from the longitudinal concentration profile of the pore-throat by solving Equation (1). The pore-scale information is implicitly considered through effective parameters: conductivities for the flow and the reaction, convection and dispersion factors for the reactive transport (Eq. 2-4). These latter parameters are calculated for each unit cell of the network as described in Section 1.2. The inputs required are the cross-sectional shapes of the porethroats and pore-bodies and the Peclet-Damköhler dimensionless number. They are given as analytical functions in Algive et al. (2009). The conductivities are calculated according to Poiseuille's law (Laroche and Vizika, 2005). The cross-section shapes of the pore-throats and pore-bodies are taken to be equilateral triangle and square, respectively.

Pe
vl D PeDa l D = = and κ

Change in Petrophysical Properties
In the pore network approach, the evolution of the petrophysical properties can easily be determined from the pore-bodies and pore-throats size distribution modifications caused by the formation of precipitation layers or the dissolution. According to the kinetics of the surface reaction (Eq. 7), its stoichiometry and the molar volume V m of the mineral, the solid/fluid interface motion is given by: where F is the wall displacement in its normal direction, ν and ν s are the stoichiometric coefficients of the solute and the mineral, respectively. This deformation is consequently proportional to the chemical disequilibrium at the rock/fluid interface but this one can be expressed as a function of the mean disequilibrium of an element on the basis of the reaction effectiveness γ' (Sharatt and Mann, 1987): (10) γ' measures the decrease of the apparent kinetics due to the mass transfer by comparing its value γ* with the value γ observed in the case of a perfectly mixed system.
Consequently, the model is able to accurately simulate a precipitation/dissolution process in function of the flow (Pe) and the reaction (PeDa) regime. The new pore bodies and pore throats size distribution, determined from the wall displacement F, give a new porosity and a new permeability calculated by applying classical PNM techniques. Thus, permeability and porosity evolutions during the diagenesis can be obtained, and the relative permeabilities and the capillary pressure for various steps of the diagenetic scenario can be computed.

Workflow
The precipitation/dissolution processes are analyzed as follows (Fig. 5). For a given rock structure, the flow field is determined in step 1 by solving the linear system of Poiseuille type equations using the classical pore network approach (Laroche and Vizika, 2005). Then, in step 2, the macroscopic convectiondiffusion equation (Eq. 1) is solved by a Newton-Raphson iterative process, in order to obtain in step 3 the reactive fluxes and the motion of the fluid-solid interface (Eq. 9). Then, in step 4, the rock structure is updated by modifying the pore-throats and pore-bodies diameters according to the precipitation/dissolution fluxes (Algive et al., 2009). These processes can be iterated (step 5) and any quantity of interest, such as porosity, absolute permeability, capilary pressure and relative permeability can be calculated. The absolute permeability is found from Darcy's law when the network is fully saturated with a single phase. Capillary pressure curves are obtained by simulating quasi-static displacement: an increasing pressure is applied on the injected fluid, while the pressure of the fluid in place is kept constant. The two-phase relative permeability curves are calculated at each step of the twophase quasi-static displacements when capillary equilibrium is reached. Flow within each phase is simulated by applying a macroscopic pressure gradient across the network.

RESULTS AND DISCUSSIONS
Various diagenetic alterations (dissolution and precipitation) have been simulated using the pore network modeling. The selected scenario simulates a period of precipitation followed immediately by a period of dissolution. The precipitation and dissolution periods are chosen as follows: a fixed amount of dissolved or precipitated mineral, 25% of the initial porosity, is applied, i.e. the precipitation will be stopped when the porosity is equal to 75% of the initial one, and then dissolution will begin and stopped once the porosity reaches the original value.
However, for real cases, the successive steps of the diagenetic alterations have to be consistent with the geologist's observations, such as petrographical analyses, geochemical measurements and fluid inclusion analyses, which will provide arguments to support the selected diagenetic scenario (time frame of different diagenetic episodes and the reactive fluid involved). Such observations are based on today's increasingly significance of providing tools for quantifying diagenetic processes (e.g. 3D μCT tomography, XRD with Rietveld modeling, and image analysis techniques). Figure 6 shows that when the intensity of the chemical reaction (k in PeDa) is weak the periods of dissolution and precipitation are similar. The initial porosity and permeability are preserved. This can be explained by the fact that both dissolution and precipitation occur uniformly in the porous media when considering a weak chemical reaction. Thus the dissolution mechanism overshadows the precipitation one. To explain the heterogeneity observed in the sedimentary rocks, the intensity of at least one of the processes must be high, i.e. the solute mass transfer onto the wall has to play a significant role.
In the studied scenarios, the precipitation will always be performed at the same reactive flow regime, characterized by PeDa and Pe dimensionless numbers initially set to 1. In this way, we will be able to evaluate the effect of various dissolution regimes on the final petrophysical properties. The pore-body size distribution (Fig. 7, 8) and the pore-throat size distribution (Fig. 9, 10), plotted at different steps of the diagenetic process (original rock, end of precipitation and Final pore structure K -Φ -P c -K r Step 1 Step 2 Step 3 Step 4 Step 5 Figure 5 General scheme of the reactive transport in PNM for the diagenetic alteration. Permeability vs Porosity laws during a diagenetic cycle -Case of a slow chemical reaction (PeDa << 1). Precipitation and dissolution curves coincide perfectly. dissolution periods), explain quantitatively the different precipitation/dissolution configurations and their impact on the petrophysical property evolution (Fig. 11, 12). Various intensities of the dissolution reaction were studied by changing the value of the dimensionless number PeDa (Fig. 7, 9 , 11). The consequences of a change in the flow intensity were examined by changing the value of the Peclet dimensionless number at the beginning of the dissolution process (Fig. 8,   10, 12). As the changes in the final permeability and porosity result from a new configuration of the pore structure, it is expected that the diagenetic process also affects the diphasic properties, such as the capillary pressure and the relative permeability. Figures 13-15 show the impact of different precipitation/ dissolution scenarios on those properties.
Let us first focus on the evolution of the permeability and porosity properties when the rocks are submitted to diagenetic Pore-body size distribution evolution during a diagenetic cycle (semi-logarithmic chart); Pe = 1 -Influence of the reaction intensity (PeDa = 1, 0.1 and 10 for reaction intensity equal to 1, 0.1 and 10 times k, respectively).
alterations (precipitation -dissolution). During all precipitation steps, the reactive flow regime is established to Pe = 1 and PeDa = 1. From a macroscopic point of view, the choice of this regime leads to a deposition along the main flow pathway (Algive et al., 2009), for which the conductivity (throats radii) and the connectivity are higher (Fig. 7). At pore scale, this regime leads to a preferential deposition in pore-throats: their radii get smaller faster than the ones of their adjacent pores (Fig. 9). The dissolution step can take place in a different way, since the reactive flow regime may change.
The variation of the reactive flow regime, once starting the dissolution step, can be simply due to the variation of the pore structure caused by the precipitation step. In fact, the mean velocity (i.e. Pe number) decreases because of the reduction of the conductivity (decrease of the throat radii) during the precipitation step. This specific deposition affects Permeability vs Porosity laws during a diagenetic cycle; Pe = 1 -Influence of the reaction intensity (PeDa = 1, 0.1 and 10 for reaction intensity equal to 1, 0.1 and 10 times k, respectively).  Figure 12 Permeability vs Porosity laws during a diagenetic cycle; PeDa = 1 -Influence of the flow regime (Pe = 1, 0.1, 2 and 10 for pressure drop equal to 1, 0.1, 2 and 10 times dp, respectively). Drainage capillary pressure vs Water saturation; PeDa = 1 -Modification due to the diagenesis -Influence of the flow regime (Pe = 1, 0.1, 2 and 10 for pressure drop equal to 1, 0.1, 2 and 10 times dp, respectively). Drainage oil relative permeability vs Water saturation; PeDa = 1 -Modification due to the diagenesis -Influence of the flow regime (Pe = 1, 0.1, 2 and 10 for pressure drop equal to 1, 0.1, 2 and 10 times dp, respectively). also the initial balance between reaction and diffusion mechanism (decrease of the PeDa). Consequently, the dissolution will essentially take place in the regions of large pore/throat radii, presenting a more important chemical disequilibrium, and will not be extended to the other areas of the pore structure (Algive et al., 2007). Furthermore, because of the unstable character of the dissolution process (Bekri et al., 1995), the dissolved pores or throats increase progressively in size, as can be observed from pore/throat size distributions (Fig. 7, 9). From a petrophysical point of view, such dissolution causes only a minor increase of the permeability. Thus at the end of the diagenetic cycle and for the same initial value of the porosity, the permeability of the altered rock is lower than the one of the original rock (Fig. 11). This can be explained by the fact that the amount of mineral deposited during the precipitation step was not removed during the dissolution: only some pores and throats were affected by both precipitation and dissolution.
This deformation mode, which is related to a macroscopic transport limiting process, can be amplified by increasing the reaction rate constant κ. In this case, the PeDa dimensionless number will be higher (case 6, Tab. 1) which means that the chemical disequilibrium will remain important in few areas of the porous medium only (big pores). Consequently, the deformation that depends on the chemical disequilibrium will not be homogenous (Algive et al., 2009). On the contrary, if the reaction rate constant κ is divided by 10, PeDa becomes lower than 1 (case 5). Since the reaction is the slowest phenomenon for PeDa << 1, the dissolution will take place in a relatively homogeneous manner. As results, the pore and throat size distribution will be shifted regarding the end of the precipitation step of the diagenesis process (Fig. 7, 9). However, as the precipitation was not uniform, the throats of the altered rock (complete cycle of the diagenesis) will be much larger than for the initial rock, except for the smallest ones which do not (or weakly) participate in the flow. Thus, diagenetic processes (precipitation/dissolution) lead, in this case, to an increase of the permeability (Fig. 11). To summarize, due to low value of Pe during the dissolution period, it was found that for a given porosity, precipitation/dissolution can be responsible for a permeability increase if the dissolution rate is slower than the precipitation rate, or a permeability decrease if the dissolution rate is faster. These results are in agreement with the literature (Bekri et al., 1995).
The changes in the reactive flow regime can also be caused by a variation of the pressure gradient (ΔP). The decrease of the interstitial velocity leads to a macroscopic transport limited regime. Thus, for an interstitial velocity or a Pe number 10 times lower (ΔP = 0.1dp), the deformation is located principally in larger pores (Fig. 8, 10). As results, the permeability of the altered rock, at the final step of the diagenesis, will be lower than in the original rock (Fig. 12). On the other hand, the increase of the pressure gradient (increase of Pe number) leads to a relatively homogeneous macroscopic dissolution associated with preferential dissolution in the throats when looking at pore scale. This explains the increase of the permeability of the altered rock (at the end of the diagenesis cycle). Finally, for a moderate increase of the pressure gradient (Pe around 2, ΔP = 2dp), the dissolution occurs essentially on the main flow pathway which favors the development of wormholes (Fig. 8, 10). As results, the permeability of the final rock increases drastically (Fig. 12).
Concerning the capillary pressure and the relative permeabilities, the presence of larger throats in the altered rock favors the drainage process which leads to a lower irreducible water saturation. This effect will be amplified in the case of lower pressure gradient during the dissolution step.
For an equal pressure gradient (case 1, Tab. 1) or 10 times lower (case 2), the throat-radius decrease with regard to the  Drainage water relative permeability vs Water saturation; PeDa = 1 -Modification due to the diagenesis -Influence of the flow regime (Pe = 1, 0.1, 2 and 10 for pressure drop equal to 1, 0.1, 2 and 10 times dp, respectively). original rock (Fig. 10) leads to a higher drainage capillary pressure (Fig. 13). For the same reason, the necessary capillary pressure is lower when the pressure gradient is 10 times higher, ΔP = 10dp (case 3). Finally, the presence of a wormhole when the pressure gradient is moderate (multiplied by 2, case 4) facilitates an easy invasion of the non wetting fluid. However, to invade the other pore-throats, the capillary pressure must be strongly increased, since the radius of the initially smallest throats has decreased during the precipitation process.
The oil relative permeabilities (Fig. 14) of cases 1 and 2 are lower than the original ones because the throats are generally smaller, which restrains the non wetting fluid percolation. On the contrary, the wormhole fingering of the case 4 facilitates extremely the flow of the oil and thus increases its relative permeability. Indeed, the other pores participate only weakly to the flow and their invasion has no influence. For the case 3, the radius increase of the majority of the throats allows a better flow when the oil saturation is still low. However, the end point permeability (K ro max ) is lower because the smallest throats have conductance lower than their initial value. Contrary to the case 4, there is a priori no percolating pathway avoiding these elements of high resistance to flow.
The water relative permeabilities (Fig. 15) are almost identical to the original rock except for the case where wormhole fingering exists. In that case, the fast invasion of the main pathway constrains the water to flow essentially via the smallest pore-throats. This explains the dramatic decrease of the water relative permeability.

CONCLUSION
Pore network modeling approach applied to reactive flow is used to simulate diagenetic alterations (precipitation and dissolution cycle) on a homogenous grainstone characterized by a narrow pore-throat size distribution. A sensitivity study on reactive and flow regimes and their impact on rock properties evolution has been carried out, following a physical approach. The diagenetic precipitation and dissolution cycle was chosen in a way to demonstrate the capability of the model to describe the heterogeneity inherent to the diagenetic alterations commonly observed in carbonate reservoir rocks. The impact of these diagenetic processes (precipitation/dissolution cycles) on the pore structure heterogeneity and consequently on the petrophysical properties (i.e. porosity and permeability) and multiphase flow properties (i.e. relative permeability and capillary pressure) has been investigated. The permeability and porosity evolution for a diagenetic precipitation/dissolution cycle were calculated and analyzed as a function of the relevant dimensionless numbers (Peclet and Damköhler numbers) that characterize the flow and reaction regimes. The correlation between these numbers and the distribution of the dissolved/precipitated layer thicknesses was presented.
It has been observed that, for a given porosity the diagenetic precipitation/dissolution cycle with low flow rate (low Pe) leads to permeability increase if the dissolution rate (κ) is lower than the precipitation one, and to permeability decrease in the opposite case. The intensity of the flow rate also entails changes in the pores network and consequently in the petrophysical properties. The increase of the flow rate (higher Pe) allows the dissolution to occur essentially on the main flow pathway. As a result, the permeability of the altered rock (end step of the dissolution period) increases drastically.
Evolutions of two-phase flow properties (due to the new configuration of the pore structure) were calculated, applying classical quasi-static two-phase flow pore network modeling on altered and non altered porous media. The capillary pressure (P c ) and the oil relative permeabilities (K ro ) are related to the altered pore-throat size distribution. The presence of the large pores or wormhole allows an easy invasion of the non wetting fluid (oil) which leads to lower entrance P c and higher K ro at high water saturation. The water relative permeability curves are identical to the ones of the non altered rock, except when wormholing appears.