A complete workflow applied on an oil reservoir analogue to evaluate the ability of 4D seismics to anticipate the success of a chemical enhanced oil recovery process

. In an Enhanced Oil Recovery (EOR) process, one of the main dif ﬁ culties is to quickly evaluate if the injected chemical products actually improve oil recovery in the reservoir. The ef ﬁ ciency of the process can be monitored in the vicinity of wells, but it may take time to estimate it globally in the reservoir. The objective of this paper is to investigate the ability of 4D seismics to bridge this gap and to help predict the success or breakdown of a production strategy at reservoir scale. To that purpose, we consider a complete work ﬂ ow for simulating realistic reservoir exploitation using chemical EOR and 4D seismic modeling. This work ﬂ ow spans from geological description to seismic monitoring simulation and seismic attributes analysis, through geological and reservoir modeling. It is applied here on a realistic case study derived from an outcrop analog of turbiditic reservoirs, for which the ef ﬁ ciency of chemical EOR by polymer and surfactant injection is demonstrated. For this speci ﬁ c ﬁ eld monitoring application, the impact of both water ﬂ ooding and proposed EOR injection is visible on the computed seismics. However, EOR injection induces a more continuous water front that can be clearly visible on seismics. In this case, the EOR ef ﬁ ciency can thus be related to the continuity of the water front as seen on seismics. Nevertheless, in other cases, chemical EOR injections may have more moderate impacts, or the ﬁ eld properties may be less adapted to seismic monitoring. This points out the importance of the proposed work ﬂ ow to check the relevance of seismic monitoring and to design the most adapted monitoring strategy. Numerous perspectives are proposed at the end of the paper. In particular, experts of the different disciplines involved in the proposed work ﬂ ow can bene ﬁ t from the availability of a complete set of well-controlled data of various types to test and improve their own tools. In contrast, the non-experts can easily and quickly bene ﬁ t from “ hands-on ” experiments for understanding the involved phenomena. Further-more, the proposed work ﬂ ow can be directly applied to geological reservoirs all over the world. injection to estimate the ef ﬁ ciency of the chosen products at reservoir scale and to optimize the production scheme. Experiments are ﬁ rst conducted


