© D. Do Couto et al., Published by EDP Sciences 2024

1 Introduction

In mainland France, the middle Rhône valley hosts a large number of energy-producing industries and high-risk industries (nuclear power stations, hydroelectric dams, chemical factories, etc.). The Tricastin region in the middle Rhône valley has been industrialised since the 1950s following the construction of dam and major hydraulic works by the Compagnie Nationale du Rhône (Donzère-Mondragon channel), then by nuclear fuel cycle facilities in the 1960s, and finally the Tricastin nuclear power plant (cooled by channel water) in the 1970s. Even though south-east France is known to be impacted by low to moderate seismicity (Ritz et al., 2021), rare yet powerful earthquakes could occur in the study area, as shown by the Le Teil earthquake (∼5 Mw) and shallow earthquake swarms in Tricastin (∼4-4.5 Mw; e.g. 1773, 1873, 1936; Thouvenot et al., 2009; Manchuel et al., 2018; Bollinger et al., 2021). In addition, the region shows ancient surface ruptures (paleoseismic) identified on the Nîmes Fault, the Ventoux Lure thrust and the La Rouvière Fault (Ritz et al., 2021, Bellier et al., 2021), evidence for its past and probably significant seismic activity (Fig. 1).

The Le Teil earthquake, which occurred in November 2019 on the La Rouvière Fault (belonging to the Cévennes Fault system) just north of the Tricastin area (Ritz et al., 2020, 2021) shed light on the seismic hazard due to seismic amplification because of the presence of soft sediments in the Rhône valley (Gélis et al., 2022), just below the nuclear facilities. The soft sediments are mostly made of Pliocene deposits infilling a canyon carved by the paleo-Rhône during the Messinian Salinity Crisis (Hsü et al., 1973; Clauzon, 1982), a recent (<6 Ma) environmental crisis in which the Mediterranean Sea level fell hundreds of metres below its present-day level. During this short episode (<600 kyr) the paleo-Rhône incised a canyon, but the present-day morphology in the area remains elusive as no geophysical subsurface imaging enabled the production of a comprehensive 3D architecture. Modelling of site effects (seismic amplification) requires a 3D image of the canyon that is as accurate as possible, as well as information concerning the nature and thickness of the sedimentary fill (Froment et al., 2022).

In the middle Rhône valley, the Tricastin area is also known to be the locus of a slope break (known as a knickpoint) in the course of this Messinian Rhône canyon (Clauzon, 1982) that is hypothesised to mark a stagnation phase during the sea-level drop caused by the Messinian Salinity Crisis (Beaudoin et al., 1997). Evidence of this knickpoint in the paleo-river profile was mainly based on the Pierrelatte borehole that cross cut the base of the Pliocene (Ballesio, 1972). Although erosion modelling along the whole Messinian paleo-Rhône reproduced a slope break in the same area, it did not reproduce such an abrupt slope break but estimated a rather smoother transition from the upstream paleo-Rhône to the downstream part (Gargani, 2004a). Detailed analysis of seismic profiles would greatly improve our knowledge concerning the existence of this knickpoint.

In this study, we used surface data (geological maps, digital elevation model) and subsurface data (seismic profiles, borehole data) to conduct a thorough analysis and to produce a 3D reconstruction of the canyon morphology of the Messinian paleo-Rhône and its SW tributaries (paleo-Ardèche and paleo-Cèze rivers). This paper considerably improves our knowledge of the processes by which the canyon was dug and then filled in, in relation to the different stages of the Messinian salinity crisis. Although to a lesser extent, our study also advances our knowledge of seismic hazards, e.g. amplification of the seismic ground motion linked to the specific geological context, identification of any active faults that could affect the plio-quaternary filling of the canyon, which is essential in a context where critical industrial sites are present.

Here, we present the results of our analysis of high-resolution seismic profiles of the central part of the middle Rhône valley. Analysis of the seismic profiles revealed details of the Messinian erosion surface and the Pliocene filling of the Rhône, Ardèche and Cèze canyons. As part of this analysis, we also wanted to link the stages of sub-aquatic filling, which are clearly visible on the profiles, to those of the Messinian-Pliocene evolution of the aggradation of the paleo-river/river. Finally, we consider the possible connection between the Ardèche Messinian canyon and its Messinian deep karst.

The rest of the paper is organised as follows: after presenting a synthesis of the regional geological framework and describing the data used in this study, we conducted a cross-correlation of wells by revisiting data collected from old and recent boreholes we then used together with seismic interpretation to unravel the history of the Messinian-Pliocene history and to model the MES. Finally, we discuss the results of this work with respect to the course and morphology of the Rhône and Ardèche canyons, as well as their infill and structural geology.

Fig. 1 (A) Synthetic structural map of the study area (black box) showing the main tectonic domains mentioned in the text and the course of the Rhône and Ardèche rivers in dotted blue lines. (B) Seismotectonic map of the lower Rhône Valley between the Cévennes Fault and Nîmes Fault. Plotted historical earthquakes are from SisFrance Database (www.sisfrance.net: data from SisFrance, BRGM, EDF, IRSN, 2022), magnitudes from FCAT (Mw catalog − Manchuel et al., 2018). BDFA refers to the Potentially Active Fault Database (Jomard et al., 2017), paleoseismological works reported in Baize et al. (2002) and Bellier et al. (2021).

2 Geological framework

2.1 Tectonic and stratigraphic framework

The study area is located in mainland France in the middle Rhône valley (Fig. 1). It lies between the Massif Central to the west and the Alps to the east. The Rhône flows from the inner Alps, through the Jura fold belt and finally into the Tertiary basins formed during the extensions linked to the opening of the West European Rift and the Gulf of Lion (Eocene, Oligocene). To the south, the Rhône river flows into the Mediterranean Sea.

The Rhône plain is currently covered by Quaternary sediments (in grey in Fig. 2). The right bank of the Rhône mostly is composed of Lower and Upper Cretaceous series, while Miocene series cover Saint-Restitut hill on the left bank of the river (Fig. 2). Pliocene deposits (in pink in Fig. 2) are found on both sides of the Rhône river, mostly nested within Cretaceous series. After the Mesozoic Tethysian rifting responsible for the inception of the major NE trending faults (Bonijoly et al., 1996), the study area experienced three main tectonic stages, namely (i) the Pyrenean shortening from Paleocene to early Oligocene, (ii) the Oligocene-Early Miocene extension and (iii) the Miocene Alpine shortening.

The convergence between the Iberian and Eurasian plate from the Paleocene to the early Oligocene caused E-W oriented folds sometimes combined into N-verging thrusts that affected the Mesozoic series (Arthaud and Laurent, 1995; Rangin et al., 2010) as well as the reactivation of pre-existing NE-striking faults in a left-lateral movements (Arthaud and Matte, 1975) such as the Nîmes fault. The most prominent structures in the study area that resulted from this event are the Mondragon and Orange anticlines to the southeast, the Echavarelles anticline in the centre of the study area, and the Donzère anticline to the north (Fig. 2). From the Late Oligocene to the Early Miocene, the study area experienced two successive E-W extensions related to the West-European rifting, followed by a NW-SE extension concomitant with the opening of the Gulf of Lion and the rotation of the Corsica-Sardinia blocks (Séranne et al., 1995, 2021; Lacombe and Jolivet, 2005). These riftings resulted in the inception or re-activation of NE-trending normal faults, most of which occurred south-east of the study area (Ballas et al., 2014). It has been reported that the Saint-Remèze fault system may have been re-activated at that time, causing normal displacement along the fault (Pascal et al., 1989). More recent Alpine orogeny had little impact on the study area (Champion et al., 2000) as it mostly accommodated the deformation along major accidents further south, such as the Luberon thrust (Rangin et al., 2010; Clauzon et al., 2011), and along the northern part of NE-SW Cevennes-fault system (Bellier and Vergely, 1987). Vertical movements are also documented. The alpine forebulge uplift (Besson, 2005) was mainly active during the Lower and Middle Miocene and partially explains how in some reliefs, for example, in the Saint Restitut Hill, Upper Miocene conglomerates remain over Burdigalian strata. Based on numerical modelling of isostatic responses of the Rhône valley during the MSC, Gargani (2004a) evaluated a positive rebound of between 20 and 180 m, in agreement with the results of field studies (Mocochain et al., 2006b). Other post-alpine movements created a slight uplift and tilting of the Miocene, Pliocene and even Quaternary series have already been reported by other authors (Guérin, 1973; Besson, 2005). This deformation results from post alpine isostatic adjustment (due to erosion unloading; Malcles et al., 2020) or is linked to the uplift of the neighbouring Massif Central (Olivetti et al., 2016).

The general structural trend in the study area is characterised by an overall gentle monoclinal trend, slowly dipping east-south-eastwards. However, within the Rhône valley, there is an outcrop of Barremian limestone of Urgonian facies in the city of Pierrelatte. The Urgonian is also recorded at a depth of 339 m (−283 m b.s.l.) in the Pierrelatte well. To explain the rapid offset of Urgonian rocks, Debelmas et al. (2004) suggested the occurrence of downfaulting across a NE-trending structure known as the Pierrelatte fault. So far, data proving its existence at depth and its geometry are lacking.

