Matthieu Caparos

Valerie Louis

Fabrice Wendling

Jean-Pierre Pierre Vignal

Didier Wolf

Extraction of reproductible epileptic patterns on scalp EEG

Keywords: EEG, epilepsy, seizure, similarity, edit distance

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Epilepsy is characterized by the sudden abnormal synchronization between a large group of neurons [START_REF] Witte | Special issue on epileptic seizure prediction[END_REF]. Temporal lobe epilepsy (TLE), the most common form of partial epilepsy in adults, is refractory to medical treatment in 20% to 30% of the cases. Surgery, which consists in removing the epileptogenic zone (EZ), is then the last resort for about half of patients suffering from it [START_REF] Iasemidis | Epileptic seizure prediction and control[END_REF]. Localization of epileptic seizure onset is thus a really important point before considering a surgical operation [START_REF] Engel | Research on the human brain in an epilepsy surgery setting[END_REF][START_REF] Ebersole | Localization of temporal lobe foci by ictal eeg patterns[END_REF]. Actually, the base of the surgery is that EZ and the epileptic network is stable for all seizures of a patient. Many techniques are used to produce a reliable diagnosis, from clinical signs to imaging (MRI) [START_REF] Asano | Ictal localization of temporal lobe seizures with scalp/sphenoidal recordings[END_REF][START_REF] Knake | Ictal functional tcd for the lateralization of the seizure onset zone -a report of two cases[END_REF][START_REF] Lee | The lateralizing and surgical prognostic value of a single 2-hour eeg in mesial tle[END_REF]. Among these techniques, electrophysiological methods are the most adapted to better understand the organization and propagation of the brain electrical activity of a seizure (Wendling et al., 2001b). Electrophysiological methods are di-vided into two groups: scalp ElectroEncephaloGraphy (EEG) and invasive techniques StereoElec-troEncephaloGraphy (SEEG), ElectroCortiGraphy (ECoG).

Signal processing methods have been successfully used on SEEG to complete visual inspection of seizure, in several studies (Wendling et al., 2001b;[START_REF] Van Der Heyden | Non linear analysis of intracranial human eeg in temporal lobe epilepsy[END_REF][START_REF] Wu | Segmentation and classification of eeg during epileptic seizures[END_REF]. It was showed that seizure evolution was closely linked to EZ organization and allowed to objectively characterize the seizure. The implantation of electrodes directly in cerebral structures let the electrophysiological signal be recorded directly at his source and then its evolution inside temporal lobe is relatively easy to described. Few similar researches have been proceed on scalp EEG. This may be explained by at least two factors. First, scalp EEG is a global reflection of cerebral activity, so during seizure, EEG signal correspond to a projection of the brain activity on scalp electrodes. Second, surface EEG signal is often artefacted, especially during seizures while uncontrolled muscular and eyes movements may occur. The main difficulty for the signal processing is then the extraction of a non-stationary event from a non-stationary noise. In this context, classical pattern extraction (mean detectors, spike detection method, ...) that have been applied to scalp EEG, prove to be non sufficiently robust. Another approach could then be the use of the decorrelation property between epileptic event and noise. As scalp EEG data is multidimensional, relationship between EEG channel could reinforce epileptic event hypothesis.

Measure of relationships between scalp electrodes is a technique often used in neurophysiological applications to point out the way cerebral structures work together. Concerning partial TLE, several methods of relationship measure have been successfully applied to SEEG signal [START_REF] Lehnertz | Non-linear time series analysis of intracranial eeg recordings in patients with epilepsy -an overview[END_REF][START_REF] Mormann | Epileptic seizures are preceded by a decrease in synchronization[END_REF]Wendling et al., 2001b;[START_REF] Chavez | Spatiotemporal dynamics prior to neocortical seizurres : Amplitude versus phase couplings[END_REF] in order to better understand interdependencies between anatomical structures during TLE seizures. The depth synchronization may be expressed by correlations on scalp electrodes, even if the transfer function of electrophysiological activity from depth to surface is only partially defined in literature. This assumption is reinforced by [START_REF] Franaszczuk | Application ot the directed transfer function method to mesial and lateral onset temporal lobe seizures[END_REF], who shows that the cause of the rhythmic activity, observed in ECoG, is an increase of functional couplings between a certain amount of brain regions.

The aim of this study is to characterize the evolution of TLE using scalp EEGs, in order to enhance preoperative diagnosis. In literature, several parametric and non-parametric methods are described to estimate these synchronizations at seizure onset: Directed Transfer Function in ECoG by [START_REF] Franaszczuk | Application ot the directed transfer function method to mesial and lateral onset temporal lobe seizures[END_REF], the non-linear regression method introduced by Pijn (1990) and applied to SEEG by (Wendling et al., 2001b), phase synchronization developed in [START_REF] Chavez | Spatiotemporal dynamics prior to neocortical seizurres : Amplitude versus phase couplings[END_REF],.... In this study, non-linear non-parametrical methods are preferred because they take into account nonlinearity in observed electrophysiological system [START_REF] Ansari-Asl | Quantitative comparison of signal analysis methods for characterization of brain functionnal connectivity[END_REF]. In [START_REF] Wu | Segmentation and classification of eeg during epileptic seizures[END_REF], authors developed algorithms to group seizure of similar morphology using scalp EEG. They based their method on computation of the edit distance, to give a quantitative information on seizure similarity. After adequate signal coding of non-linear regression signal, the method described in [START_REF] Wendling | Extraction of spatio-temporal signatures from depth eeg seizure signals based on objective matching in warped vectorial observations[END_REF] is employed to compute the similarity.

It the second section of the paper the procedure used for data base collection is summarized. Signal processing methods used to point out non-linear correlation between electrodes is explained. Signal encoding that leads to the extraction of similarities between two seizures is shown in following subsection. In third section results intra and inter patient is shown. Finally, these results are discussed and the conclusion on this work is presented.

MATERIAL AND METHODS

Patients

Seizures of forty three patients presenting TLE are studied. Two seizure for each patients were recorded in order to be able to point out intrapatient reproducible patterns. They are aged 16-45 years old and all underwent long term preoperative EEG recording at the neurology unit of Nancy's CHU (France). A comprehensive evaluation including detailed history and neurological examination, neuropsychological testing, magnetic resonance testing (MRI) study and interictal and ictal SPECT has been performed for all of them. The diagnosis of a partial TLE has been confirmed and they all underwent successful surgical resection of the EZ.

Data selection and signal acquisition

EEG signal was recorded on a Micromed system, according to 10/20 system using a 24 channels referential EEG montage. Data are sampled at 512Hz and recorded to hard disk. For further signal processings, eight couples of electrodes of the longitudinal montage are selected with the help of the neurologist. In order not to loose any usefully epilectical information that may be contained in EEG signal no pre-processings are applied, except a digital filtering in the [0.5 -30]Hz frequency band. As frequency of interest are under 30Hz, data are digitally filtered by a [0.5 30]Hz band-pass FIR filter and downsampled to 256Hz, to improve the signal processing time. Data files are then divided into two groups:

• 27 patients (55 seizures) were used to develop our methodology and optimize parameters

Seizure comparison by similarities extraction

The different steps of the signal processing methodology that leads to the extraction of reproducible patterns between TLE seizures will be detailed in the following. The methodology may be divided into three main steps: a) First, the nonlinear correlations between some EEG channels are computed to characterize the synchronisation of cerebral regions at the seizure onset. b) Then, the evolution of the correlation is encoded. c) Finally, the similarities between encoded seizures are extracted by a dynamic programming algorithm.