Introduction
Crude oil, representing more than 30% of energy demand, is still a major energy source for the current economics of the world (Manrique et al., 2010). Hence, to satisfy the increasing global energy demand and consumption, increasing effort must be devoted by the oil industry not only to identify new oil reservoirs, but also to improve the oil recovery from mature reservoirs. Primary oil recovery, by initial reservoir pressure depletion, hardly exceeds 20% of Original Oil in Place (OOIP). Water injection is applied in a secondary oil recovery to prevent depletion of the reservoir pressure and to recover more oil. However, more than roughly 65% of OOIP remains unswept after this stage.
In the next step, often referred to as tertiary oil recovery, Enhanced Oil Recovery (EOR) or Improved Oil Recovery (IOR), a substance that is not present in the reservoir is injected in order to further increase oil recovery (e.g., Cossé, 1993;Green and Willhite, 1998;Lake et al., 2014).
Among the existing techniques, chemical based EOR including synergetic combination of surfactant, polymer and alkali has been well documented either in laboratory studies (e.g., Bourbiaux et al., 2015;Hou et al., 2005;Leray et al., 2016), in simulation studies (e.g., Douarche et al., 2012Douarche et al., , 2014Zerpa et al., 2005) or in pilot tests (Delamaide et al., 2014;Demin et al., 1997;Vargo et al., 2000). Due to the cost of injected chemical products, various tests are performed prior to field injection to estimate the efficiency of the chosen products at reservoir scale and to optimize the production scheme. Experiments are first conducted at core scale (Borozdina et al., 2019). They provide an idea of the oil recovery that can be expected at the field scale. As the reservoir heterogeneity can strongly affect the efficiency of the EOR strategy, it needs to be also estimated and optimized at larger scalespilot and field scalesprior to production. This can be achieved using numerical models representing the reservoir.
To be as representative as possible, these models should be constrained to all available data. In particular, they should reproduce the production measurements acquired at wells such as pressure and injection/production rates. However, these data may not be informative enough to properly constrain all the model parameters, such as petrophysical properties. As a result, the production forecasted by the models in response to a given development plan may differ from the one that would be observed if actually producing the reservoir following this strategy. It is thus necessary to check the in-situ efficiency of the EOR process as soon as possible after the beginning of chemical product injection to avoid useless injection of costly products. Single-well chemical tracer tests can provide information in the vicinity (few meters) of the injection wells (Cockin et al., 2000;Deans, 1971;Deans and Carlisle, 2007). Oil production measurements can also be used, but the impact on these data may take time to become clearly visible depending on the reservoir properties and well location. It would thus be interesting to identify a monitoring technique able to provide information on the efficiency of the EOR strategy away from the wells and in a short time frame after the beginning of injection. The associated measurements could also be compared to the corresponding simulated values and used if necessary to update the reservoir model. Time-lapse seismics or 4D seismics, naturally comes to mind as a potential method. 4D seismics involves acquisition, processing and interpretation of repeated seismic surveys over a field during its exploitation. The objective is to determine the changes occurring in the reservoir by comparing the repeated datasets (e.g., Calvert, 2005;Johnston, 2013). This method moved from a research experiment tool in the mid-1980s to a mature technology now widely used for instance to monitor subsurface exploitation for assisting reservoir management and field development. A detailed description of this technology is beyond the scope of the present paper and can be found in the two previous references. However here we shall introduce the three main high values of 4D seismics unambiguously demonstrated by field achievements.
The first one is the ability to locate by-passed oil and undrained reserves. This has been observed in most of the major oil fields where 4D seismics have been applied, either in the Gulf of Mexico on Hoover, Madison or Marshall fields (e.g., Helgerud et al., 2011), in the North Sea (e.g., Koster et al., 2000;Landro et al., 1999), or in offshore West Coast of Africa (e.g., Lumley et al., 1999;Onuwaje et al., 2009). Identifying undrained compartments is the second value of 4D seismics. For instance, on the Gullfaks field in Norwegian North Sea, 4D seismics helped to identify not only production related changes of seismic amplitudes but also un-swept zones in the reservoir. Later, two wells were drilled in the zones identified as undrained by 4D seismics. Not only both wells encountered oil filled reservoirs but also well data confirmed the waterfront movement predicted in the reservoir (Landro et al., 1999). The last main value of 4D seismics that we shall consider is its ability to optimize the infill drilling program. For instance, on Draugen field offshore Norway, the flow model was updated manually to predict the performance of the planned infill at several alternative locations. 4D seismics contributed to find the optimal location that differs from the originally planned location. A straightforward consequence was the resulting delay in anticipated water cut which extended the forecasted plateau production and increased the predicted recovery (Koster et al., 2000). Most of these values are confirmed in the case considered in the present work.
In this paper, we thus propose to address the question of the ability of 4D seismics to characterize the impact of a chemical product injection in an oil reservoir. More specifically, we propose to compute 4D seismic cubes for different production strategies on a given reservoir to see if characteristic features can be identified and related to each strategy. This appears as an important preliminary step before real seismic data acquisition to avoid useless and costly measurements and to design the most informative 4D survey. The computation of the seismic response is performed here with the workflow described in Figure 1. It consists of geological and reservoir modeling, flow simulation, rock physics modeling of the petroelastic properties, and seismic modeling, processing and seismic attribute analysis. This workflow was already partially considered for various applications. For instance, geological modeling and reservoir simulation are considered on a polymer flood pilot in Delaplace et al. (2013). These two steps are complemented with impedance variation computation in Roggero et al. (2012), Le Tillier et al. (2012), considering reservoirs produced by water or gas injection, and SAGD process.
Our study focuses on the application of the workflow on Ainsa-1 quarry outcrop (South-central Pyrenean foreland basin in Spain), a well-known analog of turbiditic reservoirs offshore West Africa (e.g., Arbués et al., 2007;Falivene et al., 2006b;Pickering et al., 2015). Outcrop analogs have been studied for decades to constrain the physical geological properties of subsurface reservoirs and better understand depositional environments at the origin of petroleum systems and their architectures (Pringle et al., 2006). Outcrop analog studies help hydrocarbon exploration as they provide detailed characterizations of geology and facies distribution that can be used to improve seismic dataset interpretation and fluid dynamic understanding (Massonnat et al., 2017). Models built from outcrop studies can be used to generate synthetic seismic data which can serve in turn as training tools for seismic interpretation or can be compared to real seismic-reflection data acquired on the field of interest (Armitage and Stright, 2010;Bakke et al., 2013;Bourgeois et al., 2004;Doherty et al., 2002;Falivene et al., 2010;Holgate et al., 2014;Janson and Fomel, 2011;Schmitz et al., 2014;Schwab et al., 2007;Stright et al., 2014). The quality of 3D facies models built from outcrop analogs is also evaluated through the ability of such models to reproduce or predict the dynamic flow behavior of subsurface related systems (Adams et al., 2011;Falivene et al., 2006a;Schmitz et al., 2014). Reservoir dynamic simulations performed on static models from outcrop analog studies are an excellent way to understand the impact of good-reservoir and non-reservoir rock-type distributions on sweep efficiency for instance (Adams et al., 2011). It is also a very interesting way to test the effects of different scales of modeled heterogeneity on production (Falivene et al., 2006b) and/or to optimize well placement in the related types of reservoirs (Stright et al., 2014). So our work benefited from previous studies of this outcrop (Arbués et al., 2007;Falivene et al., 2006b;Schmitz et al., 2014) and strongly emphasized the integration of the various fields of expertise involved in this workflow; for instance in Virtual analog (e.g., Schmitz et al., 2014), in Geological modeling (e.g., Gasparrini et al., 2016), in Reservoir simulation & history matching (e.g., Le , in 4D seismic inversion (e.g., Labat et al., 2012) and in Rock Physics (e.g., Zinszner, 2012, 2014).
The paper outline is as follows. First, we detail the geological and reservoir modeling steps of the proposed integrated workflow. Then, we describe the flow simulation part and the results obtained with two injection strategies: a chemical EOR process, namely Surfactant-Polymer (SP) injection, and the more conventional Water-Flooding (WF). Contributions of 4D seismics to anticipate the success or breakdown of the chemical EOR strategy are then developed. The last sections are dedicated to discussions on the achievements of the proposed workflow, and to some concluding remarks and possible extensions of this work.

Geological and reservoir modeling
The model is constructed from the numerous geological observations on an outcrop located close to Ainsa town, in southern Pyrenees, Spain. This outcrop has been extensively studied as an analog of reservoirs offshore West Africa, and more generally for the overall understanding of turbidite systems (e.g., see Arbués et al., 2007;Pickering et al., 2015). Using internal software (OpenFlow Suite, 2015 and an adapted workflow, we benefited from the extensive work of Schmitz et al. (2014) for achieving a high resolution 3D virtual reconstruction of the outcrop using photogrammetric methods. The dimension of the outcrop is about 40 m thick and 750 m wide. It is oriented NNW-SSE, oblique to the mean paleoflow estimated towards WNW.
The Ainsa-1 turbidite system is subdivided into three cycles of channel-complex development and abandonment. Channel complexes are composed of closely stacked channel forms of several meters thick to 20 m thick (Fig. A1 of Appendix A). Based on geological interpretation and outcrop studies, the outcrop section is divided into five sedimentary zones by six bounding horizons. This division is based on erosional surfaces of turbidite channels, because their filling is multistory and includes facies associations that represent various degrees of erosion and sediment bypass by turbidity currents (e.g., Arbués et al., 2007).
Five facies were identified by Arbués et al. (2007), namely from the more to the less permeable, thick-bedded sandstone, heterolithics, conglomerate, mudstone clast conglomerate, and gravelly mudstone. Conglomerate and mudstone clast conglomerate are up to 1-m thick while thick-bedded sandstone has sandstone beds thicker than 10 cm and heterolithics have finer sandstone beds (up to 10 cm). More details are given in Appendix A.
The facies model is built on the same geological grid as in Jardin et al. (2010). It is defined according to the five sedimentary zones, and consists of about 650 000 active cells of size 15 m Â 15 m in the horizontal direction. Vertical cell size is not constant and equals 0.5 m on average. The vertical resolution is thus adapted to correctly model facies distribution except for heterolithics for which hypotheses (in terms of permeability, porosity. . .) have to be made to model the impact of mudstone beds contained in this facies. Facies proportions are deduced from four vertical sedimentological sections selected across the outcrop (Appendix A) also used to build the facies model obtained by Jardin et al. (2010). To obtain more spatially distributed data, four additional synthetic vertical logs (Fig. A3) were extracted from this facies model and used together with the previous data to generate a new model in the present work. The geostatistical approaches used for this simulation are the Truncated Gaussian method (Doligez et al., 1999;Galli et al., 1994) for units C1 and C2.2, and the sequential indicator simulation for units C2.1 and C3 (Fig. A1b). These geostatistical approaches have already been applied successfully in many depositional settings, including deep offshore settings (Jardin et al., 2010;Lerat et al., 2007). Referring to Falivene et al. (2006a), the values for major, minor and vertical direction ranges of variograms were considered to be the same for all facies, 500 m, 100 m, and 0.5 m respectively, with an azimuth parallel to the paleoflow direction for both approaches.
In Schmitz et al. (2014) and Falivene et al. (2006b), constant petrophysical properties were attributed to each facies. In our study, porosity and permeability are considered as random functions defined per facies. Truncated Gaussian distributions are used for porosity and truncated lognormal distributions for permeability. The geostatistical characteristics (e.g., mean and variance as detailed in Appendix A) were chosen to be representative of analogous subsurface turbidite reservoirs while being relevant to induce impacts on chemical EOR. We wanted to induce contrasted permeabilities between the thick-bedded sandstone and the heterolithics to highlight flow layering reduction thanks to chemical product injection. Figure 2c shows the resulting 3D porosity model. Figure 3a illustrates the statistical distribution of porosity and horizontal  The most porous/permeable facies (typical porosity / > 0.10 and horizontal permeability k H > 10 mD) are the thick-bedded sandstone, the conglomerates, and, to a lesser extent, the mudstone-clast conglomerates. The other facies, namely heterolithics and gravelly mudstone, are substantially less porous/permeable, within the reservoir and in the main E-W channel. Each facies permeability is assumed isotropic except for heterolithics and mudstone-clast conglomerates for which the vertical permeability is correlated to the horizontal one using a factor of 0.1 and 0.05 respectively because of mudstone beds or inclusions. For comparison, Figure 4 shows horizontal slices of the facies, porosity, and horizontal permeability models, located approximately halfway up the reservoir. The presence of a channel filled by debris-flow deposits in the southern part is represented by a narrow zone of quite low porosity/permeability. Usually, at this stage, an upscaling of petrophysical properties from the fine geological scale to the coarse reservoir scale is performed in order to speed up flow simulations (e.g., Iske and Randen, 2006). Here, the size of the model makes it quite manageable for computation, so that no upscaling was performed. Thus, in our case, the reservoir grid was identical to the geological grid. We are aware that the flow simulations are usually performed at a coarser scale than geological modeling, mainly because of computation cost. Here we perform flow simulations at the same scale in order to point out subtle effects such as layering effect due to the heterogeneity of facies and permeability distributions in the reservoir as will be detailed further. This will allow a refined interpretation of the seismic response, and as a consequence of the seismic visibility analysis. Regarding the impact of the flow simulation on the seismic visibility analysis, we expect that the randomness and the characteristic statistical lengths of the studied media have a direct influence on the answer. In other words, more homogeneous media do not necessitate flow simulation at a fine scale. Beyond such elementary consideration, we agree that a detailed analysis of the impact of the flow simulation grid resolution on seismic feasibility analysis would be a research topic in itself but is beyond the scope of the present work. Finally, average compressibility values are assigned to each facies in order to compute the pressure field during flow simulation. Rock bulk modulus is of the order of a few GPa to tens of GPa (Rasolofosaon and Zinszner, 2012), which corresponds to compressibility of the order of 10 À4 -10 À5 bar À1 (details can be found in Appendix A). We chose an average reservoir depth of 2200 m corresponding to an average overburden weight of 600 bars, roughly three times the hydrostatic pressure. The initial fluid pressure was then 253 bars and the reservoir temperature was 83°C.
As far as fluid composition is concerned, the reservoir is naturally saturated with an oil characterized by an average density of 853.7 kg/m 3 (approximately corresponding to API gravity of 34) at surface conditions (the water-oil contact is located below the reservoir). For flow simulation (step 3 in Fig. 1), we assume a two-phase saturating fluid with oil and water (average density of 1012 kg/m 3 at surface conditions), but no gas. Furthermore, the dependence of fluid properties required for flow simulation is given in Appendix A. Relative permeability curves and capillary pressure curves are typical of water wet rocks. For all facies, water saturation varies between 30% and 68% except for sandstone in which water saturation varies from 15% to 75.5%. Although relative permeability and capillary pressure end-point values are the same for all facies, the different irreducible oil and water saturations make the sandstone facies the best one in terms of oil volume storage and facility to produce it. More details, and in particular endpoint curve values, are given in Appendix A.

Flow simulation and injection strategies
As illustrated in Figure 5, the production pattern was an inverted five-spot (e.g., Lyons and Plisga, 2011), with a single injection well surrounded by four producers. The production wells were located roughly in the corners of the model, two of them being placed south of the main E-W channel and aligned with the channel axis. The injector well, at the center of the pattern, was located to the north of the channel. Such a configuration allows us to check the possible compartmentalization of the reservoir due to the presence of the main E-W channel. Productivity indices for these wells were classically estimated by the method of Peaceman (1983). Injection and production were imposed target maximum rate constraints (Q max = 600 m 3 d À1 for the injector and 150 m 3 d À1 for each producer, see details in Appendix B) with a switch to pressure control if the pressure limit was reached (maximum pressure at the injector and minimum pressure at the producer). Pressure constraint at injectors aims at protecting both wells and reservoirs from mechanical damages while for producers, this pressure limit is a flow requirement to the production facilities. Such a switch to pressure control was observed here at the injector at the beginning of polymer injection (see Appendix B). More precisely, the injector pressure increased to the maximum authorized pressure, and was then kept below or equal to this limit pressure by decreasing when necessary the water injection rate.
The spatial distribution of oil saturation in the reservoir after 14 years of water flooding is shown in Figure 5. Even after such a long time, some areas have been hardly swept, with an oil saturation as high as 0.8. This is the case for instance to the south of the E-W channel (near producers #3 and #4) and in the vicinity of well #1. The presence of the major E-W channel contributes to compartmentalize the southern part of the reservoir. That will be confirmed by time-lapse seismics in the next section. Furthermore, long-range heterogeneity features, such as  vertical interbedded contrasts, have significant effects on flow behavior. This is visible on the top of Figure 6 showing a vertical NS cross-section of oil saturation simulated after only 42 months of water-flooding. The red ellipse points out a layering effect due to the facies heterogeneity and the associated petrophysical property distributions in the reservoir.
This layering effect and the large volume of oil remaining in the reservoir at the end of the water-flooding process gave support to our interest in studying the ability of chemical EOR to improve areal sweep efficiency. Different chemical EOR scenarios with various types of combinations of alkaline (A), surfactant (S), and polymer (P), namely P, S, SP, AS, and ASP were tested and compared to a conventional water-flooding. A detailed description of the results is beyond the scope of this work but can be found in Zarate-Rada (2016). Here, we only focused on the two injection strategies described in Figure 7: Polymer-Surfactant (SP) injection and pure Water-Flooding (WF).
In both cases, water injection started in April 2016. The EOR process started a year later, in April 2017, just after an early water breakthrough. Figure 8 compares the two scenarios in terms of cumulative oil production and water cut during the 14 years of production.
The curves corresponding to WF and SP are plotted in blue and green, respectively. As expected, the SP scenario induced a better recovery than simple water-flooding, with a strong reduction of the water cut and a larger cumulative oil production. As illustrated in Figure 6, this improvement is due to a decrease of the flow layering induced by the    higher viscosity of the injected fluid leading to a reduction of mobility.
Performing this kind of comparative numerical studies is a classical step to estimate the added value of a given chemical formulation and injection strategy compared for instance to water flooding. However, these tests are based on a numerical model that, although calibrated on the available data (oil production, water cut, bottom hole pressure. . .), may provide production forecasts that differ from the real behavior. The in-situ efficiency of the EOR process should thus be checked as soon as possible after the beginning of injection. In practice, the impact on the reservoir production observed from the wells may take time to become clearly visible depending on the reservoir properties and well locations. We thus propose to evaluate the ability of 4D seismics to assess the efficiency of the EOR strategy at an early stage of the process, close to the beginning of injection.

Synthetic time-lapse seismic response
Two specific calendar times, namely Baseline (04/04/2017) and Monitor (04/14/2018), are chosen to compute the timelapse seismic response. Baseline calendar time corresponds to the beginning of polymer injection. Monitor calendar time was chosen one year after (Fig. 7). The time-lapse seismic response is computed based on petroelastic properties derived from the pressure and fluid saturation variations. In the following sections, the impact of the SP scenario on these intermediary properties is analyzed prior to conclude on the seismic response itself. The idea is to make the link between the characteristic features identified from the reservoir properties and the seismic attributes.

Pressure and oil saturation
To obtain a numerical estimation of time-lapse seismic response, pressure, and oil saturation have first to be obtained at the two dates from the flow simulations described in the previous section. Figure 9 displays the changes in the reservoir pressure and oil saturation in response to the two injection strategies. Note that the representation in difference cubes (Monitor -Baseline) was usually preferred to the representation in monitor cubes in order to magnify the visualization of the production process signature (e.g., Calvert, 2005;Johnston, 2013). Each sub-figure shows a horizontal slice of the considered 3D difference cube, at Monitor calendar time. The upper series corresponds to WF and the lower one to SP process, while the first column corresponds to the variation of oil content DS o and the second one to the variation of fluid pressure DP.
Let us notice the complicated pattern of the spatial distribution of oil content DS o for both injection strategies, mainly due to the reservoir heterogeneity. The signature of a horizontal spatial correlation of facies, porosity and permeability roughly in the SE-NW direction seen in Figure 4 can be retrieved in the DS o distributions in Figure 9, with alternating zones of increasing (light yellow) and decreasing saturation (light green). In addition, switching from water flooding to polymer injection induced a substantial change in the spatial distribution of both DP and DS o . For oil saturation, the light green area around the injector well in the case of water-flooding corresponds to the area swept during the first year of water injection (before Baseline calendar time). This region was swept a little bit more by water during the second year of production (between Baseline and Monitor calendar times), while the area newly swept during this period appears in blue. Of course, for SP, the size of the area swept before Baseline calendar time is the same as for WF, but the oil saturation variation inside this zone becomes significantly different after polymer injection. In particular, the yellow ring near the injector (and located between the two doted circles) corresponds to a local oil saturation increase due to the sweeping of additional oil which had been mobilized in the direct vicinity of the injector since the beginning of polymer injection. The spatial distribution of the fluid pressure variation DP is substantially smoother than the one of DS o . This is a classical result due to the fact that pressure variations follow diffusion processes which induce quite smooth spatial variations (e.g., Bear, 1972;Cossé, 1993;Houpeurt, 1975;Muskat and Wyckoff, 1937). Another clear result is the positive pressure variation DP for WF, which contrasts with the negative DP for SP. This result, unintuitive at first sight, can be explained by the different water injection rates for the two scenarios (WF and SP). Indeed, as polymer is less movable than water, polymer injection induced an increase of pressure at the injector compared to water flooding. The maximum authorized pressure was reached in this case so that the water injection rate decreased for the SP scenario, while remaining constant in the WF scenario (see Appendix B). Concomitantly, the area already swept by water before Baseline calendar time became an underpressure zone due to this decrease of water injection rate without changes of the production rates. As a result, pressure in this zone decreased during polymer injection (Fig. 10), leading to a negative DP between Baseline and Monitor calendar times.
This phenomenon also explains the fact that the SP waterfront at Monitor calendar time is closer to the injector well than the WF waterfront. Nevertheless, despite the change in injection strategy between WF and SP, the yellow ring is really associated with polymer injection. In this modeling, certain areas have already been swept on the Baseline date with a strong layering effect. As the injected water flows preferentially inside the swept zones (physical effect of water/oil diphasic flow), a hypothetic reduction in the injection rate for the WF case, as observed for SP, would not succeed in mobilizing additional oil. For SP case, the polymer also flows inside these swept zones but partially obstructs them, and then helps to better sweep neighboring areas. If the injection rate was smaller, the layering effect would be smaller at the Baseline date. In this case water would sweep the reservoir better while the polymer would be less efficient (this point is discussed in Sect. 5).
We now propose to check whether these differences in pressure and saturation variations also induce different characteristics on the petroelastic attributes.

Petroelastic modeling
The petroelastic modeling (step 4 in Fig. 1) is the key element to obtain a numerical modeling of the time-lapse seismic response. During this step, fluid properties (pressure, saturations, fluid densities, fluid elastic properties) and rock properties (porosity, matrix elastic properties) are converted into simulated elastic responses (P-and S-wave velocities, impedances) which are used to build the seismic responses (seismic amplitudes and other relevant seismic attributes). Previous studies proposed seismic models from outcrops-derived 3D facies models of the Eocene Ainsa turbidite systems (Falivene et al., 2010;Schmitz et al., 2014), but since they did not couple flow and seismic modeling for monitoring purposes, constant elastic velocities and densities were assigned to each facies. In contrast, in our study, the variations of seismic properties of the reservoir at any point due to fluid substitution and pressure variations induced by production were computed using Biot-Gassmann type of petroelastic model (e.g., Rasolofosaon and Zinszner, 2012). As detailed in Appendix C, the model takes into account facies average petroelastic parameters as well as the pressure dependence of fluid and rock densities and elastic properties (see Figs. C1-C3 in Appendix C). In summary, fluid substitution (oil by water in this case) induces substantially larger variations of the fluid-saturated rock bulk modulus M (sat) than fluid pressure variations. For instance, pressure-induced variation of M (sat) hardly exceeds 5% for heterolithics, and even 1.5% for thick-bedded sandstone, in the most extreme case, namely in the case of full-oil saturation (S o = 1 À S wi ) (see Fig. C3 in Appendix C). In contrast, the variations of M (sat) induced by saturation changes due to the substitution of oil by water can reach up to 23% for thick-bedded sandstone and 20% for heterolithics.
Furthermore, the reservoir rock response to simultaneous fluid substitution and fluid pressure variation is intricate due to two competing effects. On the one hand, an increase of the fluid pressure P p , implying a decrease of the differential pressure P diff at fixed confining pressure P c , tends to soften the rock (decrease of the bulk modulus) by opening the cracks and compliant pores. On the other hand, an increase of the fluid pressure P p tends to increase the fluid bulk modulus of the saturating fluid, as illustrated by Figure C1, which stiffens the fluid-saturated rock. Depending on the level of oil saturation and on the considered rock, an effect dominates more or less the other one. In such cases, phenomena are so intricate that intuition is not sufficient. This highlights the need for a complete numerical computation in the framework of the proposed integrated workflow. This is all the more true that additional complexity induced by the presence of substantial heterogeneity is ubiquitous in the model. From Biot-Gassmann theory for fluid substitution and Hertz-type theory for pressure dependence, 3D petroelastic cubes were computed. The different inputs for the petroelastic model of the different facies are summarized in Appendix C. A cube of Acoustic Impedance (IP), a cube of Shear-wave Impedance (IS), and a cube of P-wave Velocity (VP) for Baseline and Monitor calendar times are obtained (Biot, 1941;Bourbié et al., 1987;Gassmann, 1951;Rasolofosaon and Zinszner, 2014). Figure 11 shows the impact of fluid substitution and pressure variation on I P variations (DI P ) for both injection strategies (WF and SP). The major seismic impedance variations were due to fluid substitution, whereas DI P was also sensitive to fluid pressure variation beyond the swept zone but to a lesser extent. As for DS o , the signature of a horizontal spatial correlation of facies, porosity and permeability roughly in the SE-NW direction seen in Figure 4 can be observed on the DI P . Furthermore, we can note the absence of impedance variations to the south of the E-W major channel (Fig. 4), denoting a weakly swept zone corresponding to the vanishing DP and DS o due to the compartmentalization of the reservoir.

Synthetic 4D seismic cubes and relevant attributes
From each triplet of petroelastic cubes I P , I S , and V P , and for each calendar time, three P-wave seismic cubes were computed by convolution (e.g., Sheriff, 2002) using a high frequency Ricker source (typical 150 Hz central frequency due to the small thickness of the reservoir). Figure 12a illustrates such a process for the case of normal incidence angle (h = 0°), with a 3D view of the difference between the Monitor and the Baseline seismic cubes. The conversion of the vertical scale from time domain to depth domain is a problem in itself which falls beyond the scope of this work but is a well-known topic (e.g., Bartel et al., 2006;Iversen and Tygel, 2008;Yilmaz, 2001). A very simple way to convert depth to time is to use a P-wave velocity. In our case, this P-wave velocity was computed in the petroelastic modeling step of our workflow. Nevertheless, a critical issue remains and has to be carefully handled in real 4D seismic studies: the time correspondence between all seismic cubes acquired at different calendar times. Between two consecutive seismic acquisitions, the fluid content and the pressure field are different, leading to different P-wave velocities for these two different acquisitions. As a consequence, as illustrated by Figure 13, a horizon at a given depth appears at different travel times for these two different acquisitions, impeding any direct comparison.
This direct comparison can only be done after seismic cubes have been put in time correspondence, as shown in Delépine et al. (2010). In our realistic synthetic case, the time-shift law from Baseline time to Monitor time could be computed easily from the V P cubes, which allows to put in time correspondence the Monitor cube with the corresponding Baseline cube. All the difference cubes shown on Figure 12 are in the Baseline time.  We can appreciate on Figures 12 and 14 the loss of vertical and horizontal resolution (compared to impedance or saturation variations for instance) due to the limited frequency content of the seismic source. Nevertheless, it is still possible to observe on 4D monitoring data, highlighting changes in the reservoir, the signature of horizontal spatial correlation of facies, porosity, and permeability which correspond to static characteristics of the reservoir.
At last, Figures 12b and 12c show two of the most relevant attributes of the 4D seismic data, as far as our case study is concerned. These attributes are based on the notion of "Normalized Root Mean Square" (NRMS). In seismics, the NRMS between a first seismic signal (Monitor signal) and a second signal (Baseline signal) is defined by: where RMS(Á) stands for the "Root Mean Square" of the considered signal. In 4D seismics, the value of NRMS is a well-known quantification of the repeatability of a seismic experiment (e.g., Johnston, 2013). By definition, NRMS can take a numerical value comprised between 0 and 2. In the case of two quite similar seismic signals, the corresponding NRMS value is close to 0. In contrast, in the case of two quite similar seismic signals but of opposite polarity, the corresponding NRMS value is close to 2. Figure 12b shows the NRMS of the seismic difference (Monitor -Baseline) corresponding to WF scenario, and Figure 12c shows the NRMS of a specific seismic attribute called Hilbert transform (e.g., Sheriff, 2002;Yilmaz, 2001). Hilbert transform, and more exactly the modulus of Hilbert transform of a seismic signal, is a "clean" way, from the physico-mathematical point of view, to quantify the instantaneous amplitude of the seismic wave; the square of this quantity being proportional to the instantaneous energy transported by the seismic wave. In effect, on both NRMS cubes plotted in Figures 12b and 12c, the NRMS value is quite small in zones not affected by fluid substitution, such as areas in the southern part of the E-W channel or in the overburden. In contrast, NRMS has significant values in areas where seismic properties were affected by injection, either due to fluid substitution or to pressure effect. Note that all these RMS were computed on a time window of fixed width vertically sliding with overlap. The choice of the width of the time window for computing the RMS is important. A too narrow time window can strongly affect the statistical representativeness of the computation. In contrast, a too broad time window contributes to critically decrease the vertical resolution  Horizontal cross-sections halfway up the reservoir for WF (top line) and EOR injection (SP) (bottom line) process on 4D seismic response (1st column), and NRMS of seismic cube difference (2nd column). Note that the loose term "Dseimic" stands for the difference between the seismic data in the amplitude domain. of the analysis. Here, the width of the time window was close to the period of the seismic source signal.
A horizontal slice of NRMS of the seismic difference cubes (Fig. 12b) corresponding to WF is provided on Figure 14 and can be compared to the same attribute computed for SP scenario. In the case of WF scenario, we can see, beyond the waterfront, the impact of the pressure increase on NRMS. For SP scenario, the NRMS of the seismic difference nicely delineates the swept zone. Indeed, the polymers enhance the viscosity of the injected fluid which induces a slower displacement of the waterfront. The layering effect of water in permeable facies is reduced, so that the waterfront is less scattered horizontally from one layer to the other and appears easier to detect on 4D seismics (less noisy). In this case, the visibility of the waterfront on the real seismic cubes may thus be used to estimate the efficiency of the EOR process.

Impact of the injection rate on the seismic attributes
For a given EOR strategy, the injection scenario can also have a non-negligible impact on the seismic attributes. To illustrate this, we considered a second SP injection scheme that differs from the previous one by the imposed target maximum fluid injection rate. It was taken equal to Q max = 350 m 3 d À1 instead of 600 m 3 d À1 . The corresponding variations of oil saturation DS o , pressure DP and seismic response between the Monitor and Baseline calendar times are given in Figure 15. The pressure constraints discussed above and detailed in Appendix B, were unchanged. We can see that both DP and DS o were substantially influenced by this change in Q max . For the low rate case (Q max = 350 m 3 d À1 ), high DP values close to the injection well location were observed, while such a behavior did not occur for the high rate case (Q max = 600 m 3 d À1 ). This is again due to the pressure constraint at the injector. Indeed, at Baseline calendar time, the pressure around the injection well was lower for Q max = 350 m 3 d À1 than for Q max = 600 m 3 d À1 , but when polymer was injected, the well bottom-hole pressure reached the maximum authorized pressure of 500 bar whatever the injection rate scenario (see Fig. B2 in Appendix B). That is the reason why DP values were higher for the low rate case than for the high rate case. As far as variation of oil saturation is concerned, a low injection rate led to a better oil sweep (light blue to blue zones) in the area swept between Baseline and Monitor calendar times. The area swept with water before Baseline calendar time is also better swept using Q max = 350 m 3 d À1 and the layering effect observed at the Baseline date is lower than the one observed for Q max = 600 m 3 d À1 . This behavior is not general and depends on the reservoir properties and injectionproduction strategy and control. If we consider, for this period, the effective water injection rate for the two injection scenarios (Fig. B2 in Appendix B), the scenario Q max = 350 m 3 d À1 was the one for which the larger amount of fluid was injected, inducing a better oil recovery. As a consequence, the polymer has a lower and more diffusive impact, which induces a lower quantity of remobilized oil and a finer and less continuous yellow ring than for Q max = 600 m 3 d À1 .
On Figure 15, we can see that the seismic signature was sensitive to this peculiar behavior and, even if the scenario Q max = 350 m 3 d À1 led to a better sweep efficiency between Baseline and Monitor calendar times, the 4D seismic data seemed to be more noisy. In this case, it may thus be more difficult to interpret these data in terms of efficiency of the process, what emphasizes the importance of estimating the potential of 4D seismics prior to any acquisition.

Discussion
First of all, our integrated simulation approach helps to rank the impact of different parameters and effects (reservoir heterogeneity, fluid substitution, and pressure effects) and their unintuitive combined effects on seismic response. Regarding the contribution of time-lapse seismics to monitor chemical EOR, our integrated experiment once again clearly illustrates the evolution of 4D seismics. More precisely, in the 1990s the technique was considered as a sophisticated geophysical interpretation tool and was used by a handful of specialists. Now the technique is commonly used and helps in practical reservoir management. First of all, we emphasize the necessity of clean and appropriate processing of the data, as recommended in classical references on the topic (for instance Johnston, 2013). All these requirements were met in our study since we were dealing with direct problem simulation in perfect seismic acquisition and processing conditions. We are perfectly aware of the importance of correcting, or compensating for the variations in acquisition (e.g., overburden issues, seasonal effect or more general external condition variations not due to the subsurface exploitation, source/receiver positioning. . .), processing parameters, or the background noise as detailed in the last reference.
No noise equivalent to acquisition noise, background noise. . ., contributing to non-repeatability of the seismic acquisition, was numerically added to our synthetics. However noise is "intrinsically" present in the synthetics due to the randomness of the media induced by the geostatistical simulation characterized by Table A1. This also induces randomness of fluid and pressure distributions. All of these factors add noise to the synthetic data. Variations of both fluid saturation and fluid pressure due to reservoir exploitation are illustrated by Figure 9. In each of the two considered cases (Water flooding [WF] and EOR injection [SP]), one can roughly separate two spatial zones. A first zone bounded by a pseudo elliptical frontier plotted in blue on the figure where variations of both fluid saturation and fluid pressure are observed. We coarsely call this region the "swept zone". Outside this zone, even far from the boundary of the swept zone, mainly variation of fluid pressure is observed. The detection of the frontier between these two zones by 4D seismics is central in the present paper. According to Figure 14 the maximum value of the 4D NRMS due to 4D effect in the "swept zone" is around 20% both for Water Flood (WF) and EOR injection (SP). This quantity is called by some authors (e.g., Johnston, 2013) the detectability, that is to say the magnitude of the seismic response to production changes in the swept zone. Furthermore, Lumley et al. (1997) detailed the reservoir, fluid and seismic properties that affect the detectability factor and influence the chance of success of 4D seismics. In our case the observed detectability value is consistent with typical values of 4D NRMS for offshore survey reported in the literature, varying from 10% to 25% (e.g., Kragh and Christie, 2002). Note exceptional values as high as 35% in Draugen field measured by Koster et al. (2000) for water replacing oil. In contrast the 4D NRMS noise away from the swept zone is substantially smaller than the NRMS in the swept zone. A stronger 4D NRMS noise outside the swept zone is observed in the WF case (roughly 8% in average) than in the SP case (roughly 5% in average), both spatial dependent. In both cases, the 4D NRMS noise is clearly weaker than the measured detectability of the swept zone 20% except in specific directions, roughly E-W direction, and only in the case of WF as illustrated by Figure 14. All these ensure the detectability of the boundary of the swept zone.
A more complete study of the effect of "extrinsically" added noise on the non-repeatability of seismics and its implication on 4D interpretation, such as those of Houck (2010), will be part of further studies.
With the proviso related to the absence of "extrinsic" noise, we demonstrated that 4D seismics can contribute to the assessment of chemical EOR efficiency. First, as in the case of other types of subsurface exploitation, one of the major contributions of the technique is to enable the identification of water front around the swept zone and also the identification of unswept zones and flow barriers. The example of the E-W channel to the south is convincing, with no substantial variations in fluid pressure and in fluid saturation, inducing negligible seismic difference response to the south of the channel, as illustrated by Figure 11. This information is of value to understand the reservoir flow behavior and the hydraulic connectivity, and can assist in improving reservoir management. However, due to the multiplicity of the factors and their combined effects, 4D seismics may fail to detect chemical impacts on some geological reservoirs or for some injection strategies. That is the reason why such kind of integrated workflow is very valuable to help to decide whether to monitor a field or not, and to design the most relevant monitoring strategy depending on the objective (fluid saturation quantification, fluid movement estimation, connectivity. . .). Frequently repeating seismic acquisition on few specific 2D lines may also be chosen thanks to such feasibility study. As proposed by Hatchell (2015) for deep-water fields, a lower cost solution for 4D monitoring could be to reduce monitored field surface and to shoot small target subsets taking benefit of the scalability of OBN solution. This would enable frequent monitoring of specific areas of interest which could be chosen thanks to this kind of integrated workflow as well.
These conclusions on the contribution of 4D seismics to the present work may appear rather less impressive compared to what has been achieved in other contexts with this technique, for instance on Sleipner field (North Sea) for the geological storage of carbon dioxide (CO 2 ) as described in Dubos-Sallée and Rasolofosaon (2010). In the present study, except for the high data quality, we faced multiple unfavorable conditions: (i) the reservoir rock was roughly less compliant in average, (ii) the saturating fluids were of comparable compressibility, (iii) fluid pressure effect could not be neglected, and (iv) the reservoir exhibited rather strong facies and permeability distribution heterogeneity.
Despite this, and with the proviso of the relevant seismic processing and choice of attributes, important geometrical details of the swept and unswept zones were obtained. To go beyond this mere geometrical approach, a first step would be to use seismic cubes at zero-offset but also at small-offset and at large-offset, or even both P-wave data and S-wave data, at least in order to discriminate between pressure effect and fluid substitution effect (e.g., Landrø, 2001). An even further step would be to use full waveform modeling instead of 1D-convolution, with the resulting further complication of the inversion and interpretation of the data (Bourgeois et al., 2004). However, such a more complete modeling is more realistic and can contain more information on both the 3D geometry and the elastic properties of the reservoir.
For real field cases, several reservoir models may reproduce the available data (logs, production data) but provide various production forecasts for a given development plan.
If available, this uncertainty should be taken into account in the proposed workflow and seismic data should also be used to constrain the model through history-matching.
Finally, our workflow can be of interest for communication between geophysicists and reservoir engineers involved in the same operational team and could facilitate a virtuous circle between them. With this workflow, reservoir engineers can test the optimal range of production conditions leading to the best oil recovery. Using this information, geophysicists can then infer, and share with reservoir engineers, in which production conditions quantitative 4D seismic interpretation would be optimal. Then reservoir engineers have all necessary information to determine the best production strategy.

Conclusion
This paper investigates the ability of 4D seismics to capture the chemical impacts during a tertiary oil recovery phase and to identify characteristic features that can be related to the success of the EOR strategy. Seismic cubes were computed on a realistic case based on an outcrop analog of turbiditic reservoirs, considering two production strategies: the first one consists in a surfactant-polymer injection, and the second one in a classical waterflooding. The computation of the seismic responses was performed based on a workflow that spans from the geological description to seismic monitoring simulation and seismic attributes analysis, through geological and reservoir modeling, flow simulation, and rock physics modeling. The surfactant-polymer injection strategy induces changes in the dynamic behavior of the reservoir that are visible on the seismic attributes. Despite the loss of resolution, it may also lead to characteristic features on the seismic response. In particular, a clear delineation of the waterfront compared to classical water flooding can appear. In this case, the real seismic cubes may be used to estimate the efficiency of the process. In other cases, the seismic response may however be too noisy to justify the acquisition process. This underlines the necessity to test the ability of 4D seismics to detect a chemical impact before monitoring a field.
With the workflow considered here, each expert of the different fields can benefit from the availability of a complete set of well-controlled data of various types to test and improve his own tools. In contrast, the non-experts can easily and quickly benefit from "hands-on" experiments to understand the involved phenomena. As a direct consequence, it is easy to understand that any output of our workflow can be used for educational purpose towards a well-documented full case study. The new digital tools and techniques provide a new relation to training and can be a very powerful way to increase the learning capacity of students or professionals. These even allow to better train them directly in the field and in the laboratory (e.g., Joseph, 2017;Marçal et al., 2017).
Regarding future work, the completeness of the workflow and the manageability of the modest-size model open the door to several extensions of the present work. A short-term future work could be the study of alternate production strategies and well locations on this synthetic case. Complementary tests have already been performed in such direction, but only limited illustrative examples have been reported in the present paper for brevity. A detailed description of complementary results can be found in Zarate-Rada (2016). The proposed workflow can also be directly applied to geological reservoirs all over the world. There is no methodological barrier that would avoid our approach to be adapted even to quite different geological contexts and fields.

Elements for facies modeling
Four vertical logs taken from the outcrop were used to model facies distribution, extracted from sedimentological sections 5, 8, 12, and 14 (pseudo-well locations displayed in Fig. A1).
As outcrop thickness layers (~0.15 m) are smaller than vertical mesh size discretization (~0.5 m), an upscaling step is used to import facies logs to the Ainsa reservoir grid (Fig. A2).
In addition, facies proportions are deduced from four vertical sedimentological sections selected across the outcrop and also used to build the facies model obtained by Jardin et al. (2010). To better constrain a new facies model, four synthetic vertical logs (Fig. A3) were extracted from the facies model simulated in Jardin et al. (2010). The approaches considered for simulation are the Truncated Gaussian method (Doligez et al., 1999, Galli et al., 1994 for units C1 and C2.2 and the sequential indicator simulation for units C2.1 and C3 (Fig. A1b). These geostatistical approaches have already been applied successfully in many depositional settings, including deep offshore settings (Jardin et al., 2010, Lerat et al., 2007.

Associated petrophysical parameters
Petrophysical properties are generated using a truncated Gaussian distribution for porosity and a truncated lognormal distribution for permeability. The characteristics of the property distributions depend on facies as detailed in Table A1.

Rock properties for flow simulation
We used the following values for average compressibility: 2.96 for the thick-bedded sandstone, 10 for the gravelly mudstone, 90 for heterolithics, 12.4 for mudstone-clast conglomerates, and 6.51 for conglomerates, all in units of 10 À5 bar À1 .

Fluid properties
The pressure dependence of the oil formation volume factor, which is the ratio of the oil volume at reservoir (in-situ) conditions to that at stock tank (surface) conditions, and the viscosity of the oil phase are plotted in Figure A4. The corresponding parameters for water are identical to that tabulated by Batzle and Wang (1992). Average oil density and water density are 853.7 kg/m 3 and 1012 kg/m 3 at surface conditions.
The Kr-Pc model Different numerical reservoir production modeling show that thick-bedded sandstone and heterolithics have a strong impact on flow profile in the reservoir. For these two facies  the Relative Permeability-Capillary Pressure (K r À P c ) curves are given in Figures A5a and A5b respectively. Both figures are organized following the same plot convention, with capillary pressure as a function of water saturation in the top subfigure and the dependence of Relative Permeability to Oil (KRO) (red curves) and Relative Permeability to Water (KRW) (blue curves) as functions of water saturation at the bottom. The saturation range for which the two fluids are mobile is limited by the Irreducible Water Saturation (SWI) and the Residual Oil Saturation (SORW) (e.g., Zinszner and Pellerin, 2007). The end-points of the K r curves in this saturation range are marked by black solid disks in both figures. For convenience, the curves are linearly extrapolated to the expected extreme bounds, namely (KRO = 1, S w = 0) and (KRO = 0, S w = 1) for the relative permeability to oil KRO, and (KRW = 0, S w = 0) and (KRW = 1, S w = 1) for the relative permeability to water KRW. These relative permeability and capillary pressure curves are typical of water wet rocks. For all facies, water saturation varies between 30% (SWI) and 68% (1 À SORW) except for sandstone whose saturation varies from 15% to 75.5%. Although extreme K r and P c values are the same for all facies, the broader range of water saturation variation for the thick-bedded sandstone makes it the best oil reservoir facies in terms of oil volume storage and facility to produce it (Tab. A2). We may also notice that flow velocity is controlled by both permeability and relative permeability. As a consequence, even if the K r À P c model is the same for some facies, flow velocity within these facies will be different due to the differences in the permeability distributions.

Appendix B Production strategy
The conditions of exploitation were chosen as realistic as possible. Target injection and production rates were imposed (Tab. B1) with pressure constraints at the injector well (maximum pressure) and producer wells (minimum pressure). As a consequence, if injector pressure increases to the maximum authorized pressure, the well bottom hole pressure is maintained to this limit pressure and the water injection rate decreases: the injector well behavior changes (a) (b) Fig. A5. Capillary pressure and Relative permeability curves for (a) thick-bedded sandstone facies and (b) heterolithics facies. SWI: irreducible water saturation; SORW: residual oil saturation to water; KROWM: relative permeability to oil at minimum water saturation; KRWSOR: relative permeability to water at residual oil saturation.  B1. Water injection rate taking into account max pressure limit for SP injection scenario in green and WF scenario in blue (Q max = 600 m 3 d À1 ).

Fig. B2
. Water injection rate taking into account max pressure limit for SP injection scenario and the two target injection rates Q max = 600 m 3 d À1 and Q max = 350 m 3 d À1 (respectively in continuous and dashed lines).
S o and S w : oil saturation and water saturation.

Pressure effect on elastic properties of rocks
Up to now we have ignored any pressure-dependence of any kind in the model. In fact some of the parameters introduced in equations (C.1)-(C.4) depend either on the confining pressure P c and/or the pore pressure P p : the fluid bulk modulus M (fluid) depends on P p (Batzle and Wang 1992), while the bulk and shear moduli of the dry rock M (dry) and l (dry) both depend on P c and P p . However, it can be shown experimentally (e.g., Zinszner, 2009, 2012) that these moduli essentially depend on the differential pressure P diff (P diff = P c -P p ) and not in an independent way on P c and on P p . We use the simplest model of pressure dependence assuming a power law for the dependence of M (dry) and l (dry) with the differential pressure P diff : with P diff1 and P diff2 corresponding to two differential pressures and h M and h l being the Hertz exponents of the bulk and shear moduli (experimental measurement of Hertz exponents in various types of rocks are reported by Rasolofosaon and Zinszner, 2012). Note that we neglect any pressure dependence of the rock porosity / and of the bulk modulus of the grain constituent M (grain) .

P-and S-wave velocities and impedances of fluid-saturated rocks
Then P-and S-wave velocities V P (sat) and V S (sat) , and impedances I P (sat) and I S (sat) of the fluid-saturated rock are given by: where the density q (sat) of the fluid-saturated rock is given by: with q (oil) , q (water) , and q (grain) respectively designate the densities of oil, water, and grain constituent.

Application to our synthetic case study
For our monitoring feasibility study, we associated to each facies the average petroelastic parameters listed in Table C1. Using equations (C.1)-(C.8) we were able to infer the impact of fluid substitution and pressure changes on the behavior of facies elastic moduli. First of all, Figure C1 gives the pressure dependence of the density and the bulk moduli of oil and water. As expected, water and oil exhibit a slight density increase with pressure, oil being less dense than water. Note the comparable value of the density versus pressure gradient. Similar observations can be made for the bulk modulus with the proviso that the modulus versus pressure gradient is twice larger for oil than for water: roughly Symbol M designates the bulk modulus, l is the shear modulus, q is the density, m is the Poisson's ratio; h M is the Hertz pressure-dependence coefficient for M, h l is the Hertz pressure-dependence coefficient for l. Grain and dry respectively stand for the grain constituent and for the rock skeleton. TkS: thick-bedded sandstone; gM: gravelly mudstone; H: heterolithics; McC: mudstone-clast conglomerates; C: conglomerates.
14 Â 10 À4 GPa/bar for oil against 7 Â 10 À4 GPa/bar for water, oil being expectedly more compressible than water. The systematic increase of fluid bulk modulus with fluid pressure is one of the numerous manifestations of nonlinear elastic behavior of fluids (e.g., Beyer, 1965). The fluid pressure dependence of the shear modulus of thick-bedded sandstone and heterolithics is illustrated by Figure C2. As expected the result is independent of the saturating fluid since fluids cannot support shearing stress, as long as the fluid viscosity is not too strong, typically smaller than 10 3 centipoise as demonstrated experimentally by Rasolofosaon and Zinszner (2012).
As widely known, an increase of the fluid pressure P p , implying a decrease of the differential pressure P diff at fixed confining pressure P c , induces a softening of the rock (decrease of the shear modulus shown on Fig. C2) by opening the cracks and compliant pores (e.g., Bourbié et al., 1987;Mavko et al., 1998). Note that the shear modulus vs pressure gradient is roughly four times larger for the thick-bedded sandstone (23 Â 10 À4 GPa/bar) than for the heterolithics (5 Â 10 À4 GPa/bar). Figure C3 shows the fluid pressure dependence of bulk moduli of fluid-saturated thick-bedded sandstone and heterolithics, for different levels of oil saturation S o . The results are less obvious for the bulk modulus than for the shear modulus. More precisely, in the case of the thick-bedded sandstone a slight increasing trend of the bulk modulus with fluid pressure is observed at large oil saturation (typically for S o > 0.75), whereas a slight decreasing trend