Figure 3 shows a stratigraphic column of the whole study area based on Champenois et al. (1971), Masse et al. (1980), Pascal et al. (1989) and Debelmas et al. (2004). It starts from the Hauterivian marl and limestone that can vary in thickness from 50 to 750 m (Pascal et al., 1989). As detailed later, this formation was not distinguished on the seismic reflection image, so only its top is indicated on the lithostratigraphic log. The most prominent and thickest carbonate strata outcrop in the Saint Remèze plateau and along the Ardèche Gorges (Ferry and Rubino, 1989): they are made of ammonite- and rudist-rich limestone of Barremian-Lower Aptian age known as the Urgonian facies, and are approximately 550-600 m thick. Locally, the thickness of the Barremian can exceed 700 m, as recorded in the Mondragon well considering the local dip of layers (BSS002CLTG in the borehole subsurface database BSS from BRGM: https://infoterre.brgm.fr). The Urgonian terminates by a mid-Aptian marl and limestone succession known locally as the Bedoulian facies (Fig. 3) belonging to the Chabert formation (Pictet and Delanoye, 2017). The Aptian ends with ∼100 m of dark grey marls and glauconite bearing sandy limestone, known locally as the Gargasian-Clansayesian sandstones facies. The Albian mainly comprises sandstone/calcareous sandstones and marls, about 80-90 m thick, topped by the uppermost Albian Marls (Vraconian flooding). The Cenomanian evolves from lignites and sandstones (quartzarenites) to the south and southeast of the study area (Uchaux area; Malartre, 1994), towards marly limestone intercalated with calcareous sandstones in the northeast (Clansayes area, Fig. 2). The stratigraphic column continues with ∼200 to 250-m thick glauconitic sandstones, as well as sandy to micritic limestones Turonian in age, outcropping in the Mondragon anticline (Fig. 2). The Upper Cretaceous then evolves into Coniacian sandy limestones, whose thickness ranges from 30 to ∼150 m and grades into Santonian sandstones with rudist intercalations.

The Upper Cretaceous is unconformably overlain by Eocene and then Oligocene deposits. The Eocene is marked by red and white lateritic sands and quartz sandstones, often found as lenses tens of metres wide. At its base, the Oligocene is also marked by an unconformity, which mainly outcrops on the northeastern side of the study area, and mostly comprises lacustrine deposits, sometimes brackish, composed of marls and limestones evolving into red marls towards the top of the series. As a whole, the Oligocene is ∼50 to 70 m thick and dips slightly towards the E-SE.

The Miocene is unconformably deposited above previous series. In our study area, the Miocene mainly consists of marine Burdigalian strata in Saint-Restitut hill, made of approximately 120 m thick molasses (conglomerates, sandstones evolving to marls and carbonates, partly bioclastic at the top of the series; Lesueur et al., 1990). It outcrops at an altitude of approximately 180 m on the western side of Saint-Restitut hill, with a mean 3° dip towards the E-SE. Younger Miocene series, Langhian to Messinian in age, outcrop further east of the study area in the Miocene Valréas Basin (Demarcq, 1970; Rubino et al., 1990), as well as further south fringing the Mondragon anticline (Fig. 2).

The Pliocene is unconformably nested within earlier series, mostly in contact with Cretaceous strata (Fig. 2). Because Pliocene series infill Messinian valleys and canyons, its thickness varies markedly from one place to another. Thanks to descriptions of both outcrops and boreholes like the Pierrelatte well (Fig. 4), the Pliocene in the Tricastin area is known to be made up of four subunits (Ballesio, 1972), from the base to the top:

• basal pebble-sized conglomerates reworking lower Miocene material, and lenticular breccia bodies lying on top of slopes reworking hillside strata such as Miocene deposits; the thickness of this subunit can reach 50 m locally;

• continental fluvio-lacustrine deposits made of variegated marls followed by azoic greyish to yellowish clays/marls and topped by brackish layers containing bivalves (such as in Saint-Marcel d’Ardèche);

• marine bluish to greyish marls/clays with a small proportion of sand, locally containing coastal fauna (Cerithium and Pectinidae species);

• fluvio-lacustrine sands, conglomerates and marls, exhibiting some lignite layers and containing freshwater fauna.

On top of this succession lie unconformably Quaternary terraces (Mandier, 1989), however these are not dated.

Fig. 2 Simplified geological map of the studied area (modified from Champenois et al., 1971; Pascal et al., 1989; Masse et al., 1980; Debelmas et al., 2004) showing location of seismic profiles used in this study and the main boreholes mentioned in the text.

Fig. 3 Synthetic stratigraphic log of the study area, from the Hauterivian up to Quaternary. Seismic units interpreted on seismic profiles are reported on the right side of the log.

Fig. 4 Cross-correlations of the main boreholes intersecting the base of the Messinian erosional surface in the Tricastin area containing lithological data.

2.2 The Messinian Salinity Crisis in the Rhône Valley

The Messinian Salinity Crisis is an outstanding paleo-environmental event characterised by a two-stage sea-level fall (Clauzon et al., 1996) followed by a two-stage sea-level rise (Bache et al., 2012; review in Roveri et al., 2014). The first sea-level drop, occurring at 5.97 Ma (Manzi et al., 2013) was moderate in amplitude (150-200 m in Clauzon et al., 2015b), it was responsible for the precipitation of evaporites in the peripheral shallow basin around the Mediterranean Sea. At the end of the Messinian period, at 5.6 Ma (CIESM, 2008), a huge drop in the sea level of the Mediterranean Sea estimated at between 800 and 1500 m (Bache et al., 2009, Urgeles et al., 2011) led to the incision of deep canyons such as the Rhône canyon, down to 1300 m below present-day sea level (Guennoc et al., 2000) and to more than 1,000 m below the Camargue (Clauzon, 1982; Clauzon et al., 1992). Reflooding following the drawdown took place in two stages, with a first slow sea-level rise, followed at 5.46 Ma by a rapid sea-level rise (Bache et al., 2012). Sea-level oscillations during the MSC also led to the development of a karstic network along the Ardèche river (Mocochain et al., 2006b, 2009).

The incision resulting from the fall in the level of the Mediterranean Sea is highlighted by the deeply unconformable nature of the Pliocene series over all the preceding formations. Along the Rhône valley, this specific unconformity has been known since the pioneer work of Fontannes (1882), see Clauzon (1974) for a historical review. Previously, the estimated incision depth in our study area was based on (1) the identification of a pre-evaporitic abandonment surface (i.e. the last and topmost gently dipping surface of the upper Miocene alluvial plain), which recorded the position of the Rhône just before the MSC (i.e. at an altitude of 360 m on the Saint-Remèze plateau; Clauzon, 1982, Martini, 2005, 2022), and (2) the interpretation of the Pierrelatte well (Demarcq, 1960) suggesting that the base of Pliocene was located at -236 m below present sea-level (b.s.l). These observations allowed Clauzon (1982) to estimate the depth of incision at 580 m in the Tricastin area. Clauzon (1982) also drew the Messinian Rhône canyon profile, based on a compilation of wells that crossed the base of Pliocene deposits, highlighting the occurrence of a knickpoint south of the Tricastin area. North of this knickpoint, the Messinian Rhône appeared to have an average slope of ∼1.8° (∼3%), which increased southwards with an average slope of ∼9° (∼16%).

The question of the junction between the Ardèche valley and the Messinian Rhône canyon was enigmatic due to the lack of any way of estimating the depth of incision of the tributary. This uncertainty left room for numerous interpretations concerning the hypothetical digging of a Messinian canyon by the Ardèche river. Some authors considered that the Ardèche had not incised its valley or had only very slightly incised it during the Messinian (Sadier, 2013; Martini, 2019) while others did not separate the river from its base level, which is the Rhône Messinian canyon (Baulig, 1953; Belleville, 1985; Mocochain et al., 2006a, 2006b, 2009; Tassy et al., 2022).

2.3 Regional uplift since the Messinian

It should be borne in mind that the current surface position of the markers derived from the MSC is not that at their original altitudes. For example, the current transition from marine to continental sedimentation in the Pliocene is observed at an altitude of ∼130 m in several places surrounding the Middle to Lower Rhône valley (Ballesio, 1972; Mocochain et al., 2006a, 2006b, 2009; Tassy et al., 2022) and was confirmed by our own field study. According to eustatic seal-level charts (Lisiecki and Raymo, 2005; Miller et al., 2020), this marker should be located between 0 and 40 m above present sea level, indicating a regional uplift of around 90-130 m that occurred later than the date of emplacement of this marker (estimated at ∼4.7 Ma − Mocochain et al., 2009) as already proposed by Denizot (1952) and Mandier (1988). This uplift could be linked to the isostatic rebound following glacial erosion of the Alps and/or the associated mantle upwelling, as well as to Pliocene volcanism, as proposed by Olivetti et al. (2016) who suggest an uplift of around 200-300 m in the French Massif Central. Gargani (2004b) estimated a post-Messinian uplift due to an isostatic rebound of 20-170 m in the Tricastin area. Malcles et al. (2020) also suggest that part of this uplift in the Cevennes area for the last 4 Ma, is due to isostatic readjustment caused by erosion.

3 Dataset and methodology

The dataset used for this study, which is composed of 11 seismic profiles, has been the property of IRSN since 2020 (Fig. 3), it comprises nine multi-trace onshore profiles (P1, P2, P3, P4, P5, P6, P7, P9 and P10) acquired with vibroseis and geophones; plus two marine-type seismic profiles (P8C16 and P8D16) acquired with an airgun system and a streamer on the Donzère-Mondragon channel. Onshore profiles were recorded up to 3s in two-way travel time (TWT), the two channel seismic profiles up to 2s TWT. A complementary seismic profile (94 MAR01) was also included in this study. The complementary profile is one of the profiles acquired by ANDRA in the 1990s and reprocessed more recently (Hollender et al., 2015a; Hollender pers. Comm.). As we tried to unravel the recent history of the Rhône valley, we focused on the upper part of the seismic profiles, up to 0.75s twt, where the seismic signal is well recorded. We used the standard seismic-stratigraphic analysis principles presented by Mitchum and Vail (1977) to characterise and define the seismic unit architecture together with the geometry of their bounding surfaces (onlap, toplap, downlap and truncations) and vertical stacking patterns. The seismic units defined on the seismic profiles were then linked to previous studies and compared to those in a deep well at Codolet (BSS002CMTA) located farther south (Ferry et al., 1997) near the 94 MAR01 profile (Fig. 3). We performed the seismic stratigraphic analysis on all the twt seismic profiles of the Tricastin area and afterwards depth-migrated the various horizons based on the P-wave velocities retrieved from the depth-converted P4 profile (see explanation below) which we considered to be representative of the local geology of all seismic horizons that could exist in this area.

Depth conversion of seismic units requires the use of velocity data. Seismic reprocessing of seismic profile P4, intersecting two canyons and featuring an intermediate island formed by Cretaceous-age series, as demonstrated by the Grand-Malijac, Mirabelle and Préférence boreholes (see Figs. 4 and 5), was pivotal. We used an original methodology (following techniques used and described in Beccaletto et al., 2011) aiming at specifying the shape of the Ardèche and Rhône paleo-canyons, the P-wave velocities of the Messinian canyon filling units and the enclosing formations. This was achieved using a 2-step process (velocity, deep static guiding) detailed in Supplementary material. As a result of this conversion, the mean P-wave velocity of the Pliocene filling is ∼2200 m/s, in line with velocity of the Pliocene deposits further south (Schlupp et al., 2001), in the vicinity of Roquemaure-Pujaut.

Fig. 5 Seismic profile P4, intersecting the Rhône canyon and a tributary of the Ardèche canyon, in two-way travel time (twt) on top, and converted in depth at the bottom. The Grand Malijac borehole is projected on the depth converted seismic profile for reference. Vp velocities obtained from the depth-processing are drawn on the section for each formation. The MES is drawn in dotted red line. See Figures 4 and 5 for stratigraphic legend and text for more details. RC: Rhône canyon; ATrC: Ardèche tributary canyon.

4 Results, interpretation and MES surface model

4.1 Well correlation

A stratigraphic correlation of six selected boreholes, from Pierrelatte (BSS002BNBB) to Codolet (BSS002CMTA), is given in Figure 4. The deepest well in the study area is the Mondragon well (drilled to 1838 m). It reached Hauterivian marls and crossed a major thrust below which more than 1000 m of Urgonian facies were revealed. The other wells drilled through sedimentary sequences ranging from 70 m (Grand Malijac; BSS002BMYX) up to 834 m in thickness (Saint-Paul-Trois-Châteaux; BSS002BNWH). Both the Pierrelatte and Saint-Paul-Trois-Châteaux wells reached the Urgonian facies (Barremian in age), while those at Grand Malijac and Saint Paulet de Caisson (BSS002CLEB) ended in the Cenomanian sandstones. Above the Cretaceous series, only Pliocene strata were found, the base of which marked the Messinian Erosional Surface (MES). The thickness of the Pliocene deposits varies considerably depending on whether the well is located close to or far from a Messinian talweg axis. The Pierrelatte well revealed the occurrence of 269 m of Pliocene strata comprised of 16 m of basal conglomerates followed by 100 m of continental fluvio-lacustrine deposits, 130 m of azoic grey clays and marls and finally 23 m of marine bluish marls/clays (Ballesio, 1972). Until now, this was the deepest base of Pliocene known in the study area and was thought to be relatively close to the Messinian talweg (Clauzon, 1982; Debelmas et al., 2004). The Saint-Paul-Trois-Châteaux well, drilled for hydrological purposes in 2003, revealed (from bottom to top): 264 m of Barremian-Aptian limestones, 108 m of late Aptian (Gargasian) dark-grey marls and clays, 27 m of azoic greyish sandy clays, 75 m of azoic grey clayey sands, 350 m of Pliocene dark-grey clayey marls and finally 10 m of alluvium. According to the succession in this borehole, the base of the Pliocene is at -302 m b.s.l., 89 m lower than in the Pierrelatte borehole. It also means the Saint-Paul-Trois-Châteaux well is located closer to the talweg of the Messinian Rhône than the Pierrelatte well. The 102 m azoic interval made of sands and clays was not assigned to a specific interval (between 360 m and 462 m in depth), because of the lack of biostratigraphic markers, hence it could belong to either the Cretaceous series (Turonian sandstones) or the Pliocene. The interpretation of seismic profile P8D16 will play a key role in overcoming this uncertainty.

Further southwest, the Grand Malijac well crossed a very thin (∼3 m) Pliocene strata made of blue marls, whose base is at an elevation of +30.7 m a.s.l. Although located close to Grand Malijac, neither the Preference nor the Mirabelle (BSS002BNPG; Fig. 4) wells record any Pliocene. The area surrounding the village of Lapalud suggests a paleo-interfluvial area. Further south, the Saint Paulet de Caisson well, located on the present-day Ardèche river, revealed the occurrence of 142 m of Pliocene grey marls with a few sandy intercalations. Here, the MES stands at -100m b.s.l. Ninety-three metres of Pliocene marls were recorded in the Mondragon well, the basal 29 m being the richest in sands and lenses of lignite. The Codolet well is located further south of the Mondragon anticline, which crosses more than 443 m of Pliocene strata without reaching the MES. Four Pliocene units, hereafter named Pl1 to Pl4 from bottom to top, were recorded in the Codolet well: Pl1 comprises 94 m of a sandy unit with some intercalations of lignite-bearing blue clays, Pl2 is a 100 m thick unit richer in clay, Pl3 follows with a 127 m thick unit mainly comprised of clays but still fewer sandy beds, and finally the Pl4 unit made of blue marls and clays (Ferry et al., 1997). Based on the succession of sedimentary facies described in Codolet, we tentatively propose a correlation of Pl1 to Pl4 units in the Tricastin area.

4.2 Seismic units

Each seismic unit was first described using seismic profile P4, which is available in twt and depth (Fig. 5) and then expanded to the other profiles and correlated with the well records. Seismic interpretation was performed using IHS Kingdom suite seismic interpretation software. Seismic units are grouped with respect to their acoustic signal (reflector amplitude and continuity).

4.2.1 K1, K2 and K3 seismic units

At the bottom of the P4 seismic profile is a seismic unit named K1, marked by rather transparent facies, evolving towards the top to higher amplitude and more continuous reflectors, and is topped by the most prominent feature made of three couplets of very high amplitude continuous reflectors. The base of K1 unit was not interpreted in the seismic profiles. Aside from the base the K1 seismic unit dips slightly eastwards in all the E-W seismic profiles, and is close to the surface at the western end of profiles P1, P2 and P3 (Fig. 6). The depth converted model of profile P4 (Fig. 5) shows a deep domain corresponding to the K1 seismic unit with fairly high velocities (between 3,500 m/s and 4,500 m/s, increasing downwards). This K1 seismic unit almost certainly corresponds to the Lower Cretaceous level in Urgonian facies outcropping in Pierrelatte and close to the surface in profiles P6 and P7 (Fig. 7). In addition, along the P8D16 seismic profile, it was crossed by both the Pierrelatte and Saint-Paul-Trois-Châteaux wells (Figs. 8 and 9).

The K2 seismic unit, visible in the centre of the P4 seismic profile, is composed of lower amplitude yet relatively continuous reflectors. The top of the K2 seismic unit is marked by a strong couplet of reflectors. In more detail, the K2 seismic unit is made of two subunits of equal thickness in P4, the lower one characterised by lower amplitude seismic facies. The K2 seismic unit tends to become thinner towards the north of the study area, decreasing from ∼200 ms in P4 to ∼100 ms in profile P1 (Figs. 5 and 6) and thickens towards the south (profile P7; Fig. 7). Like the K1 unit, K2 dips eastwards. The depth converted model of profile P4 shows mean interval velocities of 2550 m/s in unit K2. This interval, crossed by many wells (Fig. 4) is made of Aptian, Albian and Cenomanian marls, sandstones and limestones. The low-amplitude reflectors observed on the lower part of unit K2 (Figs. 5 and 6) could correspond to grey Aptian marls (Gargasian level).

The overlying K3 unit is characterized by discontinuous, low- to moderate-amplitude reflectors, although it displays occasional high-amplitude reflectors locally. The thickness of this unit varies considerably as it is truncated by the Messinian erosion. For instance, the unit is missing in the Pierrelatte, Saint-Paul-Trois-Chateaux and Saint-Paulet-de-Caisson wells (Fig. 4). K3 is missing in most boreholes, except on (1) Grand Malijac where it corresponds to Turonian sandstone and (2) in the Mondragon well, which exhibits a thick sequence made of Albian calcareous marls and coarse sandstones evolving to sandy-clayey limestones of the Cenomanian. The K3 seismic unit is generally better preserved from erosion toward the south of the study area (see profile P6 and P7, Fig. 7) where its internal seismic facies feature high-amplitude reflectors toward the top. The depth converted model of seismic profile P4 shows that the K3 unit is characterised by an interval velocity of ∼2200 m/s.

Fig. 6 Seismic profiles P1, P2, P3 and P5 interpreted in twt. Faults are drawn in solid black lines. The MES is drawn in dotted red line. See Figures 4 and 5 for stratigraphic legend and text for more details. PF: Pierrelatte Fault; RC: Rhône canyon; AC: Ardèche canyon; ATrC: Ardèche tributary canyon.

Fig. 7 Seismic profiles P6 and P7 interpreted in twt. Faults are drawn in solid black lines. The MES is drawn in dotted red line. See Figures 4 and 5 for stratigraphic legend and text for more details. PF: Pierrelatte Fault. AC: Ardèche canyon; ATrC: Ardèche tributary canyon.

Fig. 8 Seismic profiles P8D16 and P8C16 interpreted in twt. Pierrelatte and Mondragon wells are projected on the seismic profiles. Faults are drawn in solid black lines. The MES is drawn in dotted red line. See Figures 4 and 5 for stratigraphic legend and text for more details. PF: Pierrelatte Fault; RC: Rhône canyon.

Fig. 9 Comparison of lithologies, Gamma-Ray (GR) signals and seismic facies of the Pliocene filling units Pl1 to Pl4 within the Saint-Paul-Trois-Châteaux and Codolet wells. GR curves are redrawn from printed electric logs. Seismic facies SF1 to SF5 are detailed in the text. Profile P8D16 is shown in Figure 8 and profile 94MAR01 in Figure 9.

Fig. 10 Seismic profiles P9, P10 and 94MAR01 interpreted in twt. The MES is drawn in dotted red line. See Figures 4 and 5 for stratigraphic legend and text for more details. RC: Rhône canyon; CC: Cèze canyon; TC: Tave canyon.

4.2.2 Messinian Erosional Surface (MES)

In the Tricastin area, the MES most often cuts through K2 and K3 units, and reaches the K1 unit in the southern part of the region (Figs. 5, 6, 8). Three main incisions were identified and correlated across seismic profiles: the Rhône canyon (RC), the Ardèche canyon (AC) and a tributary of the Ardèche (ATrC) (Figs. 5–8 and 10). Along the E-W seismic profiles, the main and simultaneously the deepest incision is located at the eastern edge of the present-day Quaternary valley, west of the Donzère-Mondragon channel. From the P2 seismic profile to the P4 seismic profile, the acoustic depth of the deepest incision deepens southwards from ∼500 ms to ∼620 ms (Figs. 5, 6 and 10). This incision is, without doubt, the Messinian Rhône canyon (RC in Figs. 5–8 and 10).

Unlike the other profiles, the P5 seismic line shows two adjacent canyons (Fig. 6): a western, narrow, 500 ms-deep canyon and a shallower one (300 ms deep). The most westerly canyon is also visible on the southern part of the P7 seismic profile, at a depth of ∼425 ms (Fig. 7). Being located at the present-day exit of the Gorges of the Ardèche, this canyon is very probably the Messinian Ardèche canyon (AC in Figs. 6 and 7). The other shallower canyon can be interpreted as the course of a tributary of the paleo-Ardèche river (ATrC in Figs. 5 and 6).

In seismic profile P7, the Ardèche tributary canyon appears as a shallow incision, north of the Ardèche canyon (Fig. 7). A canyon is visible at a shallow level west of the P4 seismic profile, reaching 250 ms in depth (Fig. 5). This canyon may be the upstream part of the tributary canyon of the Ardèche, with a N-S oriented course, east of the P7 seismic profile and the Lapalud area.

South of Mondragon, the 94MAR01 seismic profile, located alongside the Rhône, shows the occurrence of a deep incision by the Messinian Cèze river (CC in Fig. 10), reaching a depth of ∼530 ms. Immediately to the south of the Cèze canyon, another incision reaching a depth of ∼390 ms occurs on the seismic profile, most probably corresponding to the Tave canyon (TC in Fig. 10) that today joins the river Cèze just west of the seismic profile.

4.2.3 Pliocene units (Pl-r, Pl1, Pl2, Pl3 and Pl4)

Up to five successive seismic units of Pliocene age fill the canyon. They mostly onlap incisions, paleo-topographies and anticlines (Fig 7, 8). Some of the Pliocene infill in the central part of canyons displays a syncline-like geometry, that we interpret as the result of differential compaction.

On the eastern side of seismic profile P4 (Fig. 5), the basal unit is made of chaotic reflectors and presents an irregular top, pinched towards the western border of the canyon. This unit is interpreted as a Mass-Transport Deposit (MTD), and is referred to as Pl-r1 (r stands for reworked sediments and the 1 means unit Pl1 onlaps onto it). This Pl-r1 unit is not a stratigraphic unit but rather a seismic facies unit. On profile P4, it covers the western side of the Messinian canyon down to the centre, and also appears on seismic profiles P4, P2, P9 and P10 (Fig. 5, 6 and 10) reaching up to 150 ms twt (∼165 m with the overall Pliocene filling velocity of 2200 m/s) in seismic profile P2, south of Pierrelatte. The lateral extension of this unit appears to be controlled by the paleo-topography of the Messinian incisions (Fig. 10). Interestingly, a coeval basal Pl-r1 unit also appears on the flank of the Messinian canyon in both the Cèze canyon (Fig. 10) and the Rhône canyon where it cuts through the Mesozoic carbonates. In the Tricastin area, this unit is mostly onlapped by unit Pl1 and to a lesser extent by unit P2 (Fig. 10). Two other reworked series are visible higher up in the Pliocene filling and are described hereafter with respect to the seismic unit onlapping onto it.

The Pl1 unit onlaps the Pl-r1 unit and, overall, fills the base of the incisions. It mainly presents two seismic facies: seismic facies 1 (SF1 in Fig. 9), made of relatively continuous and moderate to high amplitude parallel reflectors, and seismic facies 2, made of moderate to low amplitude discontinuous reflectors (SF2 in Fig. 9). It should be noted that seismic facies of Pl1 unit differs on seismic profiles shot in the Donzère-Mondragon channel (see P8D16 seismic profile in Fig. 8) from those on the other seismic profiles, probably due to the difference in acquisition techniques. The Pl1 unit was drilled in Saint-Paul-Trois-Châteaux and Pierrelatte wells and is characterised by sands and clay/marls (Fig. 4). In Saint-Paul-Trois-Châteaux, this interval, which is interpreted as Pl1 (Fig. 8) corresponds to the 102 m azoic interval made of sands and clays which were not stratigraphically assigned by the drilling. Profile P8D16 clearly shows that the upper high amplitude reflectors in the K2 seismic unit are eroded by the MES. The Pl1 unit is also the thickest unit filling the canyon, where its thickness reaches 300 m in the Rhône valley (Fig. 5) and approximately 330 m in the Ardèche canyon (Fig. 6, profile P5). Higher amplitude reflectors might reflect sand bodies since in Saint-Paul-Trois-Châteaux, the Gamma-Ray (GR) log exhibits low values typical of sandy intervals, like in Codolet well in the Cèze canyon (Fig. 9).

The Pl2 unit which conformably overlies the Pl1 unit, is relatively thin compared to Pl1, with a thickness ranging from ∼80 to ∼120 m, considering a mean velocity of 2200 m/s. The Pl2 unit is characterised by two seismic facies: seismic facies 3, made of not wholly continuous parallel to divergent reflectors, moderate in amplitude (SF3 in Fig. 9) and seismic facies 4, marked by moderate amplitude discontinuous reflectors topped by higher amplitudes doublets (SF4 in Fig. 9). Pierrelatte and Saint-Paul boreholes show that unit Pl2 is more shaly and less sandy than unit Pl1. In the Cèze canyon, unit Pl2 is marked by many lignite layers (Fig. 4). The mean thickness of the Pl2 unit is relatively constant throughout the Tricastin area (calculated as ranging from ∼80 to ∼120 m), although with slight thickening southwards (profile P8D16, Fig. 8). The GR trend in Saint-Paul-Trois-Châteaux borehole shows a sharp transition to higher values marking the occurrence of a clay-richer interval (Fig. 9). In Codolet, the GR trend also shows higher values with spikes marking layers of blue clays and lignites (Fig. 9). A local MTD occurs on the northern flank of the Ardèche canyon (light-grey Pl-r2 unit on profile P5; Fig. 6). More interestingly, the top of unit Pl2 occurs at the steep shoulder of the Messinian canyon, where the canyon widens significantly (Figs. 6-8 and 10).

Pl3 conformably overlies Pl2, and onlaps onto the Messinian erosional surface. The contact between Pl2 and Pl3 is marked by a high amplitude couplet of reflectors, which are more visible on profiles P2, P4, P5 and P10. Thicker than unit Pl2, unit Pl3 is marked by three seismic facies: SF2 and SF3, as described for units Pl1 and Pl2, plus seismic facies 5, which is made of moderate to high amplitude divergent reflectors (SF5 in Fig. 9). Unit Pl3 was crossed by many wells (Fig. 4) and reveals a composition of grey clays and marls with local sand, silt or gravel layers. GR logs in both Codolet and Saint-Paul-Trois-Châteaux wells also show higher values, in agreement with the shaly content of the Pl3 unit (Fig. 10). On the western side of seismic profile P4 (Fig. 5), a local MTD (Pl-r3) is onlapped by unit Pl3. In the Cèze canyon, this layer is marked by the occurrence of thin layers of lignites. Even though unit Pl3 mostly onlaps older structures (Figs. 6 and 7), in seismic profile P8D16 (Fig. 8), unit Pl3 appears to be tilted onto the hanging wall of the Mondragon frontal blind thrust (see fault F2 in Fig. 10 as well). A detailed view of this sector is provided in Figure 11 and will be discussed later.

The shallow Pl4 unit is poorly imaged on the seismic reflection because processing did not focus on this interval. One could however observe that this unit appears to be conformably deposited above unit Pl3 (Fig. 5). In terms of internal structures, unit Pl4 is marked by low amplitude and discontinuous, sometimes transparent, reflectors, whereas in seismic profiles shot in the Donzère-Mondragon channel, unit Pl4 is characterised by high amplitude reflectors (Fig. 8). When drilled, unit Pl4 is marked by the well-known bluish clays and marls of the Pliocene, often found outcropping on the field (Ballesio, 1972). GR trends in Saint-Paul-Trois-Châteaux and Codolet both show high values consistent with the clay-rich filling (Fig. 9). It should be noted that the top of the seismic unit P4 is very difficult to interpret and that contact with the Quaternary is not clear.

4.3 Construction of the MES surface model

We used the Topo-to-Raster tool implemented in spatial analyst extension of Arcmap Desktop 10.6 version (ArcMap, ESRI) to construct a 3D morphology of the Messinian Erosional Surface. This interpolation method makes it possible to create DEMs with consistent hydrological networks (e.g. Hutchsinson et al., 2011). The Topo-to-Raster tool interpolates elevations, and considers constraints including isolines, measured elevations, and streams, for example. The output is a topography that is consistent with drainage resulting from river erosion, which is the case for the Messinian Erosional Surface.

In order to add constraints at depth, we carefully investigated the BRGM borehole subsurface database (https://infoterre.brgm.fr) and selected all boreholes that either crossed the base of the Pliocene or those that did not identify any Pliocene deposits but instead older strata below the Quaternary. The altimetric base of the Pliocene was then extracted to create a table containing X, Y and Z (Tab. S1 in supplementary material), which was then used to interpolate the MES in depth using the Topo-to-Raster tool. To this end, a combination of (1) the depth-converted SEM seismic horizon considering a constant VP velocity of 2200 m/s as explained above, (2) the depth of the Pliocene recorded in the boreholes located within the study area, (3) the altitude of the Pliocene outcrop boundaries taken from geological maps (see dashed line in Fig. 2), (4) the path of the Messinian talwegs deduced from our interpretation of reflection seismic data and hypothesised for the area for which no data are available (area between Malataverne and the Donzère-Mondragon Channel and the part of the Ardèche located north of Saint-Paulet-de-Caisson). (5) Some specific altimetric points were added downstream (the Rhône canyon at a depth of around -800 m south of Mondragon, Hollender pers. comm.; the Pujaut area where the canyon talweg is known to be at ∼ 900 m; Schlupp et al., 2001, Hollender et al., 2015b) in order to constrain the river flow in the GIS tool. Figure 12 shows the resulting digital model together with the input data. The model is a smoothed topography with a 50 m grid. The model is the result of combining the current topography retained outside the limits of the Pliocene filling of the canyon and the topography constructed from the various sources of information derived from our analysis as described previously. This topography can be compared with the most recent proposed by Roure et al. (2009) (see Fig. S6 in Supplementary material for more information). From south to north, the following topographic structures are observed: on the right bank, the Cèze canyon (imaged by ANDRA in the 1990s; Roux and Brulhet, 1997), on the left bank, the Aygues canyon (Bailly et al., 2015).

In the Tricastin area, the model differs from that of Roure et al. (2009) in that the outlet of the paleo-Ardèche is significantly offset to the south and connects to the paleo Rhône canyon south of Lapalud island. This path is controlled by seismic profiles P7 and P5. The confluence zone is not imaged, but its location to the south of the study area is partly constrained by the narrowness of Rhône canyon in the corridor between the Cretaceous shoulders of the Mondragon cluse and seismic profile P8C16 that shows no further incision toward the south.

East of Lapalud, the position of the Rhône canyon is constrained by outcrops of Pliocene marine sands and marls in the city of Bollène. Further information was available thanks to the presence of a few boreholes showing the hectometric thickness of the Pliocene. The canyon widens throughout the area between Lapalud and Pierrelatte, and has an asymmetrical transversal profile (gentle western edge and steep eastern edge). This widening seems to be associated with erosion facilitated in this area by subcropping soft Aptian to Cenomanian formations. Further north, between Pierrelatte and La Garde Adhémar, the canyon again becomes entrenched in the Cretaceous limestone, which forms a cluse at this point. The Z-shaped talweg is constrained by its repeated crossing of line P8D16. Further north, due to the absence of constraints, we chose to draw a straight line towards Malataverne Pass. Finally, the Paleo-Rhone diverges from the present-day river that flows in the “Donzère Défilé”; the canyon cut through the Malataverne Pass, as shown by detailed geological outcrop and subcrop mapping (Camus, 2003). Further north, currently no data are available to define the exact path, except those from scattered drill cores.

Fig. 11 Detailed view of the northern Mondragon anticline (seismic profile P8C16) showing the dip of Pliocene Pl3 seismic unit discussed in the text. Apparent depth (left side of the figure) and apparent dip (right side of the figure) are estimated given seismic velocities retrieved from depth conversion of the P4 seismic profile (Fig. 6).

Fig. 12 Input data (left) and resulted MES map (right) constructed by interpolating multiple parameters (current topography, depth of Pliocene basement in boreholes, depth from seismic profiles, watercourses imposed by visual analysis of seismic profiles). The topography presented is relative to sea level (hypsometry).

5 Discussion

5.1 Course, depth and morphology of the Messinian canyons

Seismic line P3 is the only one that crosses the Messinian Rhône canyon close to perpendicular. Line P3 shows an asymmetrical shape including: 1) a gently dipping western flank (∼10° or ∼18%, based on the velocities used in the depth conversion of profile P4) that mostly follows Cretaceous structural surfaces (even though slightly erosive), and 2) an eastern flank with an abrupt slope (∼16° or ∼29%). Thanks to the depth converted seismic profile P4, we were able to estimate the slope of the canyon to be ∼14° (∼25%) from the shoulder of the canyon, down to its bottom (Fig. 5). This asymmetry suggests that during Messinian period, the Rhône river progressively eroded the clastic Cretaceous series until it reached the more resistant Urgonian, controlled by the east-dipping structural slope.

The talweg of the Messinian canyon was not directly imaged north of Line P1. However, a past drilling campaign revealed the presence of the Messinian canyon below the Malataverne area (Fig. 12; Camus, 2003). The ultimate depth of the Messinian talweg remains unknown, but interpolation using data from the boreholes suggest a talweg at a depth of more than -220 m b.s.l. From that point southwards, seismic reflection depth-conversion suggests a talweg at a depth of -460-470 m b.s.l. at P2, -580/590 m b.s.l. at P3 and at least ∼700 m at P4. The Messinian canyon is not perpendicularly crossed by seismic profiles further south in the Tricastin area. Along the P5 seismic profile, the Ardèche canyon reaches ∼600 m b.s.l., meaning the Rhône canyon was even deeper close to Mondragon. Depth estimates made it possible to compute an average slope of the Messinian talweg of 2.8% (or 1.6°) between P2 and P5.

The Tricastin area is known to be the locus of an important increase in slope along the Rhône canyon (Clauzon 1982; Beaudoin et al., 1997) with an upslope profile estimated at between 0.1 and 0.2% (0.06 to 0.1°) compared to the downslope profile, which is estimated to be close to 1% (0.6°). The assumption of a knickpoint in the talweg resulted from an incomplete set of depth control data: the most important available data came from the Pierrelatte well crossing the base of the Pliocene series at -235 m b.s.l. However, recent seismic reflection profiles show that the Pierrelatte well did not reach the base of the canyon, which is much deeper, in the range of -460-470 m b.s.l. along P2. Our results imply a deeper MES in the Tricastin area, and consequently, a less steep slope of the talweg downstream of the study area. Considering a -700 m b.s.l. depth of the canyon at P4, -460-470 m b.s.l. at P2, we estimated a slope of between 0.6% and 0.9% in the Tricastin area. Our modelling of the MES gives an ∼1% (0.6°) slope from Malataverne to the confluence of Ardèche and Rhône rivers to the south. The results of this present study reveal a more regular Messinian profile of the canyon, at least up to Malataverne in the north. Additional subsurface data north of Malataverne would help identify the suspected knickpoint north of the study area.