Characterization of correlations Before a seizure, structures of the brain organize themselves in synchronized neural network (Wendling et al., 2001b;[START_REF] Van Der Heyden | Non linear analysis of intracranial human eeg in temporal lobe epilepsy[END_REF]. The seizure activity can then be seen as a flow of energy from one part of the brain to another [START_REF] Franaszczuk | Application ot the directed transfer function method to mesial and lateral onset temporal lobe seizures[END_REF]. To measure the relationship that is established in such cases, two main types of techniques have been used:

• Techniques based on parametric methods (for example AR-modelization of the system), • Techniques based on non parametric methods in the correlation space.

Both kind of methods are joining themselves in the aim, which is identifying interdependency of groups of neurons that are working together, through the exchange of energy flow. ARmodeling techniques present the main disadvantage to first search a model that suits the signal, implicating all problems such as the choice of model order, error in the parameters identification, definition of a time of quasi-stationary duration, etc... . All these not-well-defined parameters tend to make validation of the method more difficult. In this paper non parametric method are then preferred.

In (Wendling et al., 2001b), the non-linear correlation method introduced in 1989 by (Lopes da [START_REF] Lopes Da Silva | Interdependence of eeg signals: Linear vs. nonlinear associations and the significance of time delays and phase shifts[END_REF], was successfully used to identify coupling between cerebral region on SEEG. The basic principle of this method is the identification of a non-linear dependency between two signals X and Y.

The non-linear correlation coefficient h 2 XY is obtained by considering the amplitude of signal Y (t + τ ) as a perturbed function of the amplitude of the signal X(t). The variance of this perturbed function is assumed to be the conditional variance of Y (t + τ ) given X(t) normalized by the variance of Y (t + τ ). This quantity is subtracted from 1 to obtain the quantity so called h 2 XY :

h 2 XY = 1 - var[Y (t + τ )/X] var[Y (t + τ )] (1)
where

var[Y (t+τ )/X] = argmin(E[Y (t+τ )-g(X(t))] 2 )
g is the approximation of the mapping function from x to y.

The computation of h 2 Y |X , is reiterated for different values of a time shift τ between x and y. The use of this time shift gives the information of "delay" between both signals.

h 2 * XY = max τmin<τ <τmax [h 2 Y |X (τ )], (2)
The values of h 2 * XY are comprised between 0 (Y is independent of X) and 1 (Y is linearly or nonlinearly dependent of X). In case of linear relationship between X and Y , h 2 * XY reduce to common linear correlation coefficient r 2 XY . Generally, h 2 * XY is not symmetric and △h 2 = h 2 * XY -h 2 * XY = 0. As described in [START_REF] Arnhold | A robust method for detecting interdependences : application to intracranially recorded eeg[END_REF], causality relations between X and Y can be deduced from △h 2 and τ giving and indication on coupling direction.

In theory, the energy flow between two parts of the brain, generated during epileptic seizure is then well characterized by its level (correlation) and direction (coupling direction). That would be right if the stationarity hypothesis was verified on X(t) and Y (t). In practice h 2 XY is computed on a 5s sliding window, as X(t) and Y (t) are considered quasi stationary, with an overlap of 0.5s. h 2 XY (t, τ ) was calculated for all values of delay τ included in an interval [-40 40]ms around current time t, and only the maximum value for h 2 XY (t, τ ) was kept (see fig. 1).

Signal coding

The goal of the method is the identification of reproducible epileptic patterns in h 2 XY signal. Before, these patterns may be highlighted by a similarity extraction method, described in next subsection, the first step is the encoding of data.

An important synchronization (respectively desynchronization) between temporal lobe structures is traduced on scalp EEG by a significatively high (respectively low) correlation coefficient h 2 * XY . The level of correlation is considered significatively important if the signal oversteps a predefined threshold Γ XYH . Its choice should answer to two major points:

• allowing the detection of high correlation coefficients corresponding to seizures • encode as less artifacts as possible

The definition of an absolute threshold for all channels does not suit because the level of correlation is highly dependent on channel pairs. Moreover, artefacts sometimes generate higher correlations than those observed at the seizure onset. To adapt to each channel, Γ XYH is chosen through analyse of its statistical distribution. The threshold Γ XYH is then defined as follows:

Γ XYH = h 2 XY (γ) ( 3
)
where γ is given by

γ i=1 N i N p = α% (4)
With N p the total number of samples of X or Y and N i the points number in each class i of the histogram of h 2 XY . Definition of the low level threshold Γ XYL is similar, the difference is in the choice of α. Once Γ XYH and Γ XYL are set, encoded signal S eXY (n) is defined as follows:

       S eXY (n) = a if h 2 XY (n) > Γ XYH , S eXY (n) = x or y if Γ XYH > h 2 XY (n) > Γ XYL , S eXY (n) = c if Γ XYL > h 2 XY (n).
(5) a, x, y and c are characters of the coding alphabet in which h 2 XY is transcript. In the description of an epileptic process, only high and low correlations are important. The extraction of similarities should then only occur on high and low coding levels. A different coding (x or y) is then needed in both seizures, in order not to extract similarities due to their preponderance in the h 2 XY coding matrices intermediate level.

Extraction of similarities

The goal of the similarity extraction between two seizures is to find epileptic patterns in coding matrices. As shown by [START_REF] Wendling | Extraction of spatio-temporal signatures from depth eeg seizure signals based on objective matching in warped vectorial observations[END_REF], a particularly well adapted algorithm is the Wagner and Fisher [START_REF] Wagner | The string to string correction problem[END_REF] extension to multichannel case. This algorithm allows the determination of the "distance between two strings of symbols as measured by the minimum cost sequence of edit operations needed to change one string into the other".

In the multichannel case, presented here, the operations are done on matrices, coded as sequence of vectors, instead of strings. Three edit operations are defined: insertion, deletion and substitution. All these operations have a certain cost chosen as follows:

• deletion and insertion: 1 • substitution: 2 if characters to substitute are different, 0 if they are the same

The edit distance is found through the computation of a cost matrix by a dynamic programming algorithm which iteratively estimate the lowest trace operation cost. A second algorithm, is then charged to extract the trace of minimal cost from the computed cost matrix. The characters that are found similar (substitution with cost 0) are outlined as shown in figure 2. A freedom degree may be let on extracted patterns as shown. For example, on fig. 2, pattern (4-5) lasts 2 samples in second matrix, but 5 samples in the first matrix.

Fig. 2. Example of similarity extraction by the use of Wagner and Fisher algorithm on simulated matrices. One can notice, that both matrices don't need to be the same length, and that some time distortion is allowed in subsequences recovering.

In figure 2 application of the Wagner and Fisher algorithm is presented on simulated matrices.

Only the costless result is shown, but some other possible sequence may be obtained with a change in the costs of the basic edit operations.

RESULTS

For all the results presented here, value used for α in equation 4 is chosen with the a priori information that temporal lobe seizure lasts for about 50s on the 500s of the file duration. To compute high threshold, α is set to 90% in order to keep only 10% of the highest correlations. As low correlations are rare threshold α is set to 5%.

In figure 3, result of Wagner and Fisher algorithm on two real seizures of a same patient is presented. Similarities can be extracted on eleven out of twenty-eight channels. Similarities are also extracted before the seizure but are caused by artefacts. What make artefacts and seizure different after similarity research is the length of the longest subsequence. Another important point, is that onset of the seizure indicated by neurologist does not exactly correspond to the same instant in both patterns extracted from each seizure.

To avoid subjective decision, concerning pattern recognition, one can use an objective index based previous remark. It is defined by the ratio of similarities observed in a seizure by those observed everywhere else defined as:

Sim index = channels

DISCUSSION

An important point to notice, is that correlation, previously applied to SEEG signal (Wendling et al., 2001a;[START_REF] Chavez | Spatiotemporal dynamics prior to neocortical seizurres : Amplitude versus phase couplings[END_REF] can be used on scalp EEG with good results. That was unsure because scalp EEG signal is the recording of a global information whereas SEEG is the recording of a focal activity.

The similarity method described here to characterize seizure onset proves to be efficient as more than 83% of the seizures presented a high level of similarities. Only 17% of seizures (15 out of 87) presented a really poor similarity index (less than one). Comparing these seizures with clinical files of the patients, one can notice that for the most part of them, they are either multifocal or bilateral.

In [START_REF] Caparos | Automatic laterlization of tle based on non-linear correlation on scalp eeg: Study on 43 patients[END_REF], authors showed that most important information on seizure were contained in its early seconds. Results obtained with a [-10,+20] second window seem to confirm it as longest windows give worst results. For example with a [-30,+30] s window only 60% of the seizures have a similarity index greater than 2. This percentage quickly fall to less than 50% with the augmentation of the window size. With [-150,+150]s window results obtained have no sense as practically all the file is used to compute index. This information let then the pattern discrimination be easier for short windows.

Results obtained for the measure of the edit distance computed intr-patients gives more difficult to analyse results, as mean of the edit distance is 0.95 (s.d. 0.03). Patient with a low edit distance (under mean-std) are not obviously patients that present a high Sim index . Actually, data with many artefacts have generally a low edit distance, that don't reflect what happens during the seizure. Result obtained for Wagner and Fisher algorithm reflect the similarity of the whole data file. To overcome this problem a more local edit distance should maybe be used. Smith and Waterman algo-rithm [START_REF] Smith | Identification of common molecular subsequences[END_REF], for example may be able to localize more local similarity zones.

As reproducible patterns can be extracted at the onset of seizures, the extraction of similarities may reinforce the medical diagnosis. This assumption is the base of the surgery that only works because the EZ and the neuronal network involved in several seizures are stable for a given patient. The next step of similarities study on scalp EEG is to realize a seizure classification based on the coupling of intra cerebral structure measured between surface electrodes.

  Fig. 1. Evolution of h 2 coefficient before and during seizure. Appearance of first epileptic given at 350s by a neurologist, on the figure one can see a sort of "spike" and then a sudden decrease in correlation. This is one of epileptic patterns observed.

Fig. 3 .

 3 Fig. 3. Result of the similarities extraction on encoded matrices (each line of a seizure represent a encoded h 2 XY ) by Wagner and Fisher algorithm on two seizures of the same patient. Black regions are those which were found the most similar in both seizures. Onset of the seizure indicated by neurologist is indicated by a vertical line.


Luis Fernando Camacho-Ortegón

Luis Martínez

Juan Josué Enciso-Cardenas

Arturo Bueno-Tokunaga

Jacques Pironon

Fernando Núñez-Useche

Thermal history of the Sabinas -Piedras Negras Basin (Northeastern Mexico): Insights from 1D modelling

Introduction

The Sabinas-Piedras Negras Basin (Coahuila, North-Eastern Mexico) is close to the Gulf of Mexico province (Fig. 1). The Sabinas Basin is considerate a maturity hydrocarbons basin, that gas-prone production [START_REF] Santamaría-Orozco | Evolución geodinámica de la Cuenca de Sabinas y sus implicaciones petroleras[END_REF][START_REF] González-García | Las rocas generadoras de México[END_REF]Holguín-Quiñones, 1992, 2001;[START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF][START_REF] Román-Ramos | Subsistemas generadores de la región norte de México[END_REF] with an original estimate of recoverable reserves of 766.8 BCF (Eguiluz de [START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF]. Despite the long history of exploration (in the 1930s) and production (the first productive well was drilled in 1975). Little is known about its thermal and burial history. Many hypotheses were established to explain high maturity levels of observed organic matter (OM) in this basin and its surroundings (e.g., Maverick Basin): hidden batholith-sized magmatic body inducing a regional enhanced heat flow, warm deep basin fluids expulsion [START_REF] Barker | The sacatosa coalbed methane field: a first for Texas[END_REF], or simple burial (Cuevas-Lerée and dissertation, 1985;[START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF] were discussed as possible causes.

Basin modelling has become an essential tool to study sedimentary basins' thermal and burial history and enables the verification of theories on basin evolution. For the Sabinas-Piedras Negras area, only a few studies used this approach, but only taking single wells into account [START_REF] Santamaría-Orozco | Evolución geodinámica de la Cuenca de Sabinas y sus implicaciones petroleras[END_REF]. In this paper, nine exploratory wells were studied using a basin modelling approach calibrating heat flow history and erosion to measure vitrinite reflectance as the main kinetic parameter of the maturity. The model results provide an explanation for the maturity levels reached. Hydrocarbon (HC) generation is also examined to determine the timing of hydrocarbon generation. The thermal modelling results, constrained by vitrinite reflectance, consider the maximum burial depth was reached during Paleocene -Eocene ages, taking into account that the youngest preserved rocks are of

The thermal and burial history of the Sabinas-Piedras Negras basin was studied by 1D modelling using vitrinite reflectance values measured on core samples from nine wells. Three heat flow scenarios led to model calibrations, including a simple one and a complex one integrating a thermal anomaly following the maximal burial phase. The differences in hydrocarbon generation curves due to these two scenarios were studied. The studied thermal anomaly induces two distinct phases of primary cracking of kerogen. The first phase is clearly due to burial, whereas the second seems to be due to the thermal event. The increase in the proportion of converted organic matter attributed to this second phase varies depending on the source rocks considered, e.g., for the La Casita source rock, there is a 10% TR increase, whereas that of the La Pe˜na and Eagle Ford source rocks, is 35 and 15% respectively. The thermal evolution of organic matter in the Sabinas-Piedras Negras Basin is mainly due to the burial phase. In contrast, the thermal event considered only causes a temporal shift of hydrocarbon generation. This hydrocarbon generation, primordially gas, occurred mainly before and during the compressive phase ("late Laramide" time step ) that affected the Sabinas-Piedras Negras Basin.

Maastrichtian age in the studied area. The main purposes of this study are : (a) this work aims to model the history of heat flow in the Sabinas Basin, associated with the thermal evolution of vitrinites, geodynamic events, and erosion; (b) investigate whether the evolution of heat flow in the basin influenced the second phase of cracking of organic matter; and (c) test our heat flow hypothesis by modelling the maturation of organic matter compared to what was observed within the Sabinas Basin.

Geological setting

The Sabinas-Piedras Negras Basin is located in the North-Eastern of Coahuila, Mexico (Fig. 1). Initially, it developed due to the opening of the proto Gulf of Mexico during the mid-Jurassic [START_REF] Alfonso | Geología regional del sistema sedimentario Cupido[END_REF][START_REF] Longoria | Stratigraphic studies in the Jurassic of northeastern Mexico: evidence for the origin of Sabinas Basin[END_REF][START_REF] Wilson | Basement structural controls on Mesozoic carbonate facies in Northeastern Mexico -a review[END_REF][START_REF] Santamaría-Orozco | Evolución geodinámica de la Cuenca de Sabinas y sus implicaciones petroleras[END_REF][START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF], which caused the formation of horsts and grabens, including the Sabinas -Piedras Negras basin. This basin may be considered as having evolved as a rift basin in the first part of its history. After this first phase, the Sabinas -Piedras Negras basin underwent decelerating subsidence [START_REF] Goldhammer | Middle Jurassic-Upper Cretaceous paleogeographic evolution and sequence stratigraphic framework of the Northwest Gulf of Mexico rim[END_REF]. During the Late Cretaceous and Early Tertiary, the basin is assumed to have behaved as a foreland basin during the Laramide orogeny [START_REF] Goldhammer | Middle Jurassic-Upper Cretaceous paleogeographic evolution and sequence stratigraphic framework of the Northwest Gulf of Mexico rim[END_REF][START_REF] Lawton | Stratigraphy and origin of the La Popa basin, Nuevo leon and Coahuila, Mexico[END_REF] before inversion. The sedimentary fill reaches a thickness of 5000 m and is mainly Mesozoic, and Santamaría-Orozco et al. (1991) reported a thickness of more than 6000 m in the depocenter. Three super sequences representing five depositional cycles were defined by [START_REF] Santamaría-Orozco | Evolución geodinámica de la Cuenca de Sabinas y sus implicaciones petroleras[END_REF], and adjusted by Eguiluz de [START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF] to describe the evolution of the sedimentary filling in the Sabinas-Piedras Negras Basin. The stratigraphic section for the Sabinas-Piedras Negras Basin is shown in Fig. 2.

The first supersequence is related to the initial rifting phase. It includes La Gloria and Olvido Formations, which were deposited during the Upper Jurassic, and it is mainly composed of conglomerates, salts, and evaporites (Fig. 2). The deposition of this supersequence (Cycle I according to Eguiluz de Antuñano, 2001) occurred principally in a sabkha and coastal marine environment during continuous subsidence [START_REF] Heim | Notes on the Jurassic of Tamazunchale Sierra Madre oriental, Mexico[END_REF][START_REF] Humphrey | Jurassic and Lower Cretaceous Stratigraphy and Tectonics of Northeast Mexico[END_REF][START_REF] Goldhammer | Sequence Stratigraphy and Cyclostratigraphy of the Mesozoic of Sierra Madre Oriental, Northeast Mexico: a Field Guidebook[END_REF][START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF].

The second supersequence, consisting of cycles II, III, and IV, according to Eguiluz de [START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF], is related to the thermal subsidence phase and is mainly composed of carbonates evaporites coastal siliciclastic materials deposited from Berriasian to Turonian ages (Fig. 2). According to Eguiluz de [START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF], this supersequence may be described by three depositional cycles: the first cycle (or cycle II) includes the La Casita and Menchaca -Taraises Formations. The second cycle (or cycle III) comprises the Barril Viejo, Padilla, La Mula, La Virgen, and Cupido Formations. The last depositional cycle (or cycle IV) includes the Cupidito, La Peña, Upper Tamaulipas, Kiamichi, George Town, Del Rίo, and Buda Formations. The deposition of this supersequence varies from deep marine, neritic, fan-delta, reef, deep basinal to open ramp platform environments on a passive margin [START_REF] Imlay | Geology of the middle part of the Sierra Parras[END_REF][START_REF] Márquez | Evaluación petrolera de sedimentos carbonatados del Cretácico en el Golfo de Sabinas, NE de México[END_REF][START_REF] Murillo-Muñetón | Stratigraphic Architecture, Platform Evolution, and Mud-Mound Development in the Lower Cupido Formation (Lower Cretaceous)[END_REF][START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF][START_REF] Goldhammer | Middle Jurassic-Upper Cretaceous paleogeographic evolution and sequence stratigraphic framework of the Northwest Gulf of Mexico rim[END_REF].

The last supersequence, or cycle V according to Eguiluz de Antuñano ( 2001), was deposited in a foreland setting from Turonian to Danian ages and is represented by regressive, terrigenous clastic facies. This supersequence consists of Eagle Ford, Austin, Parras, Upson, San Miguel, Olmos, Escondido Formations, Difunta -Navarro -Taylor, and Midway -Wilcox Groups. This upper age (Danian) is approximate and simply a reference. According to Eguiluz de Antuñano ( 2001), the 39.5 My B.P. unconformity represents a boundary sequence for this third supersequence, with deep erosion of pre-existing stratigraphic units. The deposits of this supersequence are characterized from ramp-slope, neritic to deltaic environments related to the initial pulses of orogenic uplift to the west [START_REF] Dumble | Notes on the geology of the valley of the middle Río Grande[END_REF][START_REF] Imlay | Geology of the western part of the Sierra de Parras[END_REF][START_REF] Kelly | Evolution of the Coahuila Peninsula, Mexico, part II, geology of the mountains bordering the valleys of acatita and las Delicias[END_REF][START_REF] Price | Stress orientations during Oligocene volcanism in Trans-Pecos, Texas: timing the transition from Laramide compression to basin and range tension[END_REF][START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF].

Structural setting

During the Laramide phase (Late Cretaceous/Early Tertiary), the region underwent an important compression phase, with a regional North-Eastward tectonic transport (Eguiluz de [START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF]. This phase formed a fold and thrust belt in the eastern region (the Sierra Madre Oriental (SMOr)) and compressive structures in the basin. In the northeastern part of the basin, the deformation was controlled by a high basement. Schematically, the Sabinas-Piedras Negras Basin is represented by narrow anticlines globally-oriented NW/SE, separated by large syncline valleys (Fig. 1). Where Jurassic salt is present, it behaves as a detachment level.

After the Laramide orogeny, a significant stress direction change occurred and decreased horizontal compression. This event could have been affected by the rollback associated with the Farallon Plate [START_REF] Bird | Stress direction history of the western United States and Mexico since 85 Ma[END_REF].

Stratigraphic formations of interest

Three main source rock intervals were identified in the Sabinas-Piedras Negras area according to the data supplied by PEMEX (Petroleos Mexicanos Company), and recent petrographic and geochemical studies [START_REF] Piedad-Sanchez | Prospection des hydrocarbures par une approche intégrée de pétrographie, géochimie et modélisation de la transformation de la matière organique: analyse et reconstitution de l'histoire thermique des Bassins Carbonifère Central des Asturies (Espagne) et Sabinas-Piedras Negras[END_REF][START_REF] Camacho-Ortegón | Origine-Evolution-Migration et Stockage, des hydrocarbures dans le bassin de Sabinas, NE Mexique: étude intégré de pétrographie, géochimie, géophysique et modélisation numérique 1D-2D et 3D[END_REF][START_REF] Camacho-Ortegón | Modelado del sistema petrolero de la Cuenca de Sabinas, Noreste de México ; Parte 1: evolución térmica, generación y migración de hidrocarburos[END_REF]. The most important one is the La Casita Formation (Kimmeridgian -Tithonian), with facies from a neritic to fan-delta complex environment, as was defined by [START_REF] Imlay | Geology of the middle part of the Sierra Parras[END_REF]), after Santamaría-Orozco et al. (1991) made the organic geochemistry analysis. Finally, Eguiluz de Antuñano (2001) and [START_REF] Goldhammer | Middle Jurassic-Upper Cretaceous paleogeographic evolution and sequence stratigraphic framework of the Northwest Gulf of Mexico rim[END_REF] abridged the geologic evolution of this formation. These sediments have high organic carbon content (Total Organic Carbon -TOC), which ranges between 0.5 and 3%, and high organic maturity (R r > 2.0%). The two other source rocks are La Peña Formation (Late Aptian; 0.4-2.9% TOC; deposited in a deep basinal environment; [START_REF] Imlay | Geology of the western part of the Sierra de Parras[END_REF][START_REF] Humphrey | Geology of the Sierra de los Muertos area, Mexico (with Descriptions of Aptian Cephalopods from the La Peña Formation)[END_REF][START_REF] Santamaría-Orozco | Evolución geodinámica de la Cuenca de Sabinas y sus implicaciones petroleras[END_REF][START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF][START_REF] Goldhammer | Middle Jurassic-Upper Cretaceous paleogeographic evolution and sequence stratigraphic framework of the Northwest Gulf of Mexico rim[END_REF], and Eagle Ford Formation (Turonian; 0.6-2.5% TOC; deposited in a neritic environment; [START_REF] Kelly | Evolution of the Coahuila Peninsula, Mexico, part II, geology of the mountains bordering the valleys of acatita and las Delicias[END_REF][START_REF] Santamaría-Orozco | Evolución geodinámica de la Cuenca de Sabinas y sus implicaciones petroleras[END_REF][START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF]; [START_REF] Goldhammer | Middle Jurassic-Upper Cretaceous paleogeographic evolution and sequence stratigraphic framework of the Northwest Gulf of Mexico rim[END_REF]. These three formations mainly contain terrigenous organic matter mixed with low proportions of marine organic matter. Locally, other formations (mainly La Virgen with a facies of sabkha-type lagoonal carbonate environment; [START_REF] Imlay | Neocomian Fauna of northern Mexico[END_REF][START_REF] Márquez | Evaluación petrolera de sedimentos carbonatados del Cretácico en el Golfo de Sabinas, NE de México[END_REF][START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF][START_REF] Goldhammer | Middle Jurassic-Upper Cretaceous paleogeographic evolution and sequence stratigraphic framework of the Northwest Gulf of Mexico rim[END_REF] may have significant organic matter content (reaching 3.6% TOC). They are then also considered as possible source rocks.

Several reservoir rocks were identified in the Sabinas-Piedras Negras Basin taking into account their geochemical and mechanical characteristics (fracturing, faulting, porosity, permeability, etc.), and their geostatistical modelling [START_REF] Márquez | Evaluación petrolera de sedimentos carbonatados del Cretácico en el Golfo de Sabinas, NE de México[END_REF][START_REF] González-García | Las rocas generadoras de México[END_REF]Holguín-Quiñones, 1992, 2001;[START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF][START_REF] Román-Ramos | Subsistemas generadores de la región norte de México[END_REF]. The main reservoir rocks are in the La Gloria (Upper Jurassic), La Casita (Upper Jurassic), Padilla (Hauterivian), La Virgen (Barremian), George Town (Albian), and Austin (Campanian) Formations [START_REF] Márquez | Evaluación petrolera de sedimentos carbonatados del Cretácico en el Golfo de Sabinas, NE de México[END_REF][START_REF] Santamaría-Orozco | Evolución geodinámica de la Cuenca de Sabinas y sus implicaciones petroleras[END_REF][START_REF] González-García | Las rocas generadoras de México[END_REF]Holguín-Quiñones, 1992, 2001;[START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF].

1D modelling

One-dimensional modelling of burial history and thermal maturity was performed on nine gas fields, primarily located in the central area of the basin. The wells' locations are shown in Fig. 3. Computations were performed using the PetroMod 1D software package. Burial history modelling provides the physical framework for thermal and hydrocarbon generation simulation (see a relationship between %R r and hydrocarbon generation, Table 1). The model was calibrated by modifying the heat flow and erosion until a satisfactory relationship between measured and simulated thermal indicators was observed. The calibration model included mainly the vitrinite reflectance and the wells' bottom temperature, paleobathymetry, and paleotemperature measurements. The software needs massive information asan input data, these are:

Input parameters

The lithological composition of the different formations was described as mixtures of pure mineralogical end-members. An example of the values considered in this study is given in Table 2 (Camacho-Ortegón, 2009). Refined end-members physical parameters are presented in Table 3.

Burial history constraints

The timing of maximum burial and beginning of inversion of the Sabinas-Piedras Negras Basin are not well defined. The youngest preserved sediments in the area of interest are of the Maastrichtian age. In this study, it is considered that the maximum burial depth was reached during Paleocene -Eocene ages and that erosion occurred from these periods until Early Neogene, as postulated by Eguiluz de Antuñano ( 2001), who used kinematic data and the critical-wedge mechanical model [START_REF] Dissertation | Mechanics of fold-and-thrusts belts and accretionary wedges: cohesive Coulomb theory[END_REF][START_REF] Davis | The role of salt in fold-and-thrust belts[END_REF][START_REF] Marret | Structure and kinematic development of the Sierra Madre Oriental fold-thrust belt, Mexico[END_REF]. These events are also confirmed by [START_REF] Gray | Thermal and chronological record of syn-to post-Laramide burial and exhumation, Sierra Madre Oriental, Mexico[END_REF], who utilized K/Ar ages on diagenetic illite, aqueous fluid inclusion homogenization temperatures, apatite fission-track, and U-Th/He periods. The same hypothesis was made when Román-Ramos and Holguín-Quiñones (2001) used biomarkers, vitrinite reflectance, and thermal kinetic modelling for this basin. The chronostratigraphic timescale used to date the different events is from Eguiluz de Antuñano ( 2001) and Cuevas-Lerée and dissertation (1985).

Extent of erosion

The most important erosional event occurred in the Early Tertiary. The extent of erosion in the Sabinas-Piedras Negras area is of significant concern because Late Cretaceous and Early Tertiary deposits are preserved only in the eastern part of the study area, reaching a thickness of 2400 m. Nevertheless, deposits and subsequent erosion should have been significantly more significant in the whole area to explain the high maturity levels observed in the basin. Cuevas-Lerée and dissertation (1985), [START_REF] Santamaría-Orozco | Evolución geodinámica de la Cuenca de Sabinas y sus implicaciones petroleras[END_REF], Román-Ramos and Holguín--Quiñones ( 2001) and Camacho-Ortegón (2009) using thermal modelling assumed that ca. 1500-1700 m of sediments were eroded in the Monclova area to explain the high maturity levels. The paleobathymetry data agrees with the sea-level fluctuations in Maastrichtian (Fig. 4; [START_REF] Vail | Seismic stratigraphy and global changes of sea-level, IV. Global cycles of relative changes of sea-level[END_REF][START_REF] Haq | Chronology of fluctuating sea levels since the Triassic[END_REF][START_REF] Camacho-Ortegón | Modelado del sistema petrolero de la Cuenca de Sabinas, Noreste de México ; Parte 1: evolución térmica, generación y migración de hidrocarburos[END_REF]; the calculated erosion is ca. 600-2350 m in the studied area. It seems thus that erosion may be regarded as an essential parameter to be taken into account during simulation. The missing Late Cretaceous -Early Tertiary deposits are assumed to be primarily terrigenous since, during this period, the shoreline migrated eastward [START_REF] Mcbride | Geologic report on upper Cretaceous coal-bearing rocks, Rio Escondido basin, Coahuila, Mexico[END_REF][START_REF] Goldhammer | Middle Jurassic-Upper Cretaceous paleogeographic evolution and sequence stratigraphic framework of the Northwest Gulf of Mexico rim[END_REF].

Thermal constraints

Surface temperatures were considered as having reached a value of ca. 25 • C during the Late Jurassic to Late Cretaceous times before decreasing to the present-day value of ca. 20 • C.

The first step of this thermal study was to define a heat flow evolution curve according to the evolution of the basin. The thermal regime refers to the data provided by [START_REF] Allen | Basin Analysis: Principles and Applications[END_REF], who estimate the [START_REF] Bird | Stress direction history of the western United States and Mexico since 85 Ma[END_REF], the average heat flow value was described referring to other studies. The average heat flow value during the foreland phase in the North Alpine Foreland basin is assumed to be close to 60 mW/m 2 by [START_REF] Schegg | Migration and accumulation of hydrocarbons in the Swiss Molasse Basin: implications of a 2D basin modelling study[END_REF] and lower than 55 mW/m 2 during the compressive phase by [START_REF] Sachsenhofer | Syn-and post-collisional heat flow in the Cenozoic eastern Alps[END_REF]. In a compressive setting, [START_REF] Gonçalves | Petroleum generation and migration in the Putumayo Basin, Colombia: insights from an organic geochemistry and basin modelling study in the foothills[END_REF] estimate heat flow values of ca 40 mW/m 2 . Considering the above, we infer that one value of 60 mW/m 2 is a maximum value for the Sabinas-Piedras Negras area during the foreland phase. These estimations lead to the evolution curve called scenario one (Fig. 5).

Thermal indicators

Thermal history reconstruction is based on random vitrinite reflectance measurements (R r ). R r remains the most commonly used thermal indicator in thermal modelling [START_REF] Handhal | Basin modeling analysis and organic maturation for selected wells from different oil fields, Southern Iraq[END_REF]. Vitrinite reflectance is a percentage measure of the incident light reflected from the surface of vitrinite particles in sedimentary rock [START_REF] Beaumont | Treatise of petroleum geology/handbook of petroleum geology: exploring for oil and gas traps. In: Chapter 6: evaluating source rocks[END_REF]. It is termed as %R o . It is optical parameters symbolized by VR or R o (reflectance in oil) [START_REF] Tissot | Petroleum Formation and Occurrence[END_REF]. Vitrinite is a coalification product of humic substances originating from the lignin and cellulose of plant cell walls [START_REF] Taylor | Organic Petrology[END_REF].

The relationship between %R r and hydrocarbon generation depends on the vitrinite's and the kerogen's chemistry [START_REF] Beaumont | Treatise of petroleum geology/handbook of petroleum geology: exploring for oil and gas traps. In: Chapter 6: evaluating source rocks[END_REF]. [START_REF] Sweeney | Evaluation of a simple model of vitrinite reflectance based on chemical kinetics[END_REF]) evaluated a simple model of R r (see EASY%Ro) based on chemical kinetics due to changes in vitrinite composition with time and temperature (Table 1).

The simplified model EASY%Ro uses an Arrhenius first-order parallel-reaction approach with a distribution of activation energies. With EASY%Ro, a vitrinite reflectance vs. time profile can be obtained for a given stratigraphic level if the time-temperature history for that level has been estimated. When applied to multiple stratigraphic levels, EASY %Ro can be used to compute profiles of the percent of vitrinite reflectance with depth for comparison with borehole data and optimize thermal history models. EASY%Ro successfully estimates vitrinite reflectance due to the thermal metamorphism of sedimentary rocks heated by igneous intrusions, geothermal fluids, and burial in various basin settings [START_REF] Sweeney | Evaluation of a simple model of vitrinite reflectance based on chemical kinetics[END_REF].

Measurements complemented the R r data provided by PEMEX and additional studies (Camacho-Ortegón, 2009) on samples from nine wells made following the established norm ISO 7404-5 (1994). When the organic matter particles were scarce, the recommendations proposed by [START_REF] Bertrand | Optical methods applied to source rock study[END_REF], [START_REF] Barker | Fluid-inclusion technique for determining maximum temperature in calcite and its comparison to the vitrinite reflectance geothermometer[END_REF], and [START_REF] Taylor | Organic Petrology[END_REF] were followed. Measurements were performed mainly on source rock samples (i.e., Eagle Ford, La Peña, La Virgen, and La Casita Formations). Measured R r values are similar to those given by [START_REF] Santamaría-Orozco | Evolución geodinámica de la Cuenca de Sabinas y sus implicaciones petroleras[END_REF], [START_REF] González-García | Las rocas generadoras de México[END_REF]Holguín-Quiñones (1992, 2001), and Román-Ramos and Holguín-Quiñones (2001), and give values higher than 0.6% for the Eagle Ford Formation, between 0.8 and 1.1% for La Peña, and higher than 2% for La Casita. Thermal history was reconstructed based on these parameters. The bottom temperature of some wells was used to compare the thermal history in 1D modelling.

Results

Present-day thermal regime

The present-day heat flow values were determined using the temperature profiles available for different wells. The discrepancy between the simulated profile and the measured temperatures is moderate (2-3 • C of difference). This result may be considered good [START_REF] Demongodin | Thermal conductivity and well logs: a case study in the Paris Basin[END_REF]. The present-day heat flow values obtained range between 50 and 68 mW/m 2 . These values support the heat flow values given by [START_REF] Bird | Deformation and uplift of north America in the Cenozoic era[END_REF] and [START_REF] Price | Stress orientations during Oligocene volcanism in Trans-Pecos, Texas: timing the transition from Laramide compression to basin and range tension[END_REF].

The present-day heat flow map (Fig. 6) established with the results obtained shows that heat flow is more significant in the center of the basin and decreases westward toward the margin of the basin. It also appears that the highest heat flow values are reached in the north of the basin.

Thermal history

R r profiles were simulated using the LLNL Easy%R o algorithm [START_REF] Sweeney | Evaluation of a simple model of vitrinite reflectance based on chemical kinetics[END_REF]. As a starting point, calculations were performed considering the hypotheses on erosion and heat flow variation during burial history (scenario 1; Fig. 5).

A good match is obtained for the PMX-129 well, considering a heat flow value of 53 mW/m 2 reached during the maximum burial phase. Nevertheless, for the other wells, it appears that the best matches are obtained considering heat flow values higher than 60 mW/m 2 at a full [START_REF] Vail | Seismic stratigraphy and global changes of sea-level, IV. Global cycles of relative changes of sea-level[END_REF][START_REF] Haq | Chronology of fluctuating sea levels since the Triassic[END_REF]modified after Camacho-Ortegón, 2009). Black line; paleobatimetry Sabinas-Piedras Negras Basin. Blu line; EXXON curve. Redline; Sea Level now (SLN). burial. Keeping the heat flow under 60 mW/m 2 leads to defective calibration (Fig. 7). It, therefore, appears necessary to project the heat flow value higher when the maximum burial phase occurs.

As for the initial heat flow scenario, assuming a maximum heat flow value of 60 mW/m 2 does not lead to a meaningful calibration with the wells, two additional heat flow scenarios were defined: a first one (scenario 2; Fig. 5), more straightforward than the first one, as it considers no heat flow limit in its variation during the evolution of the basin, and a second additional heat flow scenario (scenario 3; Fig. 5), slightly more complex than the first scenario: heat flow values are kept moderate (<60 mW/m 2 ) during the compressive phase, but a thermal anomaly is considered during the Tertiary [START_REF] Bird | Stress direction history of the western United States and Mexico since 85 Ma[END_REF][START_REF] Ewing | Review of the tectonic history of the lower Rio Grande border region, south Texas and Mexico, and implications for hydrocarbon exploration[END_REF][START_REF] Chávez-Cabello | Culminación de la Orogenia Laramide en la Cuenca de Sabinas, Coahuila, México[END_REF], after the beginning of the basin inversion. This thermal anomaly appears to be related to the intense magmatic activity observed in the Gulf of Mexico Province, known as the Eastern Alkaline Province (Bloomfield and Cepeda-Dávila, 1973). In the Sabinas area, three magmatic complexes are observed: the Nueva Rosita zone, the Monclova-Candela lineament and the Sierra de Picachos complex (see Fig. 3 for location). The origin of magmatism is currently interpreted in terms of two different tectonic environments: the Monclova-Candela lineament seems to be due to subduction tectonic [START_REF] Clark | Cenozoic mineral deposits and subductionrelated arcs in Mexico[END_REF][START_REF] Morton-Bermea | Petrologie, Mineralogie und Geochemie des Alkali-Intrusivkomplexes von Monclova-Candela (Mexiko)[END_REF][START_REF] Morton-Bermea | Magmatism of the northern Gulf of Mexico alkaline province (Sierra de Picachos and Monclova-Candela belt)[END_REF], whereas the rhyolitic nature of the Sierra de Picachos complex or other structures located close to the basin (see [START_REF] González-Partida | Fluorite deposits at Encantada-Buenavista, Mexico: products of Mississippi Valley type processes[END_REF] suggests a continental rift origin [START_REF] Morton-Bermea | Petrologie, Mineralogie und Geochemie des Alkali-Intrusivkomplexes von Monclova-Candela (Mexiko)[END_REF][START_REF] Morton-Bermea | Magmatism of the northern Gulf of Mexico alkaline province (Sierra de Picachos and Monclova-Candela belt)[END_REF]. The age of magmatism is estimated to ca 30 to 45 Ma for the subduction-related magmatism [START_REF] Clark | Cenozoic mineral deposits and subductionrelated arcs in Mexico[END_REF][START_REF] Price | Stress orientations during Oligocene volcanism in Trans-Pecos, Texas: timing the transition from Laramide compression to basin and range tension[END_REF][START_REF] Chávez-Cabello | Culminación de la Orogenia Laramide en la Cuenca de Sabinas, Coahuila, México[END_REF], but magmatism due to continental rifting appears to be younger (30-32 Ma, [START_REF] Price | Stress orientations during Oligocene volcanism in Trans-Pecos, Texas: timing the transition from Laramide compression to basin and range tension[END_REF] and temporally related to the extension phase already described (e.g., [START_REF] Bird | Stress direction history of the western United States and Mexico since 85 Ma[END_REF]. These geodynamical events could have induced a significant heat flow increase, explaining the maturity trends observed in the basin. In the model, this thermal event was considered at 38 Ma.

The results obtained for these two additional scenarios are similar, as no limit is given to the heat flow value (resulting calibration curves are presented in Figs. 7 and8). It appears quite difficult to rule out one of the two scenarios.

Effect of the uncertainty of different parameters on calibration curves

This section does not deal with the impact of the flow scenarios, as discussed earlier, but deals with other parameters, such as heat flow, erosion, or physical parameters (thermal conductivity).

Two tests were performed on the different wells, using the PetroRisk module integrated into PetroMod 1D Version 10. This module enables uncertainties to be defined and quantify these uncertainties' effect on the final results.

The first test concerns the effect of heat flow and erosion on the calibration curves. These parameters were increased and decreased in a ±10% interval. The resulting curves show an imperfect pattern for the minimum and maximum calibration curves (Fig. 9A), indicating thus that this interval may be considered the confidence interval.

Another test was performed to determine the effect of the thermal conductivity of the different formations on the calibration curves. For all the tested formations, thermal conductivity was considered to vary within a ±20% interval. The resulting curves (Fig. 9B) show that the effect of an essential uncertainty of thermal conductivity is negligible regarding the impact of heat flow and erosion variations.

This observation is confirmed by the third Test results, which integrate the three parameters types (erosion, heat flow, and thermal conductivity). The resulting curve (Fig. 9C) is similar to the curve observed when only heat flow and erosion variations are considered (Fig. 9A). Thus, it may be assumed that the values (erosion and heat flow) determined in this study have a margin of error lower than 10% and that the uncertainty of the physical parameters is of minor influence.

Discussion

Thermal and burial history

Regarding the distribution of the eroded thickness estimated from R r profiles (Fig. 8), it seems that the distribution of the deposits is more significant in the south (PMX-129 well, estimated original deposit thickness = 3500 m) than in the center of the basin, where estimated deposit thickness reaches 2100-2500 m. In the same way, an eastward decrease of the considered deposit amount in the calculations is In the south of the Sabinas -Piedras Negras basin (PMX-129 well), the difference between present-day temperature calculated assuming a present-day heat flow of ca. 60 mW/m 2 and maximum paleotemperature of burial (i.e., the highest temperature reached during burial) is ca. 120 • C. This value is slightly lower than values estimated from the curves obtained by [START_REF] Gray | Thermal and chronological record of syn-to post-Laramide burial and exhumation, Sierra Madre Oriental, Mexico[END_REF] for the La Popa Basin, with average values between 140 and 155 • C close to the front belt.

Therefore, it would appear that, in the Sabinas -Piedras Negras Basin, maximum erosion occurs close to the front belt of the SMOr and decreases with increasing distance from this limit. Fig. 10 shows the PetroMod 1D models by nine studied wells 49,62,65,68,97,109,112,and 129), and estimated burial history curve, overlay with oil, and gas window, isoline to temperature evolution, erosion, the timing of hydrocarbon generation in a systematic relationship by vitrinite reflectance maturity data (Fig. 8).

Lateral distribution of heat flow values and influence of salt structures

Heat flow anomalies calculated considering scenario three are mainly located in the northern part of the basin (Fig. 11), where most of the magmatic bodies were identified. Heat flow values associated with the thermal event vary from 100 to 75 mW/m 2 for wells, separated by a short distance, and therefore, they have a significant variation. Considering that some wells may be directly affected by heating due to the intrusions, these differences may be explained.

Salt domes were also observed in the Sabinas-Piedras Negras basin (Eguiluz de [START_REF] Eguiluz De Antuñano | Geologic evolution and gas resources of the Sabinas Basin in northeastern Mexico[END_REF]. Salt has a very high thermal conductivity and can modify the temperature distribution in sedimentary basins (e.g., [START_REF] Mello | Role of salt in restraining the maturation of subsalt rocks[END_REF]. Fig. 11 shows the distribution of salt domes and estimated heat flow during the thermal event. Note that there is no systematic relationship between thermal anomalies and the distribution of salt domes. This implies that salt does not influence observed vitrinite reflectance and temporal heat flow variations.

Effect of the thermal anomaly on hydrocarbon generation

Hydrocarbon generation was simulated in this study for the selected source rocks. Kerogen cracking may be simulated using different models [START_REF] Bouster | Study of the pyrolysis of polystyrenes: I. Kinetics of thermal decomposition[END_REF][START_REF] Braun | PMOD: a flexible model of oil and gas generation, cracking and expulsion[END_REF][START_REF] Burnham | An appropriate kinetic model for well-preserved algal kerogens[END_REF] such as T max -shift, serial first-order reaction, Bouster four-parameter, or three-parameter nucleation. However, it is commonly simulated assuming parallel reactions described with a single pre-exponential factor and a discrete or continuous distribution of activation energies (e.g., [START_REF] Ungerer | Extrapolation of kinetics of oil and gas formation from laboratory experiments to sedimentary basins[END_REF][START_REF] Pepper | Simple kinetic models of petroleum generation. Part I: oil and gas generation from kerogen[END_REF][START_REF] Vandenbroucke | Kinetic modelling of petroleum formation and cracking: implications from the high pressure/high temperature Elgin Field (UK, North Sea)[END_REF].

As shown previously, primary source rocks in this basin are mixtures of type II and III organic matter; currently available models deal with pure source rocks only, making it necessary to find a compromise. As type III organic matter is dominant in these three formations, only terrigenous organic matter is considered in this study. As no immature OM sample was available, simulations were performed using a default type III OM kinetic model because the La Casita Formation has this type of organic matter. The kinetic model used is the LLNL Type III model [START_REF] Burnham | A Simple Kinetic Model of Petroleum Formation and Cracking[END_REF]. This choice has no influence, as the aim of these simulations is not to evaluate the hydrocarbon potential of the basin but to test the relative effect of different heat flow scenarios on generation

curves.

Transformation ratio (TR) evolution was simulated considering the thermal field calculated using the two different heat flow scenarios considered possible (scenarios 2 and 3). These two thermal histories induce different shapes of TR evolution curves (Fig. 12).

CONSIDERING SCENARIO TWO, the HC generation curves show a continuous evolution curve, as the thermal field is mainly controlled by burial. Scenario three, integrating a thermal anomaly, induces a temporal shift of TR evolution curves. This delay is significantly more significant for the shallowest source rocks than the deepest (the La Casita Formation).

The thermal anomaly has another effect on TR evolution curves: it induces two distinct phases of primary cracking of kerogen. The first phase is clearly due to burial, whereas the second seems to be due to the thermal event. The increase in the proportion of converted OM attributed to this second phase varies depending on the source rocks considered. For the La Casita source rock, there is a 10% TR increase, whereas that of the La Peña and Eagle Ford source rocks is 35 and 15%, respectively.

This delay in HC generation due to the thermal event does not imply any considerable variation in the distribution of HC generated, as most of the generation occurs when maximum burial is reached. However, as explained earlier, the assumption that a thermal event induced the secondary cracking of primary components entirely agrees with petrographic observations. The generation of gas by secondary cracking should be considered when estimating the Sabinas-Piedras Negras gas potential.

Based on these results, it appears that hydrocarbon generation occurred mainly before and during the compressive phase that affected the Sabinas -Piedras Negras basin. Assuming that the folding occurred during this phase of compression and taking into account the reported data by [START_REF] Santamaría-Orozco | Evolución geodinámica de la Cuenca de Sabinas y sus implicaciones petroleras[END_REF][START_REF] González-García | Las rocas generadoras de México[END_REF], 2001), Eguiluz de Antuñano ( 2001), [START_REF] Gray | Thermal and chronological record of syn-to post-Laramide burial and exhumation, Sierra Madre Oriental, Mexico[END_REF][START_REF] Román-Ramos | Subsistemas generadores de la región norte de México[END_REF][START_REF] Román-Ramos | Subsistemas generadores de la región norte de México[END_REF], hydrocarbon traps formed while most of the hydrocarbons were generated (Fig. 13). Nevertheless, a complete study would be necessary to determine if the hypothesis of a thermal event could induce differences in modelling integrating hydrocarbon migration.

Preferred heat flow scenario

From the results obtained considering either scenario two or three, it appears challenging to rule out one of the two scenarios. Nevertheless, tectonic events such as subduction and extension are expected to induce significant increases in local heat flow, more so if magmatism is associated. This would be the first argument in favor of scenario 3.

Another argument is carried out by the petrographic study of the basin [START_REF] Piedad-Sanchez | Prospection des hydrocarbures par une approche intégrée de pétrographie, géochimie et modélisation de la transformation de la matière organique: analyse et reconstitution de l'histoire thermique des Bassins Carbonifère Central des Asturies (Espagne) et Sabinas-Piedras Negras[END_REF][START_REF] Camacho-Ortegón | Origine-Evolution-Migration et Stockage, des hydrocarbures dans le bassin de Sabinas, NE Mexique: étude intégré de pétrographie, géochimie, géophysique et modélisation numérique 1D-2D et 3D[END_REF]. Different types of bitumen and pyrobitumen were observed during the petrographic analysis of the Sabinas-Piedras Negras basin. Devolatilization vacuoles, porosity, and fracturing were also observed in macerals and could be interpreted as evidence of secondary cracking [START_REF] Van Balen | Modelling the hydrocarbon generation and migration in the West Netherlands Basin, The Netherlands[END_REF][START_REF] Zeman | Migration of fluids controlled by equidistant fracture systems: an example from Central Europe (Czech Republic, Slovakia and Austria)[END_REF][START_REF] Hasiah | Suberinite -bituminite and phlobaphinite -inertinite associations in the oil -generating coals of Northwest Borneo[END_REF][START_REF] Martinez | Evolution thermique des séries sédimentaires de la Marge ardéchoise: étude pétrographique de la matière organique et des argiles (Thermal evolution of the sedimentary series of the Ardèche Margin: petrological studies of organic matter and clays)[END_REF]. Different thermal events explained these observations. USING SCENARIO THREE, the HC generation was simulated and reconstructed oil and gas generation from a source rock, migration, and accumulation with PetroMod 2D package (see the tutorial by modelling methods with Flowpath, Darcy, Invasion Percolation, and Hybrid (Darcy + Flowpath) technology [Petromod 2D tutorial, 2007]), on the seismic line A-A ′ for the Sabinas-Piedras Negras Basin (Fig. 13A). Cautiously, the modelling results confirm the studied basin like a gas-prone basin, showing mainly four gas accumulations in regional 2D cross-section (results from the 2D modelling to section A-A') equivalent of ~154 Mm 3 of dry gas (Fig. 13B). These cautious predictions agree with the defined four gas reservoirs (Fig. 13B; 1 to 4), mainly in the PMX-20 well with a reasonable gas production rate.

Relationship between %Rr and maximum fluid inclusions temperature

The fluid inclusion microthermometric analyses have been used successfully in the petroleum exploration industry (e.g., [START_REF] Bodnar | Interpretation of microthermometric data for H2O-NaCl fluid inclusions[END_REF]Pironon, 2004 and others) for estimating the temperatures, timing, and pore fluid origins of various diagenetic minerals in sources and reservoir rocks (e.g., [START_REF] Roedder | Fluid inclusions[END_REF][START_REF] Walderhaug | A fluid inclusion study of quartz-cemented sandstones from offshore mid-Norway-possible evidence for continued quartz cementation during oil emplacement[END_REF]. [START_REF] Tobin | Multidisciplinary thermal maturity studies using vitrinite reflectance and fluid inclusion microthermometry: a new calibration of old techniques[END_REF] established fluid inclusion microthermometry from a new inorganic thermal maturity tool that can help fill the void in geologic and sampling limitations. This tool uses an empirical calibration of fluid inclusion data and vitrinite reflectance data to estimate thermal maturity (Fig. 14). Currently, the models that relate the %Rr with the homogenization temperature (Th) of fluid inclusions allow the model to approach a more acceptable thermal calibration [START_REF] Barker | The correlation of vitrinite reflectance with maximum temperature in humic organic matter[END_REF][START_REF] Barker | Fluid-inclusion technique for determining maximum temperature in calcite and its comparison to the vitrinite reflectance geothermometer[END_REF]. This correlation technique [START_REF] Barker | The correlation of vitrinite reflectance with maximum temperature in humic organic matter[END_REF][START_REF] Tobin | Multidisciplinary thermal maturity studies using vitrinite reflectance and fluid inclusion microthermometry: a new calibration of old techniques[END_REF] is beneficial in basins with no macerals [START_REF] Bourdet | Petroleum accumulation and leakage in a deeply buried carbonate reservoir, Níspero field (Mexico)[END_REF][START_REF] Blaise | Reconstruction of low temperature (<100 • C) burial in sedimentary basins: A comparison of geothermometer in the intracontinental Paris Basin[END_REF], which helps users to accurately convert Th data to vitrinite reflectance equivalents (%Rr equivalent) in thermal maturity studies.

In conclusion, the latest techniques conventional in organic maturity (e.g., [START_REF] Sweeney | Evaluation of a simple model of vitrinite reflectance based on chemical kinetics[END_REF]ISO 7404-3, 1994;ISO 7404-5, 1994) have proven to be effective in determining the evolution of diagenesis in sedimentary series and optimizing basin models in the last 50 years.

Finally, Sabinas Piedras Negras Basin is suitable for applying the [START_REF] Tobin | Multidisciplinary thermal maturity studies using vitrinite reflectance and fluid inclusion microthermometry: a new calibration of old techniques[END_REF] tool since it would allow a more detailed evaluation of the evolution of the heat flux in a larger area than that represented in this work, considering a new Th study in outcrops from the Popa Basin to the Alto de Tamaulipas area, contributing further information to the regional tectonic, thermal anomaly and diagenetic history.

Conclusion

During the basin modelling study performed on the Sabinas Piedras Negras basin, three heat flow scenarios led to similar calibration results, and it appears difficult to rule out either one.

Hydrocarbon generation curves were simulated considering these two possible scenarios. The differences between the two curves show that this thermal event would not have induced a critical secondary cracking phase. Most organic matter evolution would have occurred when full burial was reached. Geodynamic considerations and, more importantly, a previous petrographic study would favor the hypothesis of a thermal anomaly occurring after the burial phase.

The results of this modelling approach for the Sabinas-Piedras Negras basin confirm an OM with a high maturity level. Apparently, in the north of the basin, a thermal anomaly can be delineated and should be considered for the HC evaluation in the basin.

Considering the constraints for modelling the regional 2D crosssection of the study area showing the simulation of HC generation for a seismic line indicates potential accumulation to 154 Mm 3 of dry gas. The gas accumulations coincide with some defined reservoirs through the gas producer.

Each exploratory well may calculate the local erosion and heat flow to improve the results of this modelling in the Sabinas-Piedras Negras Basin. New geological data of other wells (e.g., fluid inclusions, 3D seismic, isotopes, and other geochemical and petrographic data.) in the basin are indispensable to locate new reservoirs and quantify HC reserves with the best precision in this prone gas basin.

CRediT authorship contribution statement

Luis Fernando Camacho-Ortegón:Metodology, Investigation. Luis Martínez: Software. Juan Josué Enciso-Cardenas: Investigation. Arturo Bueno-Tokunaga: Data curation. Jacques Pironon: Investigation. Fernando Núñez-Useche: Writingreview & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .

 1 Fig. 1. Location and simplified structural map of the Sabinas -Piedras Negras Basin area. Discontinuous square: location of Fig. 3.

Fig. 2 .

 2 Fig. 2. Stratigraphic section of the Sabinas-Piedras Negras area (Modified from Santamaría-Orozco et al. (1991), and Eguiluz de Antuñano (2001)).

Fig. 3 .

 3 Fig. 3. Location of the studied wells, of the magmatic bodies and salt domes in the Sabinas-Piedras Negras Basin.

Fig. 4 .

 4 Fig. 4. Chronology of fluctuating Sea Level events in the Gulf of Sabinas and paleobathymetry in the Sabinas-Piedras Negras Basin (EXXON curve data; Vail[START_REF] Vail | Seismic stratigraphy and global changes of sea-level, IV. Global cycles of relative changes of sea-level[END_REF][START_REF] Haq | Chronology of fluctuating sea levels since the Triassic[END_REF] modified after Camacho-Ortegón, 2009). Black line; paleobatimetry Sabinas-Piedras Negras Basin. Blu line; EXXON curve. Redline; Sea Level now (SLN).

Fig. 5 .

 5 Fig. 5. Heat flow scenarios considered in the thermal reconstruction of the Sabinas-Piedras Negras Basin.

Fig. 6 .

 6 Fig. 6. Present-day heat flow distribution across the basin. Legend is given in Fig. 1 and 3.

Fig. 7 .

 7 Fig. 7. Example of R r profiles simulated considering scenarios and measured values for two wells.

  observed.

Fig. 8 .

 8 Fig. 8. Calibration results for the different wells. See Fig. 3 for the location of the wells.

Fig. 9 .

 9 Fig. 9. Calibration interval obtained considering margin error on heat flow and erosion (A) and thermal conductivity (B), and the three parameters (C).

Fig. 10 .

 10 Fig. 10. Burial history for nine well's in the Sabinas-Piedras Negras Basin. The results mainly indicate a gas-prone potential agreeing with the current HC production.

Fig. 11 .

 11 Fig. 11. Calculated heat flow values associated with the thermal event. Legend is given in Figs. 1 and 3.

Fig. 12 .

 12 Fig.12. The vertical scale of the Eagle Ford Formation, whose bottom does not surpass the 1000 m in depth at the Sabinas-Piedras Negras Basin, is exaggerated to show the TR trend. Considering scenarios 2 (continuous line) and 3 (dashed line), transformation ratio evolution curves were calculated for the three primary sources of rocks. Black arrows: thermal anomaly.

Fig. 13 .

 13 Fig.13. Cross-section on a seismic line in the Sabinas-Piedras Negras Basin (A). Location of the cross-section in Fig.3 (A-A ′). The simulation of the HC generation with PetroMod shows mainly four gas accumulations in section A-A ′ with a total of ~154 Mm 3 of dry gas. Interpreted section for seismic sequence analysis in the study area (B). Currently, the PMX-20, PMX-49, and PMX-94 wells have gas production. Distances and depth are in meters.

Fig. 14 .

 14 Fig. 14. Logarithmic Crossplot of vitrinite reflectance data (% Ro) vs. homogenization temperature of fluid inclusion Th ( • C), a new inorganic thermal maturity tool for estimating thermal maturity values (VR equivalents) from fluid inclusion data that capitalizes on the tendency of aqueous inclusions in calcite to reequilibrate during burial overheating. (modified after Tobin and Claxton, 2000).

Table 1

 1 Easy vitrinite reflectance after[START_REF] Sweeney | Evaluation of a simple model of vitrinite reflectance based on chemical kinetics[END_REF].

	%Rr	HC generation
	0.25-0.55 0.55-0.70 0.7-1.00 1.00-1.30 1.30-2.00 2.00-4.00 >4.0	Immature Early oil Main oil Late oil Wet gas Dry gas Over mature

Table 2

 2 Example of lithological composition considered in this study (after Camacho-Ortegón, 2009).

Table 3

 3 Main physical parameters considered for modelling.

	Lithotype	Density (kg/m3)	Initial porosity (%)	Compressibility		Thermal conductivity	Heat capacity
				(Pa-1)		(W/m K)		(cal/g K)
				Max	Min	0 • C	100 • C	20 • C	100 • C
	Water Shale Siltstone Sandstone Limestone Dolomite Marl Coal Evaporite Organic matter		0 65 56 42 42 30 47 52 10 0	2 60,000 8000 500 300 250 940 130,000 60,000 50,000	1 10 10 10 25 10 10 10 10 10	0,6 1,98 2,14 3,12 2,83 3,81 2,23 0,5 4,69 0,5	0,68 1,91 2,03 2,64 2,56 3,21 2,11 0,46 3,91 0,45	0,999 0,213 0,201 0,178 0,195 0,202 0,208 0,204 0,194 0,200	1008 0258 0242 0209 0223 0229 0248 0248 0210 0240

Acknowledgments

The authors express their gratitude to Exploration and Production -North Region -Petróleos Mexicanos (PEMEX, Mexico), for providing data and permission for publishing this manuscript. The authors are also indebted to the Centro de Investigación en Geociencias Aplicadas (CIGA-UAdeC, México), INCAR-CSIC (Oviedo, Spain), and GeoRessources lab Université de Lorraine (Nancy, France) for help with the analyses. The authors are also indebted to the IES, Germany support team for helpful technical support. R.F. Sachsenhofer is gratefully acknowledged for his helpful comments. The authors offer special thanks to the reviewers of this article, especially to Ph. D. Francisco J. Vega and Ph.D. Claudio Bartolini for his suggestions.