Benjamin Bruneau
email: benjamin.bruneau@upmc.fr

Benoit Chauveau

François Baudin

Isabelle Moretti

3D stratigraphic forward numerical modelling approach for prediction of organic-rich deposits and their heterogeneities

Keywords: Petroleum Source Rock, Paris Basin, Liassic, Basin Modelling

teaching and research institutions in France or abroad, or from public or private research centers.

3D stratigraphic forward numerical modelling approach

for prediction of organic-rich deposits and their heterogeneities

Introduction

Several organic-rich layers were deposited in the north-western European basins during the Early Jurassic due to their key position in an epicontinental sea at the north-western tip of the Tethys Ocean (e.g. Fleet et al., 1988;Korte and Hesselbo, 2011;Silva et al., 2011) (Fig. 1).

Since the area was rather flat, the size and geometry of the basins were strongly dependent on sea-level variations. While only a few areas remained continuously above sea level, others emerged from time to time, leading to temporary isolation of the north-western European seas from the open ocean or to the development of thresholds. The European domain was located at the western extremity of the domain affected by opening of the Tethys Ocean as well at the north-eastern margin of the future Atlantic Ocean between North America and Africa/northwestern Europe. Although this domain was only weakly affected by rifting, some small extensional faults are observed bordering the depocentres (Fig. 1). The transgressive phase of the Lower Jurassic occurred during sedimentation of the organic-rich deposits. Pre-existing crustal heterogeneities (related to age of the lithosphere and faulting) induced differential responses to the sedimentary loading [START_REF] Guillocheau | Mesocenozoic geodynamic evolution of the Paris Basin: 3D stratigraphic constraints[END_REF]. As a result, the bathymetry, connections and thresholds of the European basins evolved rapidly during the transgression, [START_REF] Baudin | Mapping and geochemical characterization of the Toarcian organic matter in the Mediterranean Tethys and Middle East[END_REF]. b. Paleogeography of the Upper Sinemurian (modified after Thierry et al., 2010). Red and green dots indicate the studied sections and wells.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Many studies have been carried out on Toarcian source-rocks in north-west Europe (e.g. Fleet et al., 1988;[START_REF] Bessereau | Distribution of organic matter in the Liassic series of the Paris basin: an exemple of organic heterogeneity in a source rock interval. In: Generation, accumulation and production of Europe's hydrocarbons II[END_REF], but less research has been focused on the several other organic-rich deposits of the Lower Jurassic. In the case of the Paris Basin, where the petroleum system is relatively prolific, a synthesis of the well data allows an overview of the thicknesses and richness of the three main organic-rich intervals [START_REF] Bessereau | Source rock occurrence in a sequence stratigraphic framework: the example of the Lias of the Paris Basin. In: Paleogeography, Paleoclimate, ands Source Rocks[END_REF][START_REF] A N U S C R I P T A C C E P T E D Accepted Manuscript 51 Delmas | Rapport regional d'évaluation pétrolière -Paris Basin Petroleum Potential[END_REF]Fig. 2). The thicknesses are computed from the cumulative thicknesses of intervals with TOC > 1 wt. % (modified after [START_REF] Bessereau | Distribution of organic matter in the Liassic series of the Paris basin: an exemple of organic heterogeneity in a source rock interval. In: Generation, accumulation and production of Europe's hydrocarbons II[END_REF][START_REF] A N U S C R I P T A C C E P T E D Accepted Manuscript 51 Delmas | Rapport regional d'évaluation pétrolière -Paris Basin Petroleum Potential[END_REF].

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

The spatial and temporal distributions and characteristics such as the Total Organic Content (TOC) of these organic-rich source rocks, result from a complex interaction between the processes that lead to organic matter accumulation in sediments, i.e.: primary production, transport, sedimentation, burial and degradation of organic matter. Previously, these processes were studied using specific approaches, by taking into account empirical equations or physical laws defined in the literature (e.g. Martin et al., 1987;[START_REF] Betts | The oxygen content of ocean bottom waters, the burial efficiency of organic carbon, and the regulation of atmospheric oxygen[END_REF]Mann and Zweigel, 2008). Tyson (1996) proposed a theoretical model including all these processes (Fig. 3). Opposite models are those which assume that processes associated with anoxia lead to organic-rich deposits, and those where primary productivity is the main factor explaining high organic matter contents. In these conceptual models, different factors could be involved such as the presence of sills (e.g. [START_REF] Pompeckj | Der Jura zwischen Regensburg und Regenstauf[END_REF][START_REF] Röhl | Lower Toarcian (upper liassic) black shales of the central European epicontinental basin: a sequence stratigraphic case study from the SW German Posidonia shale[END_REF], basin topography (e.g. Hallam and Bradshaw, 1979;[START_REF] Wignall | Model for transgressive black shales?[END_REF][START_REF] Röhl | Lower Toarcian (upper liassic) black shales of the central European epicontinental basin: a sequence stratigraphic case study from the SW German Posidonia shale[END_REF], variations of sea level and primary productivity (e.g. [START_REF] Wignall | Model for transgressive black shales?[END_REF][START_REF] Röhl | Lower Toarcian (upper liassic) black shales of the central European epicontinental basin: a sequence stratigraphic case study from the SW German Posidonia shale[END_REF] and oceanic anoxic events (e.g. Schlanger and Jenkyns, 1976). Although the conceptual model has been improved over recent decades, highlighting the relationships between the processes involved and empirical equations, there have been only few advanced numerical modelling studies (e.g. Mann and Zweigel, 2008;[START_REF] Bruneau | 3D Numerical Modelling of Marine Organic Matter Distribution: example of the Early Jurassic sequences of the Lusitanian Basin (Portugal)[END_REF]. Thus, despite all the research carried out to date, the history of organic matter deposition processes

M A N U S C R I P T A C C E P T E D
need to be integrated into a 3D forward model to allow the identification and quantification of the key factors responsible for the heterogeneity of organic-rich deposits.

For this purpose, we focus on north-western Europe where four main basins developed during the Early Jurassic (Paris Basin, Aquitaine Basin and South-East Basin) separated by thresholds (e.g. Poitou sill, Bresse-Jura area) from three main emerged land masses (London-Brabant Massif, Armorican Massif and Central Massif) (Fig. 1). This case study allows us to define various conditions for the deposition of organic-rich sediments. A 3D stratigraphic forward model of this area (1,000×1,000 km) from the Hettangian to the Toarcian was computed with a 0.1-Myr time step. We used the DionisosFlow© software [START_REF] Granjeon | Modelisation stratigraphique deterministe -conception et applications d'un modele diffusif 3D multilithologique[END_REF][START_REF] Granjeon | 3D stratigraphic modeling of sedimentary basins[END_REF][START_REF] Granjeon | Concepts and applications of a 3-D multiple lithology, diffusive model in stratigraphic modeling[END_REF][START_REF] Granjeon | 3D stratigraphic modeling in complex tectonics area[END_REF] and its new module, which includes the simulation of organic matter contents [START_REF] Granjeon | Sedimentary basin development method using stratigraphic simulation coupled with an organic matter production and degradation model[END_REF][START_REF] Bruneau | 3D Numerical Modelling of Marine Organic Matter Distribution: example of the Early Jurassic sequences of the Lusitanian Basin (Portugal)[END_REF]. The validity of the model was ensured with data compiled from IFP Energies Nouvelles and ENGIE legacy databases containing field and well studies and geochemical/mineralogical analyses.

M A N U S C R I P T A C C E P T E D

Fig. 3. Main processes affecting organic matter history, as simulated in DionisosFlow©

Organic-Rich Sediment (modified after Tyson, 1995).

The global data set shows that organic-rich deposits are not equally distributed in time. Three main periods appear to account for the deposition of about 80% of organic-rich sediments giving rise to marine source rocks (end-Devonian, Jurassic and Late Cretaceous). Global oceanic anoxic events have been proposed to explain some of these synchronous deposits (e.g. Schlanger and Jenkyns, 1976). However, as discussed by these authors (op. cit.), when using a higher resolution for the time frame, the source rocks do not appear quite so synchronous and other parameters influence the deposition and preservation of organic matter.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Geological setting

Structural and stratigraphic setting

The hinged position of Europe during the opening of the Tethys and Atlantic oceans led its structural development (Fig. 1). Many horsts and grabens define areas of subsiding basins and areas of structural highs or emerged lands (Fig. 1) (e.g. [START_REF] Bernoulli | Alpine, Mediterranean and central Atlantic Mesozoic facies in relation to the early evolution of the Tethys. In: Modern and Ancient Geosynclinal Sedimentation[END_REF][START_REF] Bassoulet | Le Toarcien inférieur: un épisode de crise dans les bassins et sur les plates-formes carbonatées de l'Europe du Nord-Ouest et de la Téthys[END_REF][START_REF] Dromart | Variation cyclique de la production carbonatee au Jurassique le long d'un transect Bourgogne-Ardèche, Est-France[END_REF]. The structural setting induced differential responses to the sedimentary loading and the soft extensional tectonic context led to the reactivation of some faults which localize the subsidence (e.g. [START_REF] Wilson | The Lusitanian Basin of westcentral Portugal: mesozoic and Tertiary tectonic, stratigraphic and subsidence history. In: Extensional Tectonics and Stratigraphy of the North Atlantic Margin[END_REF]Kullberg et al., 2001). Although some lower order cycles are defined (e.g. [START_REF] Guillocheau | Mesocenozoic geodynamic evolution of the Paris Basin: 3D stratigraphic constraints[END_REF] that are often specific to each basin, two second-order cycles can be observed in most of the basins (e.g. Hardenbol et al., 1998;[START_REF] Guillocheau | Mesocenozoic geodynamic evolution of the Paris Basin: 3D stratigraphic constraints[END_REF] (Fig. 4).

The first cycle corresponds to a regressive episode at the end of the Pliensbachian, coming to an end at the Pliensbachian -Toarcian boundary. This episode is marked by a non-deposition in the proximal parts and limestone deposits in the distal parts of the basins. The transgressive second-order cycle following this regressive episode marks the maximum flooding of the Liassic transgression and is observable throughout Europe in the lower Toarcian. This surface marks the maximum flooding of European platforms and the beginning of the Toarcian anoxic episode (T-OAE) (e.g. Hesselbo et al., 2007). A further regressive phase recorded by oolitic deposits across Europe marks the end of the Early Jurassic.

Toarcian oceanic anoxic event

The Lower Toarcian sediments of north-western Europe attain high TOC values (up to 25 wt. %), in common with other organic-rich deposits observed worldwide. Jenkyns (1988) suggested an episode of worldwide anoxic conditions named the Toarcian Anoxic Oceanic Event (T-OAE). Numerous studies have provided constraints on the timing of the T-OAE , proposing a duration of around 200 kyrs to 1 Myrs (Jenkyns, 1988;[START_REF] Mcarthur | Strontium isotope profile of the early Toarcian (Jurassic) oceanic anoxic event, the duration of ammonite biozones, and belemnite palaeotemperatures[END_REF][START_REF] Bailey | Paleoceanographic changes of the Late Pliensbachian-Early Toarcian interval: a possible link to the genesis of an Oceanic Anoxic Event[END_REF]Mattioli et al., 2004a;Kemp et al., 2005;[START_REF] Boulila | Astronomical calibration of the Toarcian Stage: Implications for sequence stratigraphy and duration of the early Toarcian OAE[END_REF] and a synchronous development at the European scale according to the biostratigraphy (Mattioli et al., 2004a and2004b;[START_REF] Van Breugel | Isorenieratane record in black shales from the Paris Basin, France: Constraints on recycling of respired CO2 as a mechanism for negative carbon isotope shifts during the Toarcian oceanic anoxic event[END_REF], except for the final stages of this event [START_REF] Baudin | Caractérisation géochimique et sédimentologique de la matière organique du Toarcien téthysien (Méditerrannée, Moyen-Orient), significations paléogéographiques[END_REF]Jenkyns et al., 2002;Mailliot, 2006;Hesselbo et al., 2007). However, the T-OAE shows a very particular context: the accumulation of organic-rich deposits occurred in epicontinental basins in Europe. Even if the T-OAE is well identified in Europe, there is a marked variability and heterogeneity of the organic matter content in sediments from one basin to another [START_REF] Baudin | Caractérisation géochimique et sédimentologique de la matière organique du Toarcien téthysien (Méditerrannée, Moyen-Orient), significations paléogéographiques[END_REF][START_REF] Baudin | Mapping and geochemical characterization of the Toarcian organic matter in the Mediterranean Tethys and Middle East[END_REF][START_REF] Bassoulet | Le Toarcien inférieur: un épisode de crise dans les bassins et sur les plates-formes carbonatées de l'Europe du Nord-Ouest et de la Téthys[END_REF]. Many models have led to hypotheses on the local and global factors explaining the T-OAE as well as other more localized anoxic events that may have occurred in Europe during the Early Jurassic (e.g. Fleet et al., 1988).

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 12 2.

Biostratigraphy

Many studies establish a detailed biostratigraphic zonation with precise resolution for this time interval (e.g. Domergues et al., 2004;[START_REF] Elmi | Toarcian Working Group 2005 Report and prospects[END_REF][START_REF] Comas-Rengifo | El Sinemuriense Superior (cronozonas Oxynotum y Raricostatum) en la región de S. Pedro de Moel (Cuenca Lusitánica, Portugal). Fácies carbonatadas ricas em matéria orgânica do Jurássico da Bacia Lusitânica[END_REF][START_REF] Duarte | Carbon isotope stratigraphy and ammonite biochronostratigraphy across the Sinemurian-Pliensbachian boundary in the western Iberian margin[END_REF].

Hence in this study, chronozone intervals are determined according to the biostratigraphy of the International Commission on Stratigraphy and from recent studies or revisions of the biochronostratigraphy (see [START_REF] Ogg | Jurassic. In: The Geologic Time Scale[END_REF]) (Fig. 5).

Sedimentary deposition model

The facies distribution of the Lower Jurassic is influenced by eustatic variations and local fault activity. The shallow marine and confined depositional environments with dolomitic limestones and evaporites are replaced by carbonate ramp facies which back-step onto emerged lands. At the scale of the whole of France, the facies are mainly composed of marl-limestone alternations. In this paper, the depositional model is built from our dataset combined with field studies (Fig. 6). The subdivision of environmental settings used here is derived from [START_REF] Burchette | Carbonate ramp depositional systems[END_REF]. Fig. 6. Carbonate ramp settings (modified after [START_REF] Burchette | Carbonate ramp depositional systems[END_REF]

Characteristics of organic-rich facies

The record of organic-rich facies is specific to each basin. The organic matter is mainly of marine type with some local influence from terrestrial inputs. In general, the organic-rich facies are attributed to episodes where marly sedimentation takes place in parallel with a substantial production of phytoplankton, in settings favourable for the development of anoxic environments. These deposits are often attributed to transgressive periods in the basins. The TOC values and mineralogical compositions within the same interval can be very variable (e.g. [START_REF] Espitalié | Géochimie organique du Bassin de Paris[END_REF]Hollander et al., 1991;Hanzo and Espitalié, 1993;[START_REF] Bessereau | Source rock occurrence in a sequence stratigraphic framework: the example of the Lias of the Paris Basin. In: Paleogeography, Paleoclimate, ands Source Rocks[END_REF][START_REF] Disnar | Sedimentation de la matière organique dans le nord-est du Bassin de Paris : conséquences sur le dépôt des argilites carbonées du Toarcien inférieur[END_REF]. We present here a summary of our dataset (see In the Paris Basin, the lowermost organic-rich level (a few metres thick) in the Hettangian is located in the eastern part of the basin and is strongly influenced by terrestrial organic matter inputs with strong variations in TOC (from 0 up to 4 wt. %). A thicker interval (up to 60 m) is observed within the depocentre in the Upper Sinemurian, with an average TOC of 3 wt. %. The base of the Pliensbachian shows a more distinct organic-rich interval attaining a thickness of 30 m in the southern part of the basin, with a TOC of up to 4 wt. %. The richest interval is located at the base of the Toarcian, in the Schistes Carton Formation that can be observed throughout the basin, reaching a thickness of 70 m with TOC values ranging from 3 to 7 wt. %. In the South-East Basin, thermal maturation of the Lower Jurassic deposits has led to the expulsion of most of the organic matter

M A N U S C R I P T A C C E P T E D
present in the source rocks, and the initial TOC is difficult to determine (e.g. [START_REF] Wannesson | Bassin du Sud-Est. Rapport régional d'évaluation pétrolière[END_REF]. However, on the northern edge of the South-East Basin, in the Bresse/Jura area, the observations are similar to the eastern part of the Paris Basin. In the Causses area, whose connection with the South-East Basin developed during the Early Jurassic and which was less deeply buried, the TOC is very variable and, within the Toarcian, it can locally reach 12 wt. % in certain intervals several meters thick. In the English Channel Basin, the main organic-rich interval is located in the Upper Hettangian -Lower Sinemurian, and may reach a thickness of 40 m. The TOC is very variable, with an average value of 4 wt. % and pulses reaching 25 wt. %. In the Aquitaine Basin, two organic-rich levels are located at the top of the Sinemurian and at the top of the Pliensbachian,

showing very variable TOC (from 1 to 5 wt. %) over an interval of around 10 m. Another interval is located in the Toarcian in the centre of the basin, yielding a narrower range of TOC values (from 3 to 5 wt. %).

M A N U S C R I P T A C C E P T E D

Dataset and methodology

Dataset

The study area covers an area of 1,000×1,000 km extending over the whole of France (Fig. 1), and lays particular emphasis on the Lower Jurassic interval from the Lower Hettangian (201.3 Ma) to the Upper Toarcian (174.2 Ma) (Fig. 5).

The study is based on data from outcrops and boreholes in the different basins (red and green dots in Fig. 1), including structural setting, biostratigraphy, sedimentary characterization and mineralogical/geochemical analyses (29 wells including wireline logging, 12 carbologs, and more than 3000 Rock-Eval and mineralogical analyses). This dataset is built from the compilation of numerous studies cited below as well as our own field and laboratory investigations. These data are summarized in the geological setting section. They are used in the model for calibration and chronostratigraphic, lithological and geochemical control.

Stratigraphic forward modelling

We use the DionisosFlow© 3D stratigraphic forward modelling software [START_REF] Granjeon | Modelisation stratigraphique deterministe -conception et applications d'un modele diffusif 3D multilithologique[END_REF][START_REF] Granjeon | 3D stratigraphic modeling of sedimentary basins[END_REF][START_REF] Granjeon | Concepts and applications of a 3-D multiple lithology, diffusive model in stratigraphic modeling[END_REF][START_REF] Granjeon | 3D stratigraphic modeling in complex tectonics area[END_REF] to reconstruct the stratigraphic and palaeo-

Basin morphology construction

Water depth varies between 0 and 250 m over the same period. The sea covers between 60 and 80% of the studied area during the Early Jurassic. The modelling starts at 200.9 Ma (Early Hettangian) and ends at 174.2 Ma (Late Toarcian), with a time step of 0.1 Ma (Fig. 5). The horizons dated at Using a backstripping methodology applied to bathymetries and thicknesses, these data allow us to compute the available accommodation. In this way, we are able to perform simulations of sediment production and diffusion. Backstripping is generally used to estimate tectonic subsidence by accounting for and removing the effects of other causes of subsidence, such as loading due to the sedimentary column. Here, we compute the evolution of accommodation space (Eustasy + tectonic subsidence) with the equations:

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 22 (1) Thickness[t] = Depth[t-1] -Depth[t] (2) DeltaAccommodation[t] = Thickness[t] + Bathymetry[t] -Bathymetry[t-1] (3) Accommodation[t] = Accommodation [t-1] + DeltaAccommodation [t]
were the thickness is corrected of the compaction effect. Sediment was assumed to be composed of three main sediment classes: sand, shale and limestone/dolostone. Compaction of each class was defined using porosity (Φ) vs. burial depth curve (z), as described in the Figure 9 for parameters of each lithology.

Sediment production and diffusion

The sediments simulated in DionisosFlow© models are defined as a function of production and diffusion processes. Sediments are divided into in six classes for this study:

1) Fine carbonate: associated with pelagic biogenic deposits formed by primary producers in verylow-energy environments, composed of mud-supported carbonate with wackestone to mudstone textures. This sediment can be found in the low-energy basin facies; 2) Bioclastic carbonate: associated mainly with pelagic biogenic deposits but also with bioclasts, composed of mud- The diffusion coefficients, Ks and Kw (Table 2), were defined using a simple geometrical rule assuming that sediment transport is ruled by a linear slope-driven and, for clastic inputs, a nonlinear water-driven diffusion law (n = 1.5):

= × + × ×
where Qs is the volume sediment flux (km 2 .kyrs -1 ), S the slope gradient, Q w the dimensionless water flow.

Fine carbonate

Bioclastic carbonate

Shallow carbo/dolostone Shale Sand Evaporite

Marine Ks (km².kyrs -1 ) 0.075 0.1 0.05 0.5 0.05 0.1

Marine Kw (km².kyrs -1 ) 1000 200

Table 2. Diffusion coefficients for the different sediments.

Facies and thickness calibration

Trial and error simulations are used to determine sediment production and to reproduce the thicknesses and facies observed in well and field data. Simulated Facies are defined for each cell of the model according to the bathymetry and percentage abundance of each simulated sediment type

M A N U S C R I P T A C C E P T E D (Table.
3), and the simulated facies are compared against the observed facies. Figure 10 presents the final sediment production history (function of time and bathymetry) used for the model calibration.

The best fit-model is presented in the 3D stratigraphic forward modelling section.

Shoreface Inner Ramp

Mid

Organic matter modelling

The stratigraphic forward modelling approach provides most of the parameters required for organic matter simulation (e.g. water depth, basin morphology, sedimentation rate). New parameters are

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

added into our best fit-model to simulate all the processes needed for organic matter modelling [START_REF] Granjeon | Sedimentary basin development method using stratigraphic simulation coupled with an organic matter production and degradation model[END_REF][START_REF] Bruneau | 3D Numerical Modelling of Marine Organic Matter Distribution: example of the Early Jurassic sequences of the Lusitanian Basin (Portugal)[END_REF]: primary productivity, carbon flux, organic matter transport, dissolved oxygen level and burial efficiency (which corresponds to degradation within the topmost metre of burial). All these processes and their parameters are based on empirical equations or observations summarized in [START_REF] Granjeon | Sedimentary basin development method using stratigraphic simulation coupled with an organic matter production and degradation model[END_REF] and in [START_REF] Bruneau | 3D Numerical Modelling of Marine Organic Matter Distribution: example of the Early Jurassic sequences of the Lusitanian Basin (Portugal)[END_REF].

Primary productivity

We define the sea-surface primary productivity near the shore as remaining constant with time at a value of 400 gC.m -2 .yr -1 , and a seaward stepwise decrease of 100 gC.m -2 .yr -1 is imposed every 50 km with a minimum value of 200 gC.m -2 .yr -1 . This seaward decrease is representative of long time scales (e.g. [START_REF] Berger | Productivity of the Ocean: Present and Past[END_REF][START_REF] Pilson | An introduction to the chemistry of the sea[END_REF]Sathyendranath and Platt, 2013).

Carbon Flux

After primary production, the organic particles sink to the seafloor (carbon flux; e.g. Martin et al., 1987;[START_REF] Berger | Productivity of the Ocean: Present and Past[END_REF]Honjo et al., 2008). The exported production (CF) reaching the sea-floor interface is determined using the Martin equation (Martin et al., 1987) (i.e. an exponential decrease with water depth):

= 0.409 . ( 100
)
where CF is the exported organic flux (gC.m -2 .yr -1 ), PP is the primary productivity (gC.m -2 .yr -1 ), z the water depth, and n is assumed to be constant and equal to 0.86.

Organic matter transport

Once the organic remaining matter reaches the seafloor in the form of flocs, it is eventually transported. The diffusion coefficient of silt is used to mimic the transport regime of flocs at the sediment/water interface (Schieber et al., 2007). In our model, diffusion coefficients of flocs (organo-mineral aggregates) are:

Marine Ks = 0.02 km².kyrs -1

Marine Kw = 500 km².kyrs -1

Oxygen level

In DionisosFlow© Organic-Rich Sediment, a mass balance equation between oxygen renewal and oxygen consumption allows us to determine the relative oxygen depletion based on the studies of Mann and Zweigel (2008). This depletion is currently estimated without considering any water current so it is merely 1D (vertical):

( ) = ( ) × !" -$ × ( )
where Ox z is the oxygen level at sea-floor (expressed in nondimensional units, with values varying from 0 for fully anoxic environment to 1 for fully oxic environment), A z (z) is the mixing rate coefficient, Ox REF the reference oxygen level and D z (z) the oxygen consumption

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

29 linked to organic matter degradation by organism. In our model, the mixing rate is constant in time and space and the dissolved oxygen level is therefore dependant only on bathymetry and primary productivity.

Burial efficiency

The amount of organic matter preserved after the first few metres of burial (i.e. during early diagenesis) is called burial efficiency; this parameter is mainly controlled by the sedimentation rate [START_REF] Betts | The oxygen content of ocean bottom waters, the burial efficiency of organic carbon, and the regulation of atmospheric oxygen[END_REF] and local oxygen level/redox conditions (Tyson, 1995) according to the equation:

%&' = -( %&' → %&' = -( %&'
where Corg is the quantity of available organic Carbone, k the degradation coefficient (itself function of the burial depth into the sediment, k = k 0 e -z/z0 ) and S the sedimentation rate.

Modelling of these processes allows us to compute a theoretical initial TOC after the first few metres of burial of the marine organic matter fraction.

Total organic carbon

The Total Organic Carbon (TOC) is eventually determined with the amount of preserved organic matter diluted depending on the sedimentation rate. It corresponds to an initial TOC after the first few metres of burial of the marine organic matter fraction. [START_REF] Bruneau | 3D Numerical Modelling of Marine Organic Matter Distribution: example of the Early Jurassic sequences of the Lusitanian Basin (Portugal)[END_REF] provided detailed sensitivity analyses of this initial TOC distribution with several sets of parameters for the processes

M A N U S C R I P T A C C E P T E D described before.
In the rest of the study, "simulated TOC" refer to TOC obtained by modelling. This TOC is calculated as an average value by cell (10×10 km) and by time step (100 kyrs). This resolution imply some caveats when comparing simulated TOC with TOC from well and field data.

However, the equations described in this section, are proper to large time and space scales, and the simulation of the processes governing organic matter distribution is relevant.

Results

Best fit-model

Following the above detailed methodology, the 3D model is computed for the basins studied here based on subsurface data, field observations and measurements (Fig. 11). The volume and the distribution of sediments in the basin are defined using trial and error simulations to reproduce observed thicknesses and facies. The simulated facies, which are used to compare observed facies, are defined for each cell of the model according to the percentage of each simulated sediment. The average sediment thickness error between simulation and real thickness is usually less than 10%

(around 50 m on the total sediment column simulated in the Paris Basin). Both vertical and lateral simulated facies variations match with interpreted data from field studies and the literature. Only few cells show discrepancies linked to the difficulty in precisely reproducing the lateral facies changes.

In such cases, cells with appropriate facies can be found within a maximum radius of 2 cells (10 km) from the calibration point.

In detail, some local discrepancies remain concerning the modelled thickness and facies compared with the present-day situation, since changes of thickness due to post-depositional deformation are not taken into account (except for effects due to compaction). In the eastern part of the South-East

Basin, the discrepancy reaches 10 % due to deformation of the Cretaceous-Cenozoic Alpine orogenic phase. Other minor differences in thickness may be due to the difficulty of assigning production values for each sedimentary facies in all the basins. In spite of this complexity, the facies and thicknesses computed from the 3D model are in agreement with all the major geological features observed in the basins.

Fig. 11. 3D block of the best-fit model of the depositional environments for the end of the

Sinemurian, calibrated using facies distribution and thickness.

Global organic matter distribution and content

Despite constant boundary conditions in terms of primary productivity and dissolved oxygen variation, the simulated TOC shows a broad range of values (from 0 to 20 wt. %). Areas with high TOC content (source rock) can be simulated using internal factors alone (e.g. basin morphology, bathymetry and sedimentation rate) without involving external factors (e.g. oceanic anoxic events,

M A N U S C R I P T A C C E P T E D

upwelling). TOC maps show organic-rich areas that are relatively limited in space even during the Toarcian. These areas are most often located in the depocentres and less frequently in the proximal areas of the basins.

Our results cannot be compared against data from areas that are nowadays eroded or highly deformed (e.g. current rim of the Paris Basin, eastern part of the South-East Basin, Bresse-Jura and Poitou threshold). However, despite the constant primary production and the absence of any variations of dissolved oxygen level linked to other factors such as primary productivity or bathymetry, the distribution of organic-rich sediments given by our model fits with available data. To highlight the reliability of the model, we present close-up views of some cross-sections and wells (Figs. 12, 13, 14 and 15).

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Nonetheless, we can note a misfit between the simulated TOC and the well data. The first difference stems from the burial effect, which is particularly evident in the Aquitaine and South-East basins where organic-rich layers can be buried under more than 10 km of Mezosoic and Cenozoic deposits (Fig. 13 and15). In these areas, our model shows thick organic-rich layers with initial TOC reaching 20 wt. % as against 5 wt. % in wells. The burial and maturation history leads to expulsion of a large proportion of the organic matter, which hinders any comparison between modelling and well data. In the rest of the model, the overall trends of organic matter content in the sediments are well reproduced (Fig. 12 and14). The 0.1-Myr time steps yield an average TOC over intervals ranging from a few cm to more than 10 m depending on the sedimentation rate. Hence, higher frequency variations are not taken into account in our modelled TOC values, and comparisons need to be performed using the trends and average measured values of each interval (the TOC curve is derived from logging data, with a value every 0.1 to 5 m, whereas the thickness of the model grid cells, based on time steps, varies from 0 to 10 m). The duration of TOC pulses and the importance of the time scale in comparing the simulated and measured TOC remains an open question. Furthermore, to achieve an accurate fit with the TOC data, we need to consider the time-evolution of primary production and dissolved oxygen levels (e.g. Toarcian anoxic event).

M A N U S C R I P T A C C E P T E D

ACCEPTED MANUSCRIPT 39

Discussion

Organic-rich sediment depocentres

The main organic-rich deposits are located in depocentres of the basins. This is particularly well seen in the Paris and South-East basins (Fig. 12, 13, 16 and 17-f.) in the models. The initial TOC values exceeds 4% and can reach 6 wt. % over intervals up to tens of metres thick. These areas are mainly composed of shaly-marly to shaly facies with high sedimentation rates (around 50 to 70 m.Myr -1 ), which are always located several hundreds of km from the coastline (from 100 to 300 km).

Therefore, in our modelling, the primary productivity at the sea surface directly above these deposits is quite low and never higher than 200 gC.m -2 .yr -1 , since it decreases as a function of distance from the coast . The depocentres also correspond to the deepest part of the basins, with minimum water depths varying from 150 m for the Paris Basin (Fig. 17 In view of these results, it seems that the accumulation of thick organic-rich layers is more closely related to the transport of shales in the deepest part of the basin. [START_REF] Passey | From Oil-Prone Source Rock to Gas-Producing Shale Reservoir -Geologic and Petrophysical Characterization of Unconventional Shale Gas Reservoirs[END_REF] has described this type of distribution as a distal parasequence on the platform ramp. When the carbonate ramp attains a greater water depth, this leads to preferential accumulation of shales and organic matter.

This setting is therefore more closely related to water depth than the development of a chronostratigraphic sequence. In the model results, we observe organic-rich depocentres in both regressive and transgressive phases. However, these phases have an effect on the distribution of the organic-rich deposits: their accumulation is always associated with the onset transgression (when bathymetry is sufficient to ensure gradation of carbonate ramp into basinal facies) and then spreads during the rest of the transgression as described by [START_REF] Wignall | Model for transgressive black shales?[END_REF] and [START_REF] Röhl | Lower Toarcian (upper liassic) black shales of the central European epicontinental basin: a sequence stratigraphic case study from the SW German Posidonia shale[END_REF]. But, as the carbonate ramps invade the emerged land, the shaly and organic-rich facies also shift towards the proximal domain. According to two main scenarios, sedimentation is a function of the slope gradient and the size of the basin (which themselves depend on sediment accumulation and subsidence): (i) if the slope gradient is steep enough, the shales and organic matter continue to be accumulated on top of the previous depocentre and can even spread toward the proximal part (e.g.

South-East basin) (Fig. 14) as described by [START_REF] Wignall | Model for transgressive black shales?[END_REF], (ii) if the slope is gentle in the deepest part, the organic matter is accumulated directly at the extremity of the carbonate ramp and a smaller proportion is transported basinwards where it is further diluted by shales (e.g. Paris Basin) (Fig. 16 and 17). This results in an organic-rich "ring" surrounding moderately rich deposits in the core of the basin (Fig. 17). Clearly, the relatively shallow epicontinental basin setting allows the development of organic-rich deposits, which might not be the case in other sedimentary basins. This phenomena also reflects the impact of sedimentation rate: in the case of the Paris Basin, under fairly well oxygenated conditions, sedimentation rates higher than 40 m.My -1 increase the preservation of organic matter due to fast burial and isolation from oxygen. However, if the sedimentation rate exceeds 60 m. My -1 , even if the preservation is increased, dilution will reduce the TOC of the sediment (Ibach, 1982;[START_REF] Zonneveld | Selective preservation of organic matter in marine environments; processes and impact on the sedimentary record[END_REF]. This optimal window allows a good preservation of organic matter in areas where dissolved oxygen levels remain high. In dysoxic to anoxic environments, the preservation is already ensured and lower sedimentation rates lead to lower dilution at higher TOC (Fig. 17).

Furthermore, a migration of organic-rich sediment can be observed across the basin. This is clearly seen in the Paris Basin (Fig. 16), where the Toarcian high-TOC deposits migrate in the same direction as the depocentre, associated with good preservation conditions. In this way, while any given organic-rich horizon can be followed throughout the basin, the deposits are diachronous as already assumed in the literature (e.g. Jenkyns et al., 2002;Mailliot et al., 2006;Hesselbo et al., 2007).

Proximal areas

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Other organic-rich deposits can sometimes be observed in proximal areas at distances of 10 to 30 km from the coastline. In such settings, the initial TOC values are around 3 to 6 wt. % and the organicrich layers are never more than a few meters thick, but can be followed for several km parallel to the shore. These deposits formed in shallow-water environments are composed of carbonates and/or clastics, with sedimentation rates of less 20 m.My -1 . The water depth can vary from a few to several tens of meters. Because of the proximity to the coastline, the sea-surface primary productivity in these areas is high (from 300 to 400 gC.m -2 .yr -1 ). Even if this high primary productivity leads to high consumption of oxygen, the shallow water depths allow an efficient renewal, and these environments are fully oxygenated. These depositional environments develop during transgressive phases, when the accumulation of thin organic-rich layers is favoured by the topography of the carbonate ramp.

Since the proximal deposits of the ramp are relatively flat, restricted depressions can form because of differential subsidence or clastic supply. In this case, even if the level of dissolved oxygen theoretically leads to the total degradation of organic matter, the accumulation of sediments is concentrated in local depressions, especially since organic matter is transported as silt-sized particles that are abundant in high primary productivity areas. The sedimentation rate can locally reach the optimal windows described above and thus enhance the preservation of organic matter.

The same process can be observed during regressive phases on confined sills (e.g. modelled Causses and Poitou sill). Here, the proximity of emerged land leads to high primary productivity (300 to 400 gC.m -2 .yr -1 ) over very large areas. TOC values higher than 5 wt. % can be found because of the high organic matter content of sediments deposited on this threshold. Organic-rich layers can be deposited

M A N U S C R I P T A C C E P T E D
even if the environment is fully oxygenated and with sedimentation rates below the optimal windows. However, more open thresholds (e.g. modelled Bresse-Jura area) do not display the same feature, even in cases where primary productivity is high and sedimentation rates are more favourable for organic matter preservation. These observations imply that, under the right conditions, an extensive sea area (at least 5,000 km²) with moderate primary productivity can lead to high TOC rather than localized areas with higher primary productivity.

Conclusion

This study aims to analyse the processes controlling organic matter distribution and heterogeneity through a numerical quantitative approach applied to the Lower Jurassic succession of north-western Europe. The 3D modelling approach allows us to incorporate a large amount of data and test multiple coupled processes. A high primary productivity as well as factors favouring organic matter preservation could theoretically create conditions for the development of high organic carbon contents in sediments. Modelling results and associated interpretations allow us to draw several conclusions concerning the evolution of organic matter and the key factors controlling its distribution and heterogeneity:

(1) Accumulation of thick organic-rich layers in depocentres are correlated with shaly-marly to shaly facies with high sedimentation rates (around 50 to 70 m.Myr -1 ) as well as areas with decreased dissolved oxygen levels due to isolation and/or consumption. These settings are developed during transgressive phases, and also during regressive phases if significant subsidence and sediment accumulation continues in the deepest part of the basins.

M A N U S C R I P T A C C E P T E D

(2) During transgression, organic matter are accumulated in depocentres on the top of previous depocentres but can also spread toward the more proximal parts of the basin.

(3) If the slope is gentle, an organic-rich "ring" is accumulated at the break-in-slope of the carbonate ramp, and a smaller proportion is transported basinward where it is further diluted by shales.

(4) Organic-rich carbonate and/or coarse clastic layers a few meters thick can be followed over several kilometres parallel to the shore; these layers can be observed during transgressive phases or during regressive phases on confined threshold. They are related to restricted depressions linked to differential subsidence and/or clastic supply in the proximal flat part of the ramp. In this setting, the water depth (from a few to several tens of meters) and sedimentation rates (less than 20 m.My -1 ) hinder the preservation of organic matter but, because of the proximity of the coast line, the primary productivity at the sea surface is high.

Several different favourable conditions are highlighted for the accumulation organic-rich deposits.

The main factors leading to thick organic-rich deposits are sedimentation rate (between 40 and 60 m.

My -1 for a good preservation) and basin geometry (depocentre position). The main factors leading to restricted and thinner organic-rich deposits in proximal settings or on thresholds are primary productivity and the existence of local depressions. These depositional environments for organic-rich sediments are compatible with many of the conceptual models currently proposed in the literature (e.g. [START_REF] Wignall | Model for transgressive black shales?[END_REF][START_REF] Röhl | Lower Toarcian (upper liassic) black shales of the central European epicontinental basin: a sequence stratigraphic case study from the SW German Posidonia shale[END_REF][START_REF] Passey | From Oil-Prone Source Rock to Gas-Producing Shale Reservoir -Geologic and Petrophysical Characterization of Unconventional Shale Gas Reservoirs[END_REF].

This study highlights the need for a relative quantification of the palaeo-environmental conditions and processes that lead to the widespread accumulation of organic matter.

  settings changes which have favoured the formation of organic-rich source rocks.

Fig. 1

 1 Fig. 1. a. Paleogeography of the western margin of Tethys Ocean during the Early Jurassic

Fig. 2 .

 2 Fig. 2. Thickness variations of the three main Lower Jurassic source rocks of the Paris Basin.

Fig. 4 .

 4 Fig. 4. Lower Jurassic sequence chronostratigraphy and eustatic variations (modified after

Fig. 5 .

 5 Fig. 5. Chronostratigraphic chart for the Lower Jurassic series of the north-western European

  Fig.6. Carbonate ramp settings (modified after[START_REF] Burchette | Carbonate ramp depositional systems[END_REF]. Inner Ramp

  the studied basins during the Early Jurassic. DionisosFlow© simulates the evolution of the basin geometry through time, using long-term and large-scale sediment transport equations. The results of this simulation allow us to model changes in facies distribution and sedimentation rate at the basin scale.

200. 9

 9 Fig. 7. Paleo-bathymetric maps of the Lower Jurassic.

Fig. 9 .

 9 Fig. 9. Porosity function of depth for the different head member sediments of the Lower Jurassic

  supported carbonates with wackestone to packstone textures. It can be found in all the facies of the ramp but mainly in the ramp-slope facies; 3) Shallow carbonate/dolostone: associated mainly with dolomite but also oolites, bioclasts, bivalves and other producer organisms in high-energy environments, composed of grain-supported carbonates with packstone to grainstone textures mainly corresponding to the inner to mid-ramp facies; 4) Shale: associated with fine clastic components of sediments associated with low-energy conditions of the basin facies; 5) Sand: associated with medium to coarse clastic component; 6) Evaporite: associated with evaporitic deposits located in confined environments.These sediments have different transport properties based on the coupling of slow downslope linear diffusion and fast fluvial nonlinear water-driven diffusion, which allows the simulation of sediment movement at the kilometre scale and over long time spans (e.g.[START_REF] Rabineau | Sedimentary sequences in the Gulf of Lion: a record of 100,000 years climatic cycles[END_REF][START_REF] Alzaga-Ruiz | Gravitational collapse and Neogene sediment transfer across the western margin of the Gulf of Mexico: insights from numerical models[END_REF]. The sediments are assumed to be composed of a given number of grain-size fractions. The M A N U S C R I P T A C C E P T E D mass conservation equation is applied to each grain-size fraction. Diffusive coefficients are discriminated according to grain size (e.g. fine-grained fractions such as shale are transported more easily because of the greater degree of sediment suspension).

Fig. 10 .

 10 Fig. 10. Sediment production evolution function of time and bathymetry.

  Fig. 12. a. Cross-section of the Paris Basin in current structural position and at the end of the

  -b.) to 400 m for the South-East Basin. Under these conditions, the dissolved oxygen level can become drastically lowered (e.g. South-East Basin) or simply depleted (e.g. Paris Basin, Fig.17-e.). This oxygen depletion is due solely to bathymetry because the high primary productivity zones are too far away to have a noticeable effect on oxygen consumption. However, the highest TOC values are not found in the area with optimum sedimentation rate (between 40 and 60 m.My -1 ), but in areas with lower sedimentation rate (>40 m My -1 , Fig.17-c.) In such environments, the preservation of organic matter is already -depleted conditions, thus the lower sedimentation rate results in less dilution of the TOC signal.In terms of sequence stratigraphy, these conditions are developed during transgressive phases, and the maximum spatial extent of these deposits is reached during or at the end of the maximum flooding tract. These criteria are also reached during regressive phases provided there is still a significant subsidence and sediment accumulation in the deepest part of the basins.

Fig. 16 .Fig. 17 .

 1617 Fig. 16. Thickness and TOC distribution of Sinemurian and Toarcian source rocks in the Paris

  For the time interval studied here, covering 25 Myr during the Early Jurassic without any major tectonic event and without assuming any variation in primary productivity or open ocean oxygen levels, the model predicts several source rocks and TOC ranging from 0 to 20 wt. %.

Table 1 )

 1

	Basin, sill or	Source Rock Age	Av TOC Min/max	Thickness

Table 1 .

 1 Source rock characteristics of the different basins. The TOC values show strong variations

from one basin to another and also over time in the same basin. The Liassic succession of the Poitou sill is nowadays eroded. * In the South-East Basin, the initial TOC values are not comparable to current TOC contents because of the strong impact of burial and maturation in this area.

ACKNOWLEDGMENTS

DionisosFlow© is developed in the framework of the IFP Energies Nouvelles project DORS ("Dionisos Organic-Rich Sediment"). We are grateful to Brigitte Doligez, Bertrand Coureaud, Jean-François Deronzier, Jean-Pierre Barde and Laurent Riquier for constructive discussions. We also thank Florence Savignac and Marie Salpin for laboratory technical assistance and BeicipFranlab for support with DionisosFlow©. We wish to thank Fadi H. Nader and two anonymous reviewers for their constructive remarks. Carpenter post-edited the English style and grammar.

Clearly, some of our hypotheses could be challenged and the approach used to study certain processes such as organic-matter transport merits further investigation and improvement. Variations of primary productivity have an impact in the proximal parts of basins and in threshold areas, while climatic change and other external causes may also have an influence.

CONFLICT OF INTEREST

No conflict of interest declared. Mann, U., Zweigel, J., 2008. Modeling source-rock distribution and quality variations: the organic facies modelling approach. Spec. Publ. Int. Assoc. Sedimentol., 40, 139-274. Martin, J.H., Knauer, G.A., Karl, D.M., Broenkow, W.W., 1987. VERTEX: carbon


Vahid Roostapour
email: v.roostapour@iasbs.ac.ir

Iman Kiarazm
email: i.kiarazm@iasbs.ac.ir

Mansoor Davoodi

Deterministic Algorithm for 1-Median 1-Center Two-Objective Optimization Problem

Keywords: computational geometry, Pareto optimal solutions, 1-center, 1-median, multi-objective optimization

k-median and k-center are two well-known problems in facility location which play an important role in operation research, management science, clustering and computational geometry. To the best of our knowledge, although these problems have lots of applications, they have never been studied together simultaneously as a multi objective optimization problem. Multi-objective optimization has been applied in many fields of science where optimal decisions need to be taken in the presence of trade-offs between two or more conflicting objectives. In this paper we consider 1-median and 1-center two-objective optimization problem. We prove that Ω(n log n) is a lower bound for proposed problem in one and two dimensions in Manhattan metric. Also, by using the properties of farthest point Voronoi diagram, we present a deterministic algorithm which output the Pareto Front and Pareto Optimal Solutions in O(n log n) time.

Introduction

When evaluating different solutions from a design space, it is often the case that more than one criterion comes into play. For example, when choosing a route to drive from one point to another, we may care about the time it takes, the distance traveled and the complexity of the route (e.g. number of turns). When designing a (wired or wireless) network, we may consider its cost, capacity and coverage. Such problems are known as Multi-Objective Optimization Problems (MOOP). Multi-objective optimization can be described in mathematical terms as follows:

S = {x ∈ R d : h(x) = 0, g(x) ≥ 0} min [f 1 (x), f 2 (x), . . . , f N (x)] x ∈ S,
where N > 1, f i is a scalar function for 1 ≤ i ≤ N and S is the set of constraints.

The space in which the objective vector belongs is called objective space. The scalar concept of optimality does not apply directly in the multi-objective setting.

Here the notion of Pareto optimality and dominance has to be introduced. In a multi-objective minimization problem, a solution s 1 ∈ S dominates a solution s 2 ∈ S, denoted by s 1 ≺ s 2 , if f i (s 1 ) ≤ f i (s 2 ) for all i ∈ {1, . . . , N }, with at least one strict inequality. A point s * is said to be a Pareto optimum or a Pareto optimal solution for the multi-objective problem if and only if there is no s ∈ S such that s ≺ s * . The image of such an efficient set, ie., the image of all the efficient solutions in the objective space are called Pareto optimal front or Pareto curve.

One of the common approaches for such problems is evolutionary algorithms [START_REF] Coello | Evolutionary algorithms for solving multi-objective problems[END_REF]. These algorithms are iterative and converge to Pareto front. However they need more time as the complexity of the Pareto front increases. Moreover, all of these approaches have major problems with local optimums. On the other hand there are some classical approaches like weighted sum and ϵ-constraint which can apply on MOOPs. Although these approaches guarantee finding solutions on the entire Pareto optimal set for problems having a convex Pareto front, they are largely depend on chosen weight and ϵ vectors respectively. Moreover, these approaches require some information from user about the solution space. Furthermore, in most nonlinear MOOPs, a uniformly distributed set of weight vectors wont necessarily find a uniformly distributed set of Pareto optimal solutions. Also there may exist multiple minimum solutions for a specific weight vector [START_REF] Deb | Multi-objective optimization using evolutionary algorithms[END_REF]. However we find the Pareto front of a MOOP with deterministic algorithm. Here we consider two famous propounded facility location problems [START_REF] Tansel | State of the artlocation on networks: a survey. part i: the p-center and p-median problems[END_REF].

k-median: In this problem the goal is to minimize summation of distances between each demand point and its nearest center. Charikar et. al proposed the first constant time approximation algorithm which its outputs is 6 2 3 times the optimal [START_REF] Charikar | A constant-factor approximation algorithm for the k-median problem[END_REF]. This improved upon the best previously known result of O(log p log log p), which was obtained by refining and derandomizing a randomized O(log n log log n)approximation algorithm of Bartal [START_REF]Probabilistic approximation of metric spaces and its algorithmic applications[END_REF]. The currently best known approximation ratio is 3+ϵ achieved by a local search heuristic of Arya et. al [START_REF] Arya | Local search heuristics for k-median and facility location problems[END_REF]. Moreover, Jain et. al proved that the k-median problem cannot be approximated within a factor strictly less than 1 + 2/e, unless NP ⊆ DTIME[n O(log log n) ] [START_REF] Jain | A new greedy approach for facility location problems[END_REF]. This was an improvement over a lower bound of 1+1/e [START_REF] Shmoys | Approximation algorithms for facility location problems[END_REF]. Using sampling technique Meyerson, et. al presented an algorithm with running time O(p( p 2 ϵ log p) 2 log( p ϵ log p)). This was the first k-median algorithm with fully polynomial running time that was independent of n, the size of the data set. It presented a solution that is, with high probability, an O(1)-approximation, if each cluster in some optimal solution has Ω( n•ϵ p ) points [START_REF] Meyerson | A k-median algorithm with running time independent of data size[END_REF]. Har-Peled and Kushal presented a (p, ϵ)-coreset of size O(p 2 /ϵ d ) for k-median clustering of n points in R d , which its size was independent of n [START_REF] Har-Peled | Smaller coresets for k-median and k-means clustering[END_REF]. Also, Har-Peled and Mazumdar showed that there exist small coresets of size O(pϵ -d log n) for the problems of computing k-median clustering for points in low dimension with (1 + ϵ)-approximation. Their algorithm has linear running time for a fixed p and ϵ [START_REF] Har-Peled | Coresets for k-means and k-median clustering and their applications[END_REF]. Moreover, using random sampling for k-median problem Badoiu et. al proposed a (1 + ϵ)-approximation algorithm with 2 (p/ϵ) O(1) d O (1) n log O(p) n expected time [START_REF] Bādoiu | Approximate clustering via core-sets[END_REF].

k-center:

In this problem the goal is to minimize the maximum distance between each demand point from its nearest center. Megiddo and Supowit proved that k-center and k-median are NP-hard even to approximate the k-center problems sufficiently closely [START_REF] Megiddo | On the complexity of some common geometric location problems[END_REF]. Hochbaum and Shmoys proposed the first constant factor approximation algorithm which its output is 2 times the optimal. It is the best possible algorithm unless P ̸ = NP [START_REF] Hochbaum | A best possible approximation algorithm for the k-center problem[END_REF]. It is shown that there is an algorithm with O(d O(d) n) time for 1-center problem [START_REF] Chazelle | On linear-time deterministic algorithms for optimization problems in fixed dimension[END_REF]. In the high dimension, Badoiu and Clarkson presented a (1 + ϵ)-approximation algorithm which find a solution in ⌈2/ϵ⌉ passes using O(nd/ϵ + (1/ϵ) 5 ) total time and O(d/ϵ) space [START_REF] Badoiu | Smaller core-sets for balls[END_REF]. Also, for problem of 1-center with outliers, Zarrabi-Zadeh and Mukhopadhyay proposed a 2-approximation one pass streaming algorithm in high dimension which for z, as the number of outliers, needs O(zd 2 ) space [START_REF] Zarrabi-Zadeh | Streaming 1-center with outliers in high dimensions[END_REF]. Moreover, Zarrabi-Zadeh and Chan presented an streaming one pass 3/2-approximation algorithm for 1-center [START_REF] Zarrabi-Zadeh | A simple streaming algorithm for minimum enclosing balls[END_REF]. Badoiu et. al for 1-center problem, extracted a coreset of size O(1/ϵ 2 ) which its solution is (1+ϵ)-approximation set of points in R d [START_REF] Bādoiu | Approximate clustering via core-sets[END_REF]. Also, for k-center they presented a 2 O((p log p)/ϵ 2 ) .dn time algorithm with (1 + ϵ)-approximation solution using previous result.

1-median and 1-center are practical problems which have not been considered as a two-objective optimization problem yet. Imagine mayor of a small city wants to build a fire station in a way that minimizes the distance between farthest building to the station, also since the number of fire engines is limited and each fire engine must return to the station after a service, it has to minimizes the total distance of station from all other buildings. As an another example, consider power distribution network. Due to the dependency of energy leakage to wire length, minimizing of the longest wire in the network would be regarded as an essential factor. Also, any decrement in total wire length of network considered as a second objective. The first objective is 1-median, M (u), the summation of distances of demand points from center u and the second objective is 1-center, C(u), the farthest input point from center p. It can be described in mathematical terms as follow: We study this problem in one and two dimensions in Manhattan metric. We assume no input points have the same x or y coordinate. This is a convex combinatorial multi-objective optimization problem which has been studied with a different approach called ϵ-Pareto. In [START_REF] Papadimitriou | On the approximability of trade-offs and optimal access of web sources[END_REF] it is shown that this approximate Pareto curve can be constructed in time polynomial in the size of the instance and 1/ϵ, but here we propose a deterministic algorithm for computing the exact Pareto curve because of specifying the problem.

This paper starts with considering 1-median and 1-center as two-objective of MOOP in one dimension. We will find the optimal of objectives and in terms of placement of optimums we will also find the Pareto set in time O(n) (Lemma 1).

(a) n is odd.

(b) n is even.

Fig. 1: 1-Median optimal.

We continue with a proof for convexity of Pareto set. At the end of second section we give an algorithm to compute Pareto optimal front of 1-median 1-center two-objective optimization problem and prove the optimality of the algorithm.

In section three the same problem considered in two dimensional space. First we find optimums of 1-median objective. After that by using the properties of farthest point Voronoi diagram we determine the optimum of 1-center. Finally after limiting the solution space to regions which Pareto set lies on, we specifically present Pareto solutions. Convexity of Pareto front is proven in Theorem 2.

One dimensional

Let P = {x 1 , x 2 , . . . x n } be a set of input points in one dimension, the goal is to minimize

M (x) = ∑ n i=1 |x -x i | and C(x) = max 1≤i≤n |x -x i |.
According to the properties of the absolute value function and some simple calculations, it is easy to see that M (x) is a continuous piecewise linear function which its minimum depends on n. The minimum can either be one point or an interval which we denote by M opt in the rest of the paper. Also without loss of generality we assume that input points are sorted increasingly. In one dimensional space,

M opt = [m i , m j ] ⊂ R for 1 ≤ i, j ≤ n such that m i = x i , m j = x j .
For odd n we have j = i and for even n, j = i + 1. Moreover, the function is strictly decreasing before its minimum and is strictly increasing after it (Figure 1). For

C(x) suppose c opt ∈ R denote the point which C(c opt ) is minimum. Obviously c opt = (x 1 + x n )/2. Similarly to M (x), C(x)
is strictly decreasing before optimal point and strictly increasing after that. Lemma 1. Pareto optimal set in one dimensional 1-median 1-center two-objective optimization problem is the smallest interval consisting of a solution with 1center optimal and a solution with 1-median optimal.

Proof. Suppose that n is even (the proof is similar for odd n). As shown in Figure 2, there are three different cases: First consider the case that c opt and M opt have an intersection (Figure 2a). In this case the intersection point is the only member of Pareto optimal solutions. Because not only it is optimal in both objectives, but also it is the only point where C(x) is optimal. So it dominates all the other solutions and no solution dominates it.

As shown in Figure 2b there are three regions in the second case. In region C both functions are strictly increasing. Therefore, c opt has the best value in both objectives. It dominates all solutions of this region. In A, C(x) is strictly decreasing, thus C(m j ) is strictly smaller than 1-center objective of all the other solutions. Moreover, M (m j ) is smaller than or equal with 1-median objective of the other solutions. Hence m j dominates all solutions of A. Finally we claim that B is Pareto set. By contradiction, suppose it is not true, then there must be a point p which dominates q ∈ B. It has to be on the left side or right side of q. Let p be on the right side, we know that M (x) is strictly increasing in this side. hence M (q) < M (p) and it contradicts with dominance of p. Similarly there is a contradiction if p lies on the left side of q, because C(x) is strictly decreasing in this side, ie. C(q) < C(p). This implies that all the solutions that lie on B are Pareto set.

The proof is similar for the third case which c opt is on the left side of M opt (Figure 2c). ⊓ ⊔ Lemma 2. Pareto optimal front of one dimensional 1-median 1-center twoobjective optimization problem forms a continuous, convex and piecewise linear function.

Proof. If there is an intersection between c opt and M opt the lemma is held. Now suppose there is no such intersection and consider c opt is on the right side of M opt (resp. on the left side of M opt ). From lemma 1 for Pareto solutions we have P s = [m j , c opt ] (resp. P s = [c opt , m i ]). Since C(x) derivation is constant and M (x) is piecewise linear in P s , the diagram of M (x)-C(x) is piecewise linear and break points are

( C(x i ), M (x i ) ) such that m j ≤ x i ≤ c opt (resp. c opt ≤ x i ≤ m i ).
The absolute value of slope of M (x) increases on each linear piece in P s . Thus Pareto optimal front is convex (Figure 3). Also we can conclude that piecewise linear Pareto front is one-to-one and invertible corresponding to Pareto solutions. ⊓ ⊔ Lemma 3. Computing Pareto front of one dimensional 1-median 1-center twoobjective optimization problem requires Ω(n log n) time.

Proof. The proof is based on reduction from sorting problem. By contradiction assume there is an algorithm which return set

O = {( C(α 1 ), M (α 1 ) ) • • • , ( C(α m ), M (α m ) )} -lexicographical
< • • • < b n+1 < a l and t = 2 • a h -b 1 + 1, then B = A ∪ {b 1 , • • • , b n+1 , t} is defined in O(n).
For the set B as input points of one dimensional 1-median 1-center two-objective optimization, 1-median optimal interval is [b n+1 , a l ] and 1-center optimal point is between a h and t. Using lemma 2 we conclude that m = n + 1 and

α 1 = a l < • • • < α m-1 = a h < α m = (b1+t) 2 .
Therefore, we can sort input points by given algorithm which implies that no algorithms with o(n log n) running time can compute Pareto front of one dimensional 1-median 1-center two-objective optimization problem.

⊓ ⊔ Note 1. If the algorithm output the Pareto optimal front as

O = {( C(α 1 ), M (α 1 ) ) - ( C(α 2 ), M (α 2 ) ) , • • • , ( C(α 2m-1 ), M (α 2m-1 ) ) - ( C(α 2m ), M (α 2m ) )}
, start points and end points of m segments, since the slope of each segment is an integer of O(n), the segments can be sorted in O(n). Therefore, we can have sorted break points of Pareto front function and the above proof holds. Proof. C(x) can be computed easily in constant time and M (x) can be computed in O(log n) using binary search, we obtain that line [START_REF] Megiddo | On the complexity of some common geometric location problems[END_REF] Due to space limitation, Algorithm 1 is just for the case that c opt is on the right side of M opt . The case that c opt is on the left side is similar. If there is an intersection, solution is obviously the intersection point.

Two dimensional

In this section we consider the problem in R 2 . The aim is to find the Pareto front and Pareto solutions in terms of M opt and C opt .

1-Median objective

For each point p ∈ R 2 we have:

M (u) = n ∑ i=1 ∥u -p i ∥ 1 = n ∑ i=1 |u x -p ix | + n ∑ i=1 |u y -p iy |
(1)

-2x + 2y + c n 2 , n 2 +2 2y + c n 2 +1, n 2 +2 2x + 2y + c n 2 +2, n 2 +2 -2x + c n 2 , n 2 +1 c n 2 +1, n 2 +1 2x + c n 2 +2, n 2 +1 -2x -2y + c n 2 , n 2 -2y + c n 2 +1, n 2 2x -2y + c n 2 +2, n 2
(a) Number of equations is even, 1-median optimal is a rectangular (blue) region.

-x -y + c n+1 2 , n+1 2
x -y + c n+3 2 , n+1 2
-x + y + c n+1 2 , n+3 2 x + y + c n+1 2 +1, n+3 2
(b) Number of is odd, 1-median is a (blue) point. We can observe that we need O(n) time to deterministically minimize equation 1. Moreover, because of the assumption that no points have same coordinate the optimal of M opt may be just a point or area of a rectangle.

In the rest of this paper we assume that n is even (all proofs and discussions are similar when n is odd.). Consider lines y = p ix and x = p iy such that 1 ≤ i ≤ n which partition the xy-plane into (n + 1) 2 cells where boundary cells are unbounded. The equation of M (p) for points in each cell is the same because of the absolute value function. Furthermore, for points in a column (resp. row) equation of

∑ n i=1 |x -p ix | (resp.
∑ n i=1 |y -p iy |) do not change but for transformation to upper (resp. right) cell coefficient of y (resp. x) increases by 2 (Figure 4).

1-Center objective

Let FVD be the farthest point Voronoi diagram of input points in Manhattan metric, also let R F VD (p) denote the region of FVD which consist of p and

, C(p) is ∥p -S FVD (R F VD (p))∥ 1 .
Besides the FVD partition the plane into at least two and at most four regions (Figure 5a).

According to the structure of FVD, it is impossible for regions A and to have a common site. However, either B (resp. D) can merge with A (resp. C) or B (resp. D) can merge with C (resp. A), ie. B and D cannot merge with a common region simultaneously.

Proposition 1. Site of region C is in A \ E.

Otherwise distances between points on segment ab and S FVD (C) are not equal and ab is not an edge of FVD. (Figure 5b)

From proposition 1 it can be concluded that C(p) for p ∈ C is equal to distance of p from segment ab add up to distance between segment ab and S F VD (C). 6a distances

Proposition 2. As shown in Figure

-m ′ 2 ∥ 1 = ∥m 2 -a∥ 1 .
The following two propositions determine the equation of C(p) in the plane and proof that it depends on which region of FVD includes p.

Proposition 3. For point p ∈ C (resp. p ∈ A), C(p) = k opt + c -p x -p y (resp. C(p) = k opt -c + p x + p y ) where c is y-intercept of ℓ 1 (Figure 6b).
Proof. Suppose equation of line ℓ 1 is y = -x + c and distance between site of C and segment ab is k opt , then projection of point p = (x, y)

on ℓ 1 is p ′ = ( c-y+x 2 , c+y-x 2
). Using proposition 1 and proposition 2 can obtain that: 7a since site of D is in hatched region or on its border, for point q ∈ D we have C(q) as the distance of point a from S F VD (D) add up to distance between point a and point q. Also since point a is an FVD vertex, we know that distance of point a from S F VD (D) is equal to its distance from S F VD (C) and as equal to k opt , hence:

C(p) = k opt + ∥p -p ′ ∥ 1 = k opt + c -p x -p y ⊓ ⊔ Proposition 4. In figure
C(q) = k opt + ∥a -q∥ 1 = k opt -c 1 -q x + q y c 1 = a y -a x
Similarly it can be proven that for q ∈ B:

C(q) = k opt + ∥b -q∥ 1 = k opt + c 2 + q x -q y c 2 = b y -b x

Corollary 2. According to propositions 3 and 4 we can conclude that points in

A and C which are on segments parallel to segment ab have the same 1-center objective value. Also for B and D these points are on segments perpendicular to ab. Moreover, points on ab are optimal of 1-center objective (Figure 7b).

Pareto Optimal Solutions

Suppose M opt and C opt are calculated. Obviously if they have intersection, it is the set of Pareto solutions. Hence in the rest of this section we assume that M opt and C opt have no intersection. Fig. 7: Computing 1-center optimal using FVD.

Possible Region for Pareto Optimal Set

Here the goal is to find the region P such that its boundary points dominate all points of the plane, ie. Pareto set is definitely in P.

According to the optimal of M (p) and C(p), three cases are possible. In the first case M opt is in regions A or C, in the second case M opt is in regions B or D and in the third case M opt intersects with the axis aligned the edges of FVD. For the first case (Figure 8a) let e be the lower left point of M opt and let ef and ed be the and horizontal segments hitting the edges of FVD. For all points u on line of ed and w on half-line segment ℓ 1 perpendicular to ℓ ed , C(u) < C(w) and M (u) ≤ M (w). Thus u dominates all points on ℓ 1 . Similarly for point q on ℓ ef and w on half-line segment ℓ 2 , q ≺ w. There are similar results for other edges of adef b which make us able to conclude that polygon adef b is P.

Second and third cases are similar and we consider them simultaneously (Figure 8b and Figure 8c). Let bcde be in region B. Obviously above discussion holds for points p, q, r, s and half-line segments ℓ 1 , ℓ 2 , ℓ 3 and ℓ 4 respectively. Moreover, a dominates all points of D, any point t on ab dominates all points on horizontal (resp. vertical) half-line segment which starts from t and pass through C (resp. A) and b dominates all points on ab. Therefore, we can conclude that points on the border of bcde dominate all points outside of it and bcde is P. It is the same when bcde is in D. ⊓ ⊔ Proposition 6. Suppose P is in A (resp. C). In P let q be a point in a cell with equation M (q) = αq x + αq y + c such that α < 0 (resp. α > 0). Suppose ℓ be the line passing through q with equation y = -x + c ′ . By extending proposition 5, for all p on ℓ or bellow (resp. on ℓ or above) we have M (q) ≤ M (p). The same result holds for B and D when ℓ is y = x + c ′ .

Pareto Optimal Solutions

The following lemma introduces special cells in P which are part of Pareto set.

In the rest of this paper we refer to them as Pareto cells. Proof. Here we assume P ⊂ A but the proof is similar when P ⊂ C. Consider p ∈ P such that M (p) = αp x + αp y + c and α < 0. Suppose q dominates p and ℓ be a line passing through p with y = -x + c ′ equation. if q is above ℓ then C(q) is greater than C(p). Therefore, q is on or bellow ℓ. By proposition 6, if q is bellow ℓ it means M (p) < M (q) otherwise q is on ℓ; but if both are in the same cell it concludes that C(p) = C(q) and M (p) = M (q), otherwise M (p) < M (q). We can obtain from these contradictions that no point dominates p. ⊓ ⊔ Similar to lemma 4, cells with M (p) = αp x + βp y + c such that β = -α and α > 0 (resp. α < 0) are Pareto cells in region B (resp. D).

For intersection of a Pareto cell with edges of FVD several cases are possible. If the Pareto cell intersects with a horizontal (resp. vertical) edge, segment from b (resp. a) to border of the Pareto cell will be the rest of Pareto solutions, we refer to this segment as Pareto segment. Suppose point q dominates p ∈ Pareto segment and let ℓ be the line passing through p and parallel to y = -x, then q must be on or below this line, otherwise C(q) > C(p). But if q is on or below ℓ, since p is in a cell that α/β > 1, M (p) < M (q). If Pareto cell intersects with ab, the part of cell which is in P is also Pareto cell (Figure 9).

Conclusion and Future Work

In this paper we introduced an important and useful multi-objective optimization problem with 1-median and 1-center in Manhattan metric as its objectives. We considered the problem in one and two dimensional space. We also determined the Pareto optimal front and Pareto set simultaneously. Furthermore we proved finding Pareto front and Pareto solution set of proposed problem is θ(n log n). In higher dimensions, considering Manhattan metric, similar to two dimensional space we can show that optimal of 1-median, ie. M (x), will be a d dimensional hypercube. Also, it can be computed in O(dn). For optimal of 1-center, ie. C(x), the propositions are not straight forward. However finding the smallest circumferential hypercube drives us to the hyperplane which is the locus of cube's center (optimal of 1-center). Moreover, it seems that farthest point Voronoi diagram has the most 2d regions. Thus we guess the Pareto optimal set is very similar to two dimensional space; ie. smallest interval of hypercubes from M (x) to C(x) which are connected by their corners in direction perpendicular to locus of optimal of C(x).

In Euclidean metric, we think this problem will be much harder and the Pareto solutions cannot be computed exactly. In this case we have to approximate Pareto solutions and Pareto front. Moreover, this approximation can be followed for harder objectives such as 2-median and 2-center.

Definition 1 .

 1 1-median 1-center two-objective optimization problem: Let P = {p 1 , . . . p n } be a set of demand points in R d . Consider functions M (u) = ∑ n i=1 D(u, p i ) and C(u) = max 1≤i≤n D(u, p i ) are the values of point u ∈ R d as a center for 1-median and 1-center objectives respectively for a certain distance function D. The goal is finding u * to minimize the objectives.

  (a) copt and Mopt have intersection.(b) copt is on the right side of Mopt.(c) copt is on the left side of Mopt.

Fig. 2 :

 2 Fig. 2: Pareto set computation in one dimension.

  ordered break points of the piecewise linear Pareto front function-besides the Pareto solutions interval in o(n log n) running time. Let A = {a 1 , . . . , a n } is the set of input values of sorting problem, l = arg min 1≤i≤n a i and h = arg max 1≤i≤n a i . Suppose b 1 , . . . , b n+1 and t are values such that b 1

Theorem 1 .

 1 Algorithm 1 compute one dimensional 1-median 1-center two-objective Pareto front and Pareto solutions interval in O(n • log n).

Fig. 3 :

 3 Fig. 3: One dimensional 1-median 1-center two-objective Pareto optimal front.

Fig. 4 : 1 -

 41 Fig. 4: 1-median optimal and equation of M (p) in middle cells.

  Farthest point Voronoi diagram regions in Manhattan metric. Possible region for site of C is A \ E.

Fig. 5 :

 5 Fig. 5: Farthest point Voronoi diagram properties.

1 (

 1 b) Projection of points to line ℓ1.

Fig. 6 :

 6 Fig. 6: Compute C(p) using FVD

  We have shown that M (p) partitions the plane to cells in which equation of M (p) is known. According to this partitioning and C opt , seven cases are possible. First three cases happen when M opt is in A or C of FVD. Next three cases occur when M opt is in B or D. Last case occurs when M opt and axis aligned edges of FVD have intersection. Mopt is on the axis aligned edge of F VD.

Fig. 8 :Proposition 5 .

 85 Fig. 8: Different cases for possible region of Pareto solutions.

Lemma 4 .

 4 Suppose P is in A (resp. C). All points like p of cells with M (p) = αp x + βp y + c such that α = β and α < 0 (resp. α > 0), are all or part of Pareto set.

Lemma 5 .Theorem 2 .Fig. 9 :Corollary 3 .

 5293 Fig. 9: Intersection of Pareto cells with edges of FVD

  is O(log n) running time. Therefore, we can conclude that Algorithm 1 is O(n • log n).

				⊓ ⊔
	Corollary 1. Pareto front of one dimensional 1-median 1-center two-objective
	optimization problem can be computed in θ(n log n).
	Algorithm 1 Compute Pareto Optimal Front
	Input: Set I s.t. |I| = n
	Output: Ps(Pareto solutions), P f (Pareto front)
	1: Sort I increasingly to {x1, x2, . . . , xn}
	2: if n is even then
	3:	b = n 2 + 1
	4: else
	5:	b = n+1 2
	6: end if
	7: Ps = [x b , (x1 + xn)/2]
	8: P f = Φ 9: Add ( C(x b ), M (x b ) )	to P f
	10: i = b
	11: while xi+1 < (x1 + xn)/2 do
	12: 13:	i = i + 1 Add ( C(xi), M (xi) )	to P f
	14: end while 15: Add ( C(x (x 1 +xn)/2 ), M (x (x 1 +xn)/2 ) )	to P f
	16: return Ps, P f

  of m 1 , m 2 ∈ A from line ℓ 1 is equal to their distances from segment ab. For point m 1 both distances are obviously the same and are equal to ∥m 1 -m ′ 1 ∥ 1 . For point m 2 we have ∆pqm ′

	2 2 and qa are 2 a as equal isosceles triangles. Therefore, segments qm ′ and ∆qm ′
	equal. Hence ∥m 2