Interestingly, on P1 and P2 (Fig. 6), the MES displays a flat surface parallel to both the top of K2 and the top of Pl2 seismic units, which are in lateral continuity across the shoulder of the Rhône and Ardèche canyons. This flat surface seems to either follow the structural surface of K2 seismic unit in seismic profile P1, or to erode through the K2 seismic unit in profile P2 (Fig. 6). To the south of the Ardèche canyon, in seismic profile P7 (south of AC in Fig. 7), the topmost surface of Pl2 seismic unit laterally corresponds to a flat surface following a structural surface on top of (or within) the K3 seismic unit. The shoulders of the Rhône and Ardèche canyons deepen toward the south: the shoulder of the Rhône canyon being approximately -30 m b.s.l. on the P1 seismic profile (given the velocities extracted from the depth converted seismic profile), -60 b.s.l. along the P2 seismic profile and the shoulder of the Ardèche canyon at approximately -100 m b.s.l. (Fig. 7). The flat surface might have been generated during a base-level stagnation, during either the fall or the rise in sea level.

In our study area, the Pl2 seismic unit, the top of which is laterally continuous with the flat surface, is made of sands and clays of continental origin (lignites in the Codolet well, Fig. 4). There is thus no direct evidence of marine deposits at that depth that would favour a wave-cut surface during the Pliocene reflooding like that seen offshore (Bache et al., 2012). It is thus very likely that this flat surface is the result of the development of an alluvial plain or a lake during a period of sea-level stagnation.

Considering the abandonment surface of the Miocene located at +360 m a.s.l. in Saint-Restitut, the difference in elevation between the flat surface and the previous abandonment surface of the Miocene is at least 400 m considering the shoulder of the Rhône canyon along the P1 seismic profile. This difference in altitude may increase to ∼460 m, considering the Ardèche canyon shoulder along P7 seismic profile. Could this flat surface have been generated during a sea-level drop of that magnitude? The drop in sea level during the first stage of the MSC is estimated to range from 400 m (Beaudoin et al., 1997) to ∼500 m (Gargani and Rigollet, 2007) in the Rhône canyon. The depth of the flat surface fits within this range. This kind of flat surface has already been observed further south in the Durance Messinian canyon (Clauzon et al., 1995) and its genesis has also been linked to the first stage of the MSC sea-level fall (Hippolyte et al., 2020).

Is it possible that this surface was generated during the sea-level rise following the deposition of evaporites in the central Mediterranean basin? The lateral continuity between the flat surface and the Pl2 unit argues in favour of a surface-forming process that would genetically link those features during the Pliocene sea-level rise. Periods of stagnation in sea-level rise during the Pliocene reflooding produced wave-cut surfaces in the Gulf of Lion and offshore Provence (Bache et al., 2012, Tassy et al., 2014), but none of those surfaces reach the altitude of the flat surface identified underneath Tricastin. The Ardèche area was subject to the development of a subsurface karstic network where the sea-level drop of the MSC was responsible for a deep flooded karst and a rise in sea level produced per ascensum terracing of horizontal networks at +90 and +130 m a.s.l. (Mocochain et al., 2006a; Arfib and Mocochain, 2022). In addition, several wave-cut surfaces have been identified 20 km south of Tricastin at +60 m and +100/105 m (Caziot, 1890). These observations indicate oscillating reflooding during the Pliocene, however the limited lateral extension of the above-mentioned sea-level stagnations suggest that none of those could have given rise to such a large flat surface as that observed on the seismic profiles in the Tricastin area. The development of such a wide flat surface at the top of the Pl2 unit almost certainly marks a base-level stagnation that lasted for a substantial period of time during the Pliocene reflooding. A forward stratigraphic modelling approach including sea-level fluctuations and stagnations during the reflooding would provide more insights into the creation of this surface.

5.2 Interpretation of the canyon infill

In the Western Mediterranean, most coastal rias Pliocene in age are filled with regressive deposits marked by a Gilbert-type fan delta superimposed over debris flows (Clauzon et al., 1995; Breda et al., 2007; Duvail, 2008). Such deposits are mostly found in areas located close to the apex of the sedimentary system where transgressive deposits due to the reflooding of the Mediterranean Sea are missing. In the Tricastin area, the successive Pl1 to Pl4 units, evolving from continental sands to marine clays and shales, mark an overall transgression during the filling of the canyons, but no clear deltaic systems such as Gilbert-type fan delta were observed in the seismic profiles. The stacking pattern of the study area starting with sands evolving to marine clays is a pattern that has already been identified in both the Pujaut graben (Ballesio, 1972) and the Cèze canyon (Ferry et al., 1997). Ferry et al. (1997) also concluded that within the Rhône canyon, the Pliocene denotes the complex transgressive system tract of the Pliocene sea level rise, which apparently did not develop in the coastal rias due to low sedimentation rate.

The significance of the various MTDs identified on the seismic profiles needs to be discussed. The basal Pl-r1 unit, observed on seismic profiles P2, P4, P9 and P10, is located on the western flank of the Rhône canyon and seems to pinch out eastwards on profile P2 (Figs. 6 and 10). More importantly, the basal Pl-r1 unit pre-dates deposition of unit Pl1 which onlaps over the MTD. On seismic profile P9 and 94MAR01 (Fig. 10), the base of unit Pl1 appears to be lower than the base of unit Pl-r1, suggesting an incision occurred between the deposition of unit Pl-r1 and Pl1. No wells have ever been drilled in this unit, however some debris flows, sometimes called block formation (“Formations à gros blocs” in the Languedoc area; Ambert, 1989), polygenic conglomerates or even olistoliths and breccias, have been found unconformably over the Miocene and below the Pliocene shales, covering the flank of the Cèze canyon (Ballesio and Truc, 1978), along the Rhône valley (Ballesio, 1972; Clauzon, 1978; Schlupp et al., 2001) and even offshore (Bache et al., 2012). Their origin along the Rhône and the Cèze canyons is usually interpreted as the result of hillslope processes (scree aprons). Ballesio and Truc (1978) also pointed to the close relationship between the occurrence of these thick bodies and the vicinity of fractures and faults, which are assumed to have been active during Alpine shortening. All these criteria lead us to assume Pl-r1 deposits were emplaced during the second stage of the MSC (5.6–5.46 Ma), i.e. during the massive drawdown in sea level. An emblematic example of such deposits is the Carros Breccia, which is located on the western flank of the Var Messinian canyon (Clauzon, 1978) between the MES and the marine Pliocene foresets. Map and facies analyses clearly show that they are scree aprons deposited during the subaerial phase of erosion at the base of an active relief in this case, a thrust front.

On seismic profile P5, the tributary of the Ardèche river (ATrC in Fig. 6) is partly filled by a MTD, sandwiched between unit Pl1 and the Pl2 unit that onlaps onto it (light-grey polygon in Fig. 6). This MTD dips toward the canyon axis, identifying its origin as the steep Turonian-Coniacian shoulder to the north. The fact that this MTD is intercalated within the Pliocene filling suggests that slope instabilities remained active at various places and times during the filling of the canyon caused by reflooding. This continuing activity is also supported by the intercalation of another MTD between the units Pl2 and Pl3 on seismic profile P4 (Fig. 5), most probably originating from the same Turonian-Coniacian series of the Lapalud sector. The generation of MTDs triggered during sea-level rise has not been recognised in seismic profiles either within the Cèze canyon or within the Pujaut graben further south (Ferry et al., 1997; Schlupp et al., 2001). However, the breccias and olistoliths pointed out by Ballesio (1972) at the surface along the Cèze canyon may be the best analogues of those higher MTDs. Similar boulders are recorded in Roussillon at the boundary of the granitic basement along the Têt river canyon (Clauzon et al., 2015a): in this case, they are imbricated with marine foresets but against the erosional surface. As a result, it is hard to decipher if they are coeval to the sub-aerial period or if the slope continued to be eroded after the submersion.

5.3 Identification of structures and deformations

Identification of faults and of their apparent vertical displacement depends on the seismic resolution (dominant frequency spectrum of our dataset between 50 and 100 Hz) as well as on processing issues (see Sect. 2). It is therefore difficult, if not impossible, to detect faults with a vertical offset of less than 50 m in the Urgonian and 4 to 10 m in the most superficial levels. However, experience shows that ruptures usually occur on pre-existing structures (as was the case on the La Rouvière fault during the Le Teil earthquake, see Ritz et al., 2020). Consequently, it is important to identify deep seated faults to identify potential fault rupture for seismic hazard assessment.

In addition, the identification of “structures” on seismic profiles should be interpreted with caution as most seismic profiles have not benefitted from static correction, particularly in our case where there is a significant “pull down" effect of reflected waves due to the presence of the Plio-Quaternary filling of the canyon, which may produce fault-like or fold-like structures (see example in Hanot and Thiry, 1999; Beccaletto et al., 2011). However, the structures described below display offsets that are large enough to be considered as faults or folds.

5.3.1 Faults and folds affecting the Cretaceous series

N-S oriented profiles P8D16, P6 and P7 cross a broad ENE-WSW oriented anticline affecting the Cretaceous series (Urgonian and locally Aptian). The maximum amplitude of this anticline is around 500 m computed with 2200 m/s for post-Urgonian sequences, for a distance of around 7 km. This structure correlates with the Echavarelles anticline (Fig. 13) and also accounts for the small outcrop of Urgonian limestone in Pierrelatte due to additional normal faulting of the anticline axial plane. This large structure is affected by shorter wavelength series of folds, ENE-WSW oriented (Fig. 13), which can be identified on profiles P6 and P7 (Fig. 7). These N80-trending structures run parallel to the Mondragon fold in the Uchaux massif, which was created during the Pyrenean shortening (between the Palaeocene and Oligocene; see Ballas et al., 2014 and references therein). At La Garde-Adhémar, the Echavarelles anticline is unconformably covered by Oligocene marls that overlie Aptian (Gargasian) marls to the south and Albian glauconitic sands to the north.

The Echavarelles anticline hinge is affected by a large fault (hereafter PF for Pierrelatte Fault, Fig. 13) which offsets the Urgonian top by around 0.2s twt (500 m considering a VP velocity of 4,500 m/s in the Urgonian) on profile P7. This offset reaches at least 0.1 s on profile P6 and is clearly visible on profile P8D16 even though its offset remains difficult to compel (Figs. 7 and 8). We also identified the Pierrelatte fault at the western end of the P2 profile (Fig. 6), but again the offset of this fault cannot be estimated due to the fault orientation running oblique to the seismic profile. The Pierrelatte fault shows a normal component, it offsets the Urgonian limestones and the overlying Aptian and is truncated by the MES. The age of this fault is post Aptian and predates the formation of the MES. Ballas et al. (2014), described normal faults of nearly the same strike in the Uchaux massif, 12 km to the south (Mondragon, Mornas fault) and attributed their activity to the Oligo-Early Miocene extension. The fault could therefore have been active during this geodynamic episode. However, the position of this Pierrelatte fault in the extrados of this regional-scale anticline could also be explained mechanically by the local extension of the anticline’s hinge during the Pyrenean tectonic episode. This fault has no clearly identified extension to the east (a few faults of similar direction are mapped to the north of Clansayes) and is not mapped on the eastern limb of the Saint Remèze dome to the south of Bourg Saint Andéol, where one would expect to find it. Apart from this major fault identified on 3 profiles, minor structures have been interpreted cutting the Cretaceous series, and do not appear affect the overlying Pliocene fill.

To the south of the Tricastin area, on profile P8C16, the deformation is clearly more pronounced. Combes and Carbon (1997) and Ballas et al. (2014) have shown that the Mondragon-Uchaux anticline is a fault-propagation fold Pyrenean in age, later intersected by ∼E-W normal faults. The Pyrenean shortening is marked in profile P8C16 by the occurrence of south-dipping blind thrusts noted F1 and F2 (Figs. 8 and 13), F2 being the northern front thrust accommodating the Mondragon fault propagation fold. Projection of the Mondragon borehole shows fault F2 is responsible of the duplication identified in the Mondragon borehole (two occurrences of the Aptian/Urgonian succession: BSS002CLTG). The F1 fault located further south, is most probably a splay of the Mondragon thrust, as proposed by Combes and Carbon (1997). P8C16 displays a 0.2s twt offset of K3 series, some 2 km north of F2. Although the seismic image does not clearly show the nature of this offset, the Oligocene Dessoulière Fault that strikes N80°E (Combes and Carbon, 1997) may correlate with this structure, noted F3 (Fig. 13).

Fig. 13 Structural map showing the main structures identified on seismic profiles. Red structures are post-Cretaceous and apparently sealed by the MES. The Pierrelatte fault (PF), which has a normal component, is mechanically consistent with the regional stress field during Pyrenean orogeny if it acts as an extrados structure along the Echavarelles anticline hinge. The F3 normal fault would have behaved as a normal fault during Oligocene extension (Combes and Carbon, 1997; Ballas et al., 2014) . Red structures (Mondragon frontal thrust F2 and F1) might be the only ones active during or after the Pliocene, see text for more details.

5.3.2 Deformation of the Plio-Quaternary canyon filling

Seismic profiles reveal no folding or reflector offset at the scale of the seismic reflection resolution. However, most profiles display thickening of the Pliocene sequence toward the canyon axis. While a “pull-down" effect is likely (thicker series are slower and therefore appear even thicker on double-time sections), the amplitude of this syncline geometry is nonetheless significant. It can also be seen on profile P4, which has benefitted from static corrections (cf. Fig. 5, profile P4 corrected for depth). We suggest that this long-wavelength deformation may be linked to the differential compaction of sediments along the canyon edges as it has been demonstrated and modelled close to escarpment (Carminati and Santantonio, 2005). Such pattern has also been identified in the Cèze canyon (Ferry et al., 1997).

Finally, above the Mondragon fault propagation fold, on profile P8C16, we observed a south-dipping bedding in the Pliocene fill (unit Pl3) on the hanging wall (black arrows in Fig. 11, with apparent dip reaching ∼15°) and north dipping bending reflector over the emergence or the Mondragon thrust (white arrows in Fig. 11, with apparent dip reaching ∼30°). In our interpretation, these layers appear to be truncated by the upper level Pl4. This suggests a Pliocene reactivation of the Mondragon thrust, during Pl3 and prior to Pl4 deposition. This north-south component shortening would be coherent with compressive Quaternary activity observed in Courthézon near the Nîmes Fault and on the La Rouvière Fault on the Cévennes Fault system (Bellier et al., 2021; Ritz et al., 2021). These observations suggest that further investigation is needed to evaluate the recent (e.g. Plio-Quaternary) tectonics in this area.

If the observed clinoforms originate from a sedimentary process (black and white arrows in Fig. 11), they could represent steeply dipping clinoforms such as foreset beds of a Gilbert-type fan delta, widely recognised along the Mediterranean coasts during the reflooding of the Mediterranean Sea. However, if this is the case, such a thick delta would require an important lateral feeder system, most probably from a tributary of the Rhône river, which is enigmatic.

5.4 The Ardèche: from the canyon to the karst

The present-day Ardèche canyon is incised into the Urgonian facies limestone and displays an average slope of 1%. Marine Pliocene sediment found in the lower part of the present-day gorges indicates that the current river course is inherited from the Messinian evolution (Ballesio, 1972; Pascal et al., 1989; Mocochain et al., 2006b).

Downstream of the gorges, the seismic profiles interpreted in this contribution clearly shows the excavation of a deep and narrow canyon, visible on profile P7 at a depth of around -460 m b.s.l. (425 ms twt, Fig. 7) and on profile P5 at a depth of around -600 m b.s.l. (500 ms twt, Fig. 6). The canyon cuts through the post-Urgonian Cretaceous formations that are richer in clastics and easier to erode than the Urgonian bio-constructed limestone. This suggests the incision was halted by the Urgonian limestone. The canyon axis interpreted on the P7 and P5 profile, which are 2300 m apart, indicates a longitudinal slope of approximately 6% towards the SE, which is equal to the regional dip of the Urgonian limestone.

The Messinian Ardèche longitudinal profile therefore displays a significant increase in slope in the area of the present-day outlet of the gorges. Incision into the silico-clastic dominated K3 and K2 sequences reached the Urgonian limestone immediately downstream of this point (see P5 on Fig. 6). This resulted in a hydrogeological window that allowed the river water flux to take a shortcut through the deep karstic network incised into the Urgonian limestone, to connect with the base level represented by the Rhône canyon. As soon as these karstic networks developed beneath the gorges, the erosion within the canyon downstream may have slowed down considerably. This decrease in subaerial erosion would explain why the Ardèche canyon did not cut deeper into the Urgonian limestone massif.

The deep karstic system ceased its activity during reflooding of the Mediterranean. The deep karst was then buried under Pliocene sediments, leaving only the upstream parts of the submerged networks, which currently form a set of Vauclusian systems (Arfib and Mocochain, 2022; Mocochain et al., 2006a, 2011).

6 Conclusions

Our study showed that the Rhône canyon in the Tricastin area is up to several hundred metres deeper than previously estimated (Clauzon, 1982; Roure et al., 2009). Depth converted seismic profiles gave us access to seismic velocities of the canyon infill and allowed us to estimate the depth of the Rhône canyon to range from -400 to -700 m b.s.l. in the area. In addition, a deep and narrow canyon was imaged downstream of present-day Ardèche river, which we interpret as the Messinian Ardèche canyon. This canyon most probably merges with the Messinian Rhône canyon south of Mondragon. Prior to this study, this deep Ardèche canyon has never been imaged (see Roure et al., 2009 map in Supplementary material). The Rhône and Ardèche canyon morphologies widen at various locations following a flat surface located between -60 and -100 m b.s.l. We interpret this relatively flat surface to be the result of a sea-level stagnation during the Pliocene reflooding.

The infill of the Rhône canyon started with the occurrence of a large basal MTD, located below the Tricastin nuclear site, certainly deposited during the sea-level fall of the Mediterranean Sea after 5.6 Ma as it is onlapped by stratified Pliocene continental deposits. The overlying Pliocene infill of the Rhône and Ardèche canyons consists of four seismic units that represent an overall transgressive trend related to the reflooding of the Mediterranean Sea at the end of the MSC. At least two other generations of MTDs took place during the rise in the level of the Mediterranean Sea. Their proximity to interfluvial hills made of Cretaceous clastics certainly explains their genesis. No regressive system tract was found underneath Tricastin, meaning it has been eroded above the present-day Rhône alluvial plain.

From a structural point of view, the northern part of the Tricastin area is marked by ENE-WSW folds affecting the Cretaceous, parallel to structures in the Uchaux massif. The seismic reflection data indicate that the main structure is the continuation of the Echavarelles anticline which is faulted in its extrados (Pyrenean faulting?) in the same direction, marking the Pierrelatte fault. The Uchaux structures can be extrapolated downwards, where they show south-dipping thrusts marking the fault propagation fold of the Uchaux massif. Pliocene reflectors appear to be deformed above the hinge zone of the Mondragon anticline, reflecting possible recent reactivation of that structure. The latter is of importance for future research on recent activity of tectonic structures as the area is subject to low to moderate seismicity.

In a broader perspective, ongoing research dealing with seismic hazard and seismic wave amplification below nuclear sites will greatly benefit from this study as we have contributed a considerably more detailed subsurface geological model.

Supplementary Material

Table S1: Depth of the base of Pliocene in boreholes or outcrops. Depth values are in meter.

Figure S1: Location of seismic lines studied on a simplified geological map, in addition to boreholes used for analysis and construction of the MES model (see Table S1).

Figure S2: Stack in time after conventional processing.

Figure S3: Final interval velocity model.

Figure S4: RMS velocity model derivated from the time based interval velocity model that has been tested in a stack process to validate its reliability.

Figure S5: Stack with « deep » statics in time → improvement of the quality of the stack under the Paleo-Ardéche canyon.

Figure S6: Comparison of the MES surface at depth from Roure et al. (2009) and our new model of the MES at depth. We used the same color-scales for both maps. Note the different course of the ardéche and Rhône rivers.

Figure S7: Uninterpreted seismic lines used in this manuscript.

Acknowledgments

This work received support from the French government under the France 2030 investment plan, as part of the Initiative d’Excellence d’Aix-Marseille Université − A*MIDEX (AMX-19-IET-012) and from the Research Federation ECCOREV (FR 3098; Aix-Marseille Univ., CNRS, INRAE, IRSN, CEA, Univ. Toulon, Univ. Avignon, Univ. Nîmes). Interpretation of seismic profiles benefited from academic licence of S&P Kingdom software provided to ISTeP (Sorbonne Université). Ludovic Peignard is also warmly thanked for providing a research licence of the EarthQuick software. The authors thank Francois Lavoué for fruitful discussions on the canyon geometry and model implementation for the numerical simulation of site effects. Sincere thanks are also due to Hervé Jomard and Stéphane Baize for discussion and field participation on the Pliocene paleo-beaches and paleo-notches. Vincent Ollivier, Nicolas Cathelin, Jean-Claude Hippolyte, Jules Fleury, Estelle Hanouz, Jean-Pierre Suc and Gilles Dromart are thanked for fruitful discussions during a workshop organized in November 2021 and various fieldtrips. Fabrice Hollender from CEA is warmly thanked for providing reprocessed data on the seismic profile 94MAR01. Sylvain Pouliquen and Christophe Vergniaux from EDF for discussion and data exchanges on the seismic processing with CDP company, Didrik Vandeputte and the city hall of Saint-Paul-Trois-Châteaux who provide us the eponymous borehole documents. Seismic data use and publication was granted by ENGIE company to IRSN in the framework a non-exclusive licence contracted in June 2020. We thank Hubert Mignot who made these data available. We would also like to pay tribute to our colleague Christophe Clément, who helped us obtain the ENGIE profiles and supported us at the beginning of this project. Michel Séranne and Ferran Estrada are warmly thanked for their constructive reviews that helped improving significantly the manuscript.

References

All Figures

Fig. 1 (A) Synthetic structural map of the study area (black box) showing the main tectonic domains mentioned in the text and the course of the Rhône and Ardèche rivers in dotted blue lines. (B) Seismotectonic map of the lower Rhône Valley between the Cévennes Fault and Nîmes Fault. Plotted historical earthquakes are from SisFrance Database (www.sisfrance.net: data from SisFrance, BRGM, EDF, IRSN, 2022), magnitudes from FCAT (Mw catalog − Manchuel et al., 2018). BDFA refers to the Potentially Active Fault Database (Jomard et al., 2017), paleoseismological works reported in Baize et al. (2002) and Bellier et al. (2021).

In the text

	Fig. 2 Simplified geological map of the studied area (modified from Champenois et al., 1971; Pascal et al., 1989; Masse et al., 1980; Debelmas et al., 2004) showing location of seismic profiles used in this study and the main boreholes mentioned in the text.
In the text

	Fig. 3 Synthetic stratigraphic log of the study area, from the Hauterivian up to Quaternary. Seismic units interpreted on seismic profiles are reported on the right side of the log.
In the text

	Fig. 4 Cross-correlations of the main boreholes intersecting the base of the Messinian erosional surface in the Tricastin area containing lithological data.
In the text

Fig. 5 Seismic profile P4, intersecting the Rhône canyon and a tributary of the Ardèche canyon, in two-way travel time (twt) on top, and converted in depth at the bottom. The Grand Malijac borehole is projected on the depth converted seismic profile for reference. Vp velocities obtained from the depth-processing are drawn on the section for each formation. The MES is drawn in dotted red line. See Figures 4 and 5 for stratigraphic legend and text for more details. RC: Rhône canyon; ATrC: Ardèche tributary canyon.

In the text

	Fig. 6 Seismic profiles P1, P2, P3 and P5 interpreted in twt. Faults are drawn in solid black lines. The MES is drawn in dotted red line. See Figures 4 and 5 for stratigraphic legend and text for more details. PF: Pierrelatte Fault; RC: Rhône canyon; AC: Ardèche canyon; ATrC: Ardèche tributary canyon.
In the text

	Fig. 7 Seismic profiles P6 and P7 interpreted in twt. Faults are drawn in solid black lines. The MES is drawn in dotted red line. See Figures 4 and 5 for stratigraphic legend and text for more details. PF: Pierrelatte Fault. AC: Ardèche canyon; ATrC: Ardèche tributary canyon.
In the text

	Fig. 8 Seismic profiles P8D16 and P8C16 interpreted in twt. Pierrelatte and Mondragon wells are projected on the seismic profiles. Faults are drawn in solid black lines. The MES is drawn in dotted red line. See Figures 4 and 5 for stratigraphic legend and text for more details. PF: Pierrelatte Fault; RC: Rhône canyon.
In the text

	Fig. 9 Comparison of lithologies, Gamma-Ray (GR) signals and seismic facies of the Pliocene filling units Pl1 to Pl4 within the Saint-Paul-Trois-Châteaux and Codolet wells. GR curves are redrawn from printed electric logs. Seismic facies SF1 to SF5 are detailed in the text. Profile P8D16 is shown in Figure 8 and profile 94MAR01 in Figure 9.
In the text

	Fig. 10 Seismic profiles P9, P10 and 94MAR01 interpreted in twt. The MES is drawn in dotted red line. See Figures 4 and 5 for stratigraphic legend and text for more details. RC: Rhône canyon; CC: Cèze canyon; TC: Tave canyon.
In the text

	Fig. 11 Detailed view of the northern Mondragon anticline (seismic profile P8C16) showing the dip of Pliocene Pl3 seismic unit discussed in the text. Apparent depth (left side of the figure) and apparent dip (right side of the figure) are estimated given seismic velocities retrieved from depth conversion of the P4 seismic profile (Fig. 6).
In the text

	Fig. 12 Input data (left) and resulted MES map (right) constructed by interpolating multiple parameters (current topography, depth of Pliocene basement in boreholes, depth from seismic profiles, watercourses imposed by visual analysis of seismic profiles). The topography presented is relative to sea level (hypsometry).
In the text

Fig. 13 Structural map showing the main structures identified on seismic profiles. Red structures are post-Cretaceous and apparently sealed by the MES. The Pierrelatte fault (PF), which has a normal component, is mechanically consistent with the regional stress field during Pyrenean orogeny if it acts as an extrados structure along the Echavarelles anticline hinge. The F3 normal fault would have behaved as a normal fault during Oligocene extension (Combes and Carbon, 1997; Ballas et al., 2014) . Red structures (Mondragon frontal thrust F2 and F1) might be the only ones active during or after the Pliocene, see text for more details.

In the text