Florent Lallier
email: lallier@gocad.org

Guillaume Caumon

Jean Borgomano

Sophie Viseur

Francois Fournier

Christophe Antoine

Théophile Gentilhomme

Relevance of the stochastic stratigraphic well correlation approach for the study of complex carbonate settings: Application to the Malampaya buildup (Offshore Palawan, Philippines)

The stochastic stratigraphic well correlation method considers the stratigraphic correlation of well data as a set of possible models to sample and manage uncertainty in subsurface studies. This method addresses the incompleteness of typical subsurface data such as limited seismic resolution, seismic blindness due to the lack of impedance contrast between distinct stratigraphic formations, borehole preferential location.

The stochastic stratigraphic well correlation method is applied to the Malampaya buildup (a well documented offshore gas field located North-West of the Palawan Island, Philippines), aged upper Eocene to lower Miocene, and developed on the crest of a tilt-block. Among the available data, ten wells, seven of which are cored, have been drilled and a high resolution 3D seismic survey was acquired by Shell Philippines in 2002.

Previous studies highlight that rock petrophysical properties are mainly controlled by diagenesis. Correlation rules are thus developed in order to adapt the stochastic stratigraphic well correlation method to the study of diagenetic units. These rules are based on wireline log shape and diagenetic units types.

Four stratigraphic correlation models are generated using the proposed correlation method: a deterministic one corresponding to the most probable model considering only well data and three stochastic ones. These correlation models are bound with geostatistical methods to build static reservoir models. Synthetic seismic profiles are computed from facies models conditioned to acoustic impedance models. It leads to comparable seismic amplitude images, highlighting the importance of considering several well correlation models for one given seismic survey. Stochastic stratigraphic correlations are shown to have a first order impact on reservoir unit characterization, rock volumes and fluid flow response on the reservoir model.

The Malampaya gas field is a well documented carbonate buildup located North-West of the Palawan Island, Philippines [START_REF] Grötsch | Integrated 3-D reservoir modeling based on 3-D seismic: The Tertiary Malampaya and Camago buildups, offshore Palawan, Philippines[END_REF][START_REF] Neuhaus | Quantitative seismic reservoir characterization of an Oligocene-Miocene carbonate buildup: Malampaya field, Philippines[END_REF][START_REF] Fournier | Paleoenvironments and highfrequency cyclicity from Cenozoic South-East Asian shallow-water carbonates: A case study from the Oligo-Miocene buildups of Malampaya (Offshore Palawan, Philippines)[END_REF][START_REF] Fournier | Development patterns and controlling factors of Tertiary carbonate buildups: Insights from high-resolution 3D seismic and well data in the Malampaya gas field (Offshore Palawan, Philippines)[END_REF][START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF](Fig. 1). Reservoir rocks belong to the upper Eocene to lower Miocene aged Nido Limestone formation [START_REF] Sales | The petroleum potential of deep-water northwest Palawan Block GSEC[END_REF]Williams, 1997). Deposition of the Nido limestone occurred on the edge of a tilted block resulting from the rifting phase of the China sea (upper Eocene to lower Oligocene). Development history, paleo-environments and diagenesis of the Malampaya buildup has been studied by [START_REF] Fournier | Development patterns and controlling factors of Tertiary carbonate buildups: Insights from high-resolution 3D seismic and well data in the Malampaya gas field (Offshore Palawan, Philippines)[END_REF][START_REF] Fournier | Paleoenvironments and highfrequency cyclicity from Cenozoic South-East Asian shallow-water carbonates: A case study from the Oligo-Miocene buildups of Malampaya (Offshore Palawan, Philippines)[END_REF]. Meter-scale, 4 th to 5 th order subtidal cycle bounded by exposure sur-Author's copy of doi:10.1144/SP370.12 published in: Garland, J., Neilson, J. E., Laubach, S. E. & Whidden, K. J. (eds) 2012. Advances in Carbonate Exploration and Reservoir Analysis. Geological Society, London, Special Publications, 370, c The Geological Society of London 2012. [START_REF] Grötsch | Integrated 3-D reservoir modeling based on 3-D seismic: The Tertiary Malampaya and Camago buildups, offshore Palawan, Philippines[END_REF] and [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF].

faces are interpreted as constitutive of the buildup. [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF] studied petrophysical properties of the Nido limestone of the Malampaya buildup and show that reservoir properties are mainly controlled by diagenetic processes affecting carbonate rocks. They defined five types of diagenetic units (Ia, Ib, II, IIIa, IIIb) characterized by: coherent set of diagenetic transformations (meteoric diagenesis and late burial cementation and leaching); specific statistical porosity and acoustic impedance distributions (Fig. 1); pore types; porosity-permeability relationships; and well-log signatures. From well-to-seismic ties, [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF] showed that the main seismic reflectors within the buildup image the boundary of diagenetic units. The 3D picking of the reflectors and well correlation of diagenetic units thus provide an envelope of geobodies characterized by specific pore type and reservoir properties. Diagenetic units also define the layering of the reservoir and impact the 3D architecture, showing in particular an alternation of porous and tight diagenetic units. Hence, it impacts the understanding of the reservoir flow be-haviour on fluid flow modelling. [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF] proposed a stratigraphic model of seismo-diagenetic units of the Malampaya buildup on the basis of well logs and interpretation of an high resolution seismic survey(Fig. 1d). However, in common subsurface studies the well-to-well stratigraphic correlation process is subject to uncertainties due to limited seismic resolution, ambiguity between rock types and well logs and inherent variability between boreholes. The objectives of this study are to question if the seismic response of the Malampaya buildup may be explained by alternative well-to-well stratigraphic correlation models of diagenetic units and to assess the impact of well correlation uncertainty on reservoir description and dynamic behaviour. Different possible sets of diagenetic units are automatically generated by a stochastic well correlation process and used to produce a static and dynamic model of the Malampaya buildup.

Different approaches and well descriptors (markers, facies, wireline logs, etc.) have been Author's copy of doi:10.1144/SP370.12 published in: Garland, J., Neilson, J. E., Laubach, S. E. & Whidden, K. J. (eds) 2012. Advances in Carbonate Exploration and Reservoir Analysis. Geological Society, London, Special Publications, 370, c The Geological Society of London 2012.

proposed to perform computer-aided stratigraphic well correlations. For instance, deterministic methods have been proposed to correlate wells on the basis of logs (Zoraster, 2004) or lithostratigraphic description [START_REF] Griffiths | Interwell matching using a combination of petrophysically derived numerical lithologies and genetyping techniques[END_REF]. Recently, [START_REF] Lallier | 3D stochastic stratigraphic well correlation of carbonate ramp systems[END_REF] have presented a stochastic procedure for generating several plausible stratigraphic well correlations from a given set of wells, described by depositional facies.

In the first part of this study, an adaptation of the stochastic method proposed by [START_REF] Lallier | 3D stochastic stratigraphic well correlation of carbonate ramp systems[END_REF] to the stratigraphic correlation of diagenetic units is presented. It relies on the correlation rules proposed by [START_REF] Fang | Computer-aided well log correlation[END_REF]. Using this proposed approach, four stratigraphic correlation models of diagenetic units are generated and used to construct cross sectional models of these units. Geostatistical simulations of porosity, permeability and acoustic impedance are then performed. The impact of stochastic stratigraphic correlation is then evaluated on the basis of synthetic seismic computation and fluid flow modelling.

Stratigraphic correlation method

Correlation rules

Two rules are used here to perform stratigraphic well correlation of diagenetic units identified by [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF] from well logs, thinsection, cuttings, cores and well logs: (1) Two units are correlated only if they are of the same diagenetic type. This rule implies that diagenetic units are correlated as geological bodies; (2) Two units are correlated if they display a comparable well log pattern as explain below.

Evaluation of the value of a correlation between two units

An automatic method of well correlation requires a mathematical formulation of correlation rules. Several methods have been proposed to evaluate the similarity between two well log curves. [START_REF] Olea | Correlator 5.2 -a program for interactive lithostratigraphic correlation of wireline logs[END_REF] used the product between standardized similarity coefficient and Pearson's correlation coefficient to evaluate similarity between two parts of a wireline log. These two coefficients suppose that for the two considered units the same number of points are sampled, hence may call for re-sampling the well log for the considered unit. Moreover in the case of lithostratigraphic correlation, the use of this correlation coefficient assumes a globally constant preservation rate between the two wells (i.e no distortion in the signal). The sedimentary record is also assumed to be quasicomplete between the two units, i.e there is no erosion on a studied stratigraphic unit which is not recorded in the other unit. In this study, we use the method proposed by [START_REF] Fang | Computer-aided well log correlation[END_REF] to evaluate similarity between well logs of two diagenetic units. This method, based on the dynamic programming approach, follows a workflow with three steps (Fig. 2): Transformation of well logs into pattern primitives. Well log curves are filtered to extract points that are representative of a log pattern. This method is called "adaptive" by [START_REF] Fu | Digital pattern recognition[END_REF] and is performed in two steps. First, a generalized difference log d is computed along the well at each point at depth i:

d 2 i = (s i-α -s i+α ) 2 (1)
where d i is the value of the generalized difference log d at depth i; s i-α the slope of the considered well log between the point at depths i and i -α; and s i+α the slope between the point at i and i+α.

The half-window α is to be defined according to the well log resolution and noise characteristics.

Then for each peak of the generalized difference log, the value of the studied log at peak depth is selected and used to build a filtered log. Finally, the obtained simplified log is transformed into a string of pattern primitives according to the dip angle of the filtered log segment. Similarity evaluation. Considering two stratigraphic units a and b transformed into two strings of pattern primitives (respectively [a 1 , ..., a i , ...a n ] and [b 1 , ..., b i , ...b m ]), the best match is searched between a and b using dynamic programming [START_REF] Levenshteiti | Binary codes capable of correcting deletions, insertions, and reversals[END_REF] where the penalty of an insertion and a deletion equals 1 and the penalty of a mismatch is equal to the difference between a i and b j . Finally, the similarity between the two units a and b has a value C which equals 0 if a and b are similar and increases with the difference between a and b:

C = T (4 × L) (2)
Author's copy of doi:10.1144/SP370.12 published in: Garland, J., Neilson, J. where L is the mean length of the strings and T the cumulative penalty of association between a and b. This method has two main advantages:

(1) it handles possible distortions of the initial signal through a string of pattern primitives; (2) it manages a missing events in a unit with regard to the other unit through the use of dynamic programming.

Most of the time, more than one well log is available. Given n wireline logs noted λ 1 to λ n (for instance, a density log, a porosity log and a sonic log), the possibility of a correlation between two considered units can be evaluated using a weighted sum of cost C λi computed for each available well log:

C λ 1 ,...,λn = n i=0 w i C λ i (3)
where w i are weights defined according to: (i) the redundancy of information contained in the wireline log; (ii) the relevance of the rock property recorded by the log.

Automatic stratigraphic correlation between two wells

Once the similarity between stratigraphic units has been computed, the global correlation between the considered wells can be built. The goal of this study is to handle uncertainties in stratigraphic well correlation and to stochastically generate correlation models constrained by well logs. [START_REF] Lallier | 3D stochastic stratigraphic well correlation of carbonate ramp systems[END_REF] proposes the use of a stochastic variation of the Dynamic Time Warping algorithm (DTW) [START_REF] Levenshteiti | Binary codes capable of correcting deletions, insertions, and reversals[END_REF][START_REF] Myers | Comparative study of several dynamic time-warping algorithms for connected-word recognition[END_REF] which automatically builds correlations between two ordered series. The correlation between two wells u and v, on which stratigraphic units u 1 to u n and v 1 to v m are defined, is determined with a recursive function D(u i , v j ):

D u i , v j = S   α = m(u i , v j ) + D u i-1 , v j-1 , β = g(u i ) + D u i-1 , v j , γ = g(v j ) + D u i , v j-1   (4) where:
(a) D(u i , v j ) is the cost of the correlation between units u 1 to u i of w 1 and units v 1 to v j of w 2 ;

(b) m(u i , v j ) is the cost of a correlation (match) between the units u i and v j ;

Author's copy of doi:10.1144/SP370.12 published in: Garland, J., Neilson, J. E., Laubach, S. E. & Whidden, K. J. (eds) 2012. Advances in Carbonate Exploration and Reservoir Analysis. Geological Society, London, Special Publications, 370, c The Geological Society of London 2012.

(c) g(u i ) is the cost of a gap of unit u i , i.e. the unit u i of the well u is not recorded on well v and end as a stratigraphic gap between the two wells;

(d) g(v i ) is a cost of a gap of units v j of v;

(e) S(α, β, γ) is a value equal to cost α (respectively β or γ) with a probability inversely proportional to the relative cost α α+β+γ (respectively

β α+β+γ or γ α+β+γ ).
Applying equation 4 from i = 1 to n and j = 1 to m minimizes the cost D(u n , v m ) of the correlation between wells u and v and thus builds correlations which are consistent with rules used to compute m(u i , v j ), g(u i ) and g(v i ). In the original deterministic version of the DTW [START_REF] Myers | Comparative study of several dynamic time-warping algorithms for connected-word recognition[END_REF], the S function is not a stochastic one, and returns the minimum value of α, β and γ ensuring that the DTW algorithm builds the minimum cost correlation between wells u and v. This correlation is called deterministic.

In this study, the correlation between two diagenetic units is considered only if they are of the same diagenetic type. Considering two units u i and v j of a diagenetic type respectively t(u i ) and t(v j ), the cost of a match is computed as following: ..,λn , as inFang et al. (1992) (5) The analysis of the match function C λ1,...,λn highlights that the cost of a match between two units has a cost lower than 1 if log patterns of u i and v j are similar and higher than 1 if no similarity appears between the considered units. Therefore, we set the gap costs g(u i ) and g(v j ) to 1 so that units with comparable log patterns have a higher chance of correlation than pinch out.

m(u i , v j ) = ∞ if t(u i ) = t(v j ) else m(u i , v j ) = C λ1,.

Results and discussions

Stratigraphic correlation models

Four stratigraphic correlation models (Fig. 3) are built between wells MA-1, 2, 5, 7 and 8: a stratigraphic correlation model (DSM) corresponding to the minimum cost correlation (S(•) ≡ M in(•)) and three stochastic stratigraphic correlation models (SSM1, 2 and 3).

Validity of the well correlation method. The confrontation of our DSM (i.e. the most probable stratigraphic correlation model regarding used correlation rules, Fig. 3a) with the stratigraphic model proposed by [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF] (Fig. 1d) serves as validation of the method. Several points can be highlighted: (i) all of the diagenetic units correlated by [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF] are also correlated in the DSM; (ii) hiatuses interpreted by [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF] are also present in the DSM; (iii) the only difference between the two models is the interpretation of the units marked with an arrow on figure 1d which are interpreted as pinching out diagenetic lenses by [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF](Fig. 1B) and top-lap in the DSM. This difference comes from a modelling decision we made: the presence of diagenetic lenses has been interpreted by [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF] based on a high resolution seismic survey on which the abrupt disappearance of seismic reflectors (labelled M20.0 and M30.0 on figure 1b andd) corresponding to the diagenetic units indicates pinch out. In our case, only well log data on which no indication of such pinch out could be identified, are considered to preform stratigraphic well correlation. As a consequence, we did not design our stratigraphic correlation algorithm to manage units that pinch out. However, the likeness between the DSM and the stratigraphic model proposed by [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF] suggest that our algorithm generates an acceptable stratigraphic well correlation.

Geometrical model building. The four possible stratigraphic correlation models are combined with top and bottom Nido horizons geometry interpreted from 3D seismic to build 3D stratigraphic grids [START_REF] Mallet | Geomodeling. Applied Geostatistics[END_REF] conforming to diagenetic units (sections of these grids are shown in Fig. 3). In addition to the uncertainties due to the correlations, there are uncertainties due to the way horizons are interpolated between the correlated markers [START_REF] Goff | Simulation of stratigraphic architecture from statistical and geometrical characterizations[END_REF][START_REF] Caumon | Elements for stochastic structural perturbation of stratigraphic models[END_REF]Seiler et al., 2010). To limit the influence of this geometric uncertainty, we focus on the neighbourhood of wells MA-1, 2 and 5 (Fig. 1a). This restriction is necessary since the internal geometry of the reservoir is only constrained by well marker positions and by the top and bottom reservoir horizons.

Author's copy of doi:10.1144/SP370.12 published in: Garland, J., Neilson, J. E., Laubach, S. E. & Whidden, K. J. (eds) 2012. Advances in Carbonate Exploration and Reservoir Analysis. Geological Society, London, Special Publications, 370, c The Geological Society of London 2012. Author's copy of doi:10.1144/SP370.12 published in: Garland, J., Neilson, J.

Seismic response to alternative stratigraphic models

The value of a stochastic correlation model can be questioned when high quality seismic data allows us to follow stratigraphic units between wells. Moreover, since the petrophysical properties distribution through the Nido Limestone is controlled by the geometry of tight and porous diagenetic units imaged by seismic survey, we could question if the seismic response of the Malampaya buildup can be explained by alternative internal architectures. Synthetic seismic computation. A realisation of acoustic impedance is presented for each correlation model in figure 4. Synthetic seismic amplitude models are computed by convolving the acoustic impedance models (DSM, SSM1, 2 and 3) with a Ricker 50Hz wavelet extracted from the seismic data by [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF](Fig. 4). A realisation of convolution is presented for each correlation model in figure 5.

Qualitative comparison of synthetic seimic. Synthetic seismic amplitude models are visually compared with MA-1 to MA-2 section of high resolution 3D seismic presented by [START_REF] Fournier | Development patterns and controlling factors of Tertiary carbonate buildups: Insights from high-resolution 3D seismic and well data in the Malampaya gas field (Offshore Palawan, Philippines)[END_REF] and [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF] (Fig. 1b). Synthetic seismic shows a good visual correlation with original reflectivity data. The major features of the original high resolution seismic data (Fig. 5A) are present on the synthetic seismic section: top and bottom reservoir reflectors, major internal reflectors and antiforms located on the eastern part of the buildup [START_REF] Fournier | Development patterns and controlling factors of Tertiary carbonate buildups: Insights from high-resolution 3D seismic and well data in the Malampaya gas field (Offshore Palawan, Philippines)[END_REF]. This visual comparison allows us to validate our reservoir model on a first order. When comparing all of the synthetic seismic models together, it is hard to discriminate which geometrical model (Fig. 3) or acoustic impedance model (Fg.4) has been used to generate the synthetic seismic model in figure 5. It appears that even when good quality reflection data are available, internal architecture of the Malampaya buildup can be explained by several stratigraphic correlation models of diagenetic units.

Quantitative comparison of synthetic seismic. The normalized root mean squared difference (NRMS) of the MA-1 to MA-2 cross-section is computed to evaluate the likeness between synthetic seismic models. Considering one reference seismic model r and one alternative synthetic seismic model a, the NRMS at a seismic trace between top (t) and bottom (b) reflectors is defined by:

N RM S = 2×RM S(r i -a i ) RM S(r i )-RM S(a i ) with RM S(x i ) = i=b i=t (x i ) 2 N (6)
where N is the number of samples between t and b.

We define the synthetic seismic computed from the convolution of one realisation of acoustic impedance on the DSM as our reference seismic. From ten simulations of acoustic impedance, ten synthetic seismic models are computed on each stratigraphic correlation model (DSM, SSM1, 2 and 3). Figure 6 shows the NRMS computed between our reference model and each of the forty alternative synthetic seismic models. The mean value of NRMS of each stratigraphic correlation model (Fig. 6b) shows that the difference between the reference seismic model and DSM-based alternative models is comparable with the difference between the reference model and SSM-based seismic models. Moreover, figure 6c Author's copy of doi:10.1144/SP370.12 published in: Garland, J., Neilson, J.

Implication on fluid flow modelling

Flow simulation are performed on the MA-1 to MA-2 cross-section model, using the GPRS flow simulator, to evaluate the impact of alternative stratigraphic frameworks on flow response. The porosity is simulated with SGS conditioned by wells MA-1, 2, 5, 7 and 8, using the same random walk as that used to simulate acoustic impedance and synthetic seismic models (Fig. 4 and Fig. 5 respectively). The permeability has been modelled using porosity-permeability crossplots presented by [START_REF] Fournier | Geological significance of seismic reflections and imaging of the reservoir architecture in the Malampaya gas field (Philippines)[END_REF]. To highlight the influence of correlation uncertainty on the recovery in a classical reservoir setting, we consider: (i) Malampaya as a black oil reservoir without an aquifer; (ii) a water injector close well to well MA-2; (iii) a productor well close to well MA-1; (iv) an initial oil saturation of 0.8. Figure 7 shows permeability models and oil saturation at four different time steps. The visual evaluation of saturation at different time steps shows that the use of alternative stratigraphic correlation models significantly impacts connectivity between porous/permeable units. This results in different drainage areas and may impact decision-making in reservoir development. The stratigraphic correlation models generated with our approach can be used to generate an a priori set of models representing reservoir stratigraphic uncertainties. A model screening method based on production data (Suzuki et al., 2008) could be then used to efficiently select acceptable stratigraphic correlation models with regard to observed reservoir flow behaviour.

Conclusions

The stratigraphic well correlation of units identified along section is a hazardous task even if high quality data are available. This study shows that different stratigraphic correlation models may produce similar seismic images, and hence several stratigraphic well correlations models may be considered from the initial data (including in, our case, 5 wells and seismic data). As shown by the computed synthetic seismic models, seismic data may not discriminate stratigraphic correlation models. However, stratigraphic correlation models have a significant impact on fluid flow. This paper demonstrates the necessity and the technical feasibility of a numerical management of uncertainty on stratigraphic well correlation. Additionally, model screening approaches could be used to select only stratigraphic models compatible with production data or well tests. Limits of the proposed approach should be pointed out; in particular, the inability of our stratigraphic well correlation method to deal with diagenetic lenses. A way to manage this issue could be found in using improvement of the DTW algorithm proposed by Waterman & Raymond Jr. (1987). In the version of the DTW we use, a stratigraphic unit can be correlated to only one other unit ("one to one correlation") excluding diagenetic lenses. The "one to many correlation" configuration proposed by Waterman & Raymond Jr. (1987) makes it possible to correlate several successive stratigraphic units of one well to one unit on another well and thus to manage lenses. However, such a configuration calls for computing a correlation cost associated to the "one to many correlation". The mathematical expression of this cost is not easy to formulate honestly and requires additional development. As highlighted by [START_REF] Bashore | Importance of a geological framework and seismic data integration for reservoir modeling and subsequent fluidflow predictions[END_REF], the choice of the well correlation strategy affects not only the geometrical model used to guide the geostatistical simulations, but also dynamic reservoir model and the associated fluid flow predictions. This study corroborates their conclusions and points out that even when stratigraphic well correlation rules are defined, several well correlation models have to be considered. In the presented case study, diagenetic units are correlated using lithostratigraphic rules. However, in other geological settings such as carbonate ramps, significant lateral facies variability may exist. In such a case, correlation rules presented in this study are unsuitable. The DTW algorithm should then be associated to correlation rules managing lateral facies transitions as proposed by [START_REF] Lallier | 3D stochastic stratigraphic well correlation of carbonate ramp systems[END_REF].

Figure 1 :

 1 Figure 1: (a) Topographic map (in meter sub-sea) of the top horizon of the Malampaya gas reservoir located offshore western Palawan island, Philippines. (b) Seismic line between wells MA-1 and MA-2 displaying Malampaya interior reflector. (c) Porosity and acoustic impedance histogram for the different diagenetic unit types. (d) Stratigraphic model proposed by Fournier & Borgomano (2007) on the basis of wells and seismic interpretation. Modified after Grötsch & Mercadier (1999) and Fournier & Borgomano (2007).

Figure 2 :

 2 Figure2: Evaluation of similarity between log responses of two units. Well logs are first filtered using generalized difference (a) and transformed into strings of pattern primitives according to dip angle of filtered log segments (c). Cost is computed as the difference between the two strings representing the two considered units (b). After[START_REF] Fang | Computer-aided well log correlation[END_REF]

Figure 3 :Figure 4 :

 34 Figure 3: Sections of the diagenetic units model built from stratigraphic correlations computed using the proposed algorithm. The deterministic model is the minimum cost one. For reference, wells MA-1, 2 and 5 are projected on the section.

  and 6d show that the NRMS variability is similar for DSM-based and SSM1-based models, both reaching minimum and maximum values of the envelope of all the realizations. Dispersion at the well location observed on figure 6 are explained by well deviation Author's copy of doi:10.1144/SP370.12 published in: Garland, J., Neilson, J. E., Laubach, S. E. & Whidden, K. J. (eds) 2012. Advances in Carbonate Exploration and Reservoir Analysis. Geological Society, London, Special Publications, 370, c The Geological Society of London 2012. Distance along seismic line (m) Distance along seismic line (m) Distance along seismic line (m) Distance along seismic line (m

Figure 5 :Figure 6 :

 56 Figure 5: MA-1 to MA-2 cross-section of synthetic seismic amplitude model computed by convolution of acoustic impedance model (Fig.4) with a 50hz Ricker wavelet. All of the four seismic sections display comparable internal reflectors and may result in similar stratigraphic interpretation.

  E., Laubach, S. E. & Whidden, K. J. (eds) 2012. Advances in Carbonate Exploration and Reservoir Analysis. Geological Society, London, Special Publications, 370, c The Geological Society of London 2012.

	IniƟal Log	Generalized difference	Filtered Log	String transform	String transform	IniƟal Log Filtered Log
	Depth

  E., Laubach, S. E. & Whidden, K. J. (eds) 2012. Advances in Carbonate Exploration and Reservoir Analysis. Geological Society, London, Special Publications, 370, c The Geological Society of London 2012.

	Property modelling
	3D stratigraphic grids are populated with porosi-
	ties property using Sequential Gaussian Simula-
	tion (SGS) conditioned by well porosity logs on
	wells MA-1, 2, 5, 7, and 8. As shown by Fournier
	& Borgomano (2007), petrophysical properties are
	controlled by diagenetic history and thus diage-
	netic unit type. Therefore, in this study, petro-
	physical properties are considered second order
	stationary on each diagenetic unit type. Univari-
	ate, bivariate and spatial analyses are performed
	for each diagenetic type. Experimental histograms
	(Fig. 1c) are used as target distributions for the
	SGS. Due to lack of information (data are pro-
	vided by five wells), variogram model can not
	be discriminated between diagenetic units as ex-
	pected theoretically. A unique spherical variogram
	model is thus built and used to perform SGS with
	a range of 950 m horizontally and 40 m vertically.
	Porosity is then simulated on each stratigraphic
	grid. The same process is used to simulate acous-
	tic impedance.

  E., Laubach, S. E. & Whidden, K. J. (eds) 2012. Advances in Carbonate Exploration and Reservoir Analysis. Geological Society, London, Special Publications, 370, c The Geological Society of London 2012.

	and grid rescaling between acoustic impedance
	simulation and synthetic seismic computation.
	From this quantitative analysis, it appears that
	synthetic seismic data computed from different
	stratigraphic correlation models are comparable.
	One of the consequences of this observation is that
	several stratigraphic correlation models should
	ideally be considered when performing seismic in-
	terpretation.

Stochastic well correlationF. LALLIER et al.

Author's copy of doi:10.1144/SP370.12 published in: Garland, J., Neilson, J. E., Laubach, S. E. & Whidden, K. J. (eds) 2012. Advances in Carbonate Exploration and Reservoir Analysis. Geological Society, London, Special Publications, 370, c The Geological Society of London 2012.

This work was performed in the frame of the Gocad research project. The companies and universities members of the GOCAD consortium are acknowledged for their support. Authors also thank CGGVeritas for providing convolution software, Paradigm for providing Gocad software and API and Stanford University for providing the GPRS flow simulator. This is CRPG contribution number 2165.

Author's copy of doi:10.1144/SP370.12 published in: Garland, J., Neilson, J. E., Laubach, S. Author's copy of doi:10.1144/SP370.12 published in: Garland, J., Neilson, J.


Florent Lallier
email: lallier@gocad.org

Christophe Antoine
email: antoine@gocad.org

Julien Charreau
email: charreau@crpg.cnrs-nancy.fr

Guillaume Caumon
email: guillaume.caumon@univ-lorraine.fr

Jeremy Ruiu
email: ruiu@gocad.org

Management of ambiguities in magnetostratigraphic correlation $

Keywords: magnetostratigraphy, computer method, uncertainty

Magnetostratigraphy is a powerful tool to provide absolute dating of sediments enabling robust and detailed chronostratigraphic correlations. It is based on the correlation of a magnetic polarity column, observed and measured in a given sediment section, to a magnetic polarity reference scale where polarity changes are well dated via other independent methods. However, magnetostratigraphic correlations are loose as they are only constrained by binary magnetic chrons (i.e normal or reversal) and their thickness, which are both defined from depth variations of the magnetic remanent directions. The thickness of a given magnetic polarity zone is a function of time and sediment accumulation rate, which may not be stationary, leading to ambiguities when performing the correlations. To address these ambiguities, a numerical method based on the Dynamic Time Warping algorithm is proposed. Magnetostratigraphic correlations are computed in order to minimise the local variation of the accumulation rate. The main advantage of the proposed method is to automatically provide a set of reasonably likely correlations. This set can then be scrutinised for further analysis and interpretation. However, the likelihood of a correlation should be handled carefully as it depends on the information content of the magnetostratigraphic section itself and remain ultimately valid by ancillary constraint. Nevertheless, the method gives consistent results on difficult synthetic cases that simulate abrupt variations of the sedimentation rate. Insights on true sections debated by previous authors are also given.

Introduction

Sedimentary basins provide key records of many earth processes including tectonic, climatic and erosional. They may also host large resources of ore deposits (e.g gold, U, etc.), hydrocarbons and fresh water. Accurate and detailed chronostratigraphy is critical to yield quantitative constraints on the wide variety of interactions that may exist between these processes. This observation holds for reservoir and basin modeling which are also based on stratigraphic correlation and request reliable chronology. Sediment dating is therefore a fundamental step for a broad community in Earth sciences. Various techniques have been developed such as radio-isotopic dating (e.g K-Ar, Fission Tracks), biostratigraphy, cosmonucleid isotopes (10Be / 26Al), astronomical records, palynology, 14C, varves etc. All these methods enable the dating of sediments in a large range of environments and for time scales ranging from years to 10 7 years. However, in continental areas, due to the paucity of biostratigraphic markers, reliable dating via classic methods remains challenging. Hence, the depositional ages are mainly inferred from loose lithostratigraphic correlation, which suffers from flaws in the presence of strong lateral facies variations and diachronous deposition (e.g. Charreau et al. (2009b)).Yet, continental sediments may record the Earth magnetic field and its polarity at the time of deposition via the Detrital Remanent Magnetization (DRM) mechanism (e.g [START_REF] Tauxe | Depositional remanent magnetization: Toward an improved theoretical and experimental foundation[END_REF]). Given a sedimentary section or a log, one can measure the variation of the remanent magnetic directions and identify the polarity changes of the Earth magnetic field. A polarity column can then be reconstructed as a succession of polarity zones (or reversal) characterised by their length in meters or duration in years and their polarity (normal or reverse). A normal polarity chron is represented in black in a polarity column and corresponds to a period where the Earth magnetic dipole is similar to the present one with the magnetic North close to the geographic North. A reverse polarity chron is shown in white and corresponds to a period where the dipole has been switched with its magnetic North close to the southern Earth's pole. Magnetostratigraphy is based on the correlation of a measured polarity column in a sediment section with the geomagnetic polarity time scales (GPTS [START_REF] Lourens | Orbital tuning calibrations and conversions for the Neogene Period[END_REF]), where polarity changes are well dated via independent methods (astronomical constraints, spreading of the oceanic floor, thermochronology). This powerful technique enables the dating of sediments in many areas where other approaches have failed. However, magnetostratigraphy is based on several assumptions: (1) the Earth magnetic field is bipolar; (2) the sedimentary section is con-tinuous with no depositional hiatus nor erosion; and (3) the accumulation rates are constant at short time scales (< few Ma). This method is particularly well suited for the Cenozoic Era and Jurassic Period in which numerous polarity changes have been documented. Nevertheless, the major drawback of this dating technique is that the correlation between the observed polarity column and the reference scale is a subjective and qualitative interpretation as only based on eyes, providing a single dating with no assessment of uncertainty. Other datings, for example from fossil or lava flows observed in the studied section, are therefore very useful to independently constrain the correlation and reduce uncertainty. In many cases, the absence of such data leads to controversial debates between authors leaving magnetostratigraphic correlation results to doubts and criticism. For example in the Subei section in Western China, three different authors have proposed different correlations with ages ranging between 33.5 and 27 Ma [START_REF] Yin | Tectonic history of the Altyn Tagh fault system in northern Tibet inferred from Cenozoic sedimentation[END_REF], 26 and 19 Ma [START_REF] Gilder | Oligo-Miocene magnetostratigraphy and rock magnetism of the Xishuigou section, Subei (Gansu Province, western China) and implications for shallow inclinations in central Asia[END_REF], and 20 and 9 Ma [START_REF] Wang | Danghe area (western Gansu, China) biostratigraphy and implications for depositional history and tectonics of northern Tibetan Plateau[END_REF], leading to conflicting interpretations about the tectonoclimatic history of the Northeast Tibet Plateau (section 2.1). A magnetic polarity column only provides a partial view of the sedimentary record. It is therefore almost impossible to find "the right" correlation model, as to do so would require a huge amount of independent data impossible to acquire in practice. Instead, we suggest to automatically generate possible correlation models by using wellestablished rules and principles in order to deal with uncertainty. In the proposed algorithm (section 1.3.1), each correlation is associated with a cost, which is used to rank the various proposed solutions according to specified rules. Moreover, the proposed approach provides an estimate on the instantaneous accumulation rate and their uncertainties on the basis of several possible solutions while classic "eye-based" correlation only gives the 90-95% confidence limit from basic linear regression built on a unique correlation. Tests carried out on synthetic polarity columns (Section 1.4) attest that the algorithm can handle sedimentary hiatus and sharp variations of the accumulation rate. The method is therefore applied to the Subei and Yaha magnetostratigraphic sections in Central Asia (section 2) where a debate still remains on the depositional ages although they represent standard sections for tectonics and climate in this region. The analysis of the 10,000 minimum cost correlations, covering many acceptable possibilities for these sections provides interesting insights about the correlations discussed in the literature.

Computer method for magnetostratigraphic correlation

Problem Settings

Computing correlations between magnetic polarity section and polarity reference scale is challenging since:(i) the magnetic record is binary information (ii) the studied section represents only a part of the whole polarity reference scale and has potentially been transformed by changes in the accumulation and compaction rates and (iii) undersampling may yield unrecorded polarity changes. Moreover, in a sedimentary section, the correlation to the reference scale strongly depends on the thickness of each polarity zone, which is a function of three parameters : (i) the duration of the polarity chrons i(t); (ii) the accumulation rate p(t) and (iii) the compaction rate c. From a mathematical analysis, [START_REF] Lowrie | Geomagnetic polarity timescales and reversal frequency regimes[END_REF] show that the geomagnetic polarity time serie is not stationary but is composed of stationary periods. They also reveal that the duration of the polarity chrons follows a Poisson law. The accumulation rate is also a non-stationary process as it is driven by the nature of the depositional environment and erosion processes, which may rapidly (i.e. few kyr) and strongly varies. However, authors often consider the accumulation rate to change stepwise between long stationary periods. For example, from magnetostratigraphic study of the Jingou He section in northern Tianshan piedmont, Charreau et al. (2009a) found that the mean accumulation rate was 0.10 mm/yr, 0.18mm/yr and 0.29mm/yr from 23 Ma to 16Ma, from 16Ma to 11Ma and from 11Ma to 1Ma, respectively. Finally, the compaction rate may vary along a given section due to rock lithology and burial leading to distortion of the observed accumulation rate. To provide robust, relevant and trustworthy correlations, computer based methods will be required to handle distortion in the time signal of the magnetic chrons induced by these processes (accumulation rate, compaction rate. . . ). (2011) proposes an interesting computer method which correlates a given polarity column to the reference scale under the following assumptions: (i) the thickness of the observed polarity zone is supposed stationary, (ii) the polarity column is supposed complete (all polarity changes have been recorded and sampled in the section), (iii) the age of the studied column is a priori roughly known and (iv) the geomagnetic reversal process is supposed to follow a Gaussian law. As pointed out above, these assumptions can be discussed when working on long sedimentary sections. First, the chron thickness and duration are mainly non-stationary. Second, due to non-depositional and/or erosional periods, several magnetic polarity changes may be missed in the sediment record. Moreover, non-amenable lithology for sampling (e.g conglomerate) may sometimes impose low sampling density insufficient to identify all magnetic reversals. And finally, in continental areas it is sometimes very difficult to assess the age of a given section from sedimentology only, due to the absence of biostratigraphic markers. For example, in Central Asia, thick conglomerate lying at the top of the basin sedimentary piles, have long been assigned to the Plio-pleistocene whereas recent magnetostratigraphic studies have shown that their deposition could start during the middle Miocene (Charreau et al., 2009). Another way to correlate a given polarity column to the reference scale was proposed by [START_REF] Man | On the identification of magnetostratigraphic polarity zones[END_REF]. This method is based on a preprocessing step which transforms both the polarity column and the reference scale to dimensionless data by calculating the logarithm of the thickness (or the duration) of a given magnetic chron minus the loga- rithm of the thickness (or the duration) of the underlying chron. Thanks to this transformation, accumulation rate and chron duration can be supposed locally stationary. The correlation between the reference scale and the studied polarity column is then computed using a cross correlation function. This method is less sensitive to the non-stationarity of the data but still requires that all reversals have been recorded and sampled in the sediments.

Previous work

Man

In this paper, we similarly use thickness ratios to relax the global assumption of stationary preservation rate.

Proposed approach

Our approach has been designed to manage non-stationary data, unrecorded chrons in the polarity column (sedimentary hiatus and/or erosion) and to support independent constraints derived from biostratigraphy or geochemical dating techniques. This method is based on the Dynamic Time Warping (DTW) algorithm described in section 1.3.1 [START_REF] Levenshteiti | Binary codes capable of correcting deletions, insertions, and reversals[END_REF]. The DTW algorithm is modified to output not only the "best" correlation between the studied polarity column and the GPTS, but also the "n best" correlations in order to manage uncertainties on sediment dating and accumulation rate. The Python implementation used in this paper is available at the address http://www. gocad.org/w4/index.php/research/free-software, at which new versions will be posted. Input data are: (i) the reference polarity section which is supposed exhaustive (i.e. there is no polarity change in the recorded section which is not in the reference scale); (ii) the thickness of polarity zones observed in the sediment section.

The Dynamic Time Warping algorithm

The Dynamic Time Warping (DTW) algorithm has been introduced by [START_REF] Levenshteiti | Binary codes capable of correcting deletions, insertions, and reversals[END_REF] and found many applications for automatic signal correlation in speech recognition [START_REF] Myers | Comparative study of several dynamic time-warping algorithms for connected-word recognition[END_REF], bioinformatic and geoscience for the correlation of magnetic susceptibility stratigraphic sections [START_REF] Hladil | The dynamic time-warping approach to comparison of magnetic-susceptibility logs and application to Lower Devonian calciturbidites (Prague Synform, Bohemian Massif)[END_REF], lithostratigraphic correlation [START_REF] Smith | New stratigraphic correlation techniques[END_REF] or for the correlation of well logs [START_REF] Fang | Computer-aided well log correlation[END_REF]. [START_REF] Lisiecki | Application of dynamic programming to the correlation of paleoclimate records[END_REF] propose to correlate paleoclimate record using a dynamic programming method comparable to the DTW algorithm.

Consider two columns a and b presented in Fig. 1, composed of polarity zones noted from a 1 to a 9 and chrons noted from b 1 to b 13 , respectively. The correlation between the two columns can be represented as a path CP (a, b) = [cp 1 , . . . , cp i , . . . , cp k ] in a 2D diagram (Fig. 1) where cp i is an element of the path. In the example presented in Fig. 1 CP (a,b) is noted:

cp 1 cp 2 cp 3 cp 4 cp 5 cp 6 cp 7 cp 8 cp 9 cp 10 cp 11 a 1 a 2 a 3 a 4 a 5 a 6 a 7 ∅ ∅ a 8 a 9 b 1 b 2 b 3 b 4 b 5 b 6 b 7 b 8 b 9 b 10 b 11 b 12 b 13
Four different configurations may arise: (1) two units are well correlated (i.e. diagonal arrow in the diagram noted cp i = (a f , b g ); (2) a polarity zone in column a is missing in column b (i.e. vertical arrow noted cp i = (a f , ∅)); (3) a polarity chron in column b is missing in column a (i.e. a horizontal arrow noted cp i = (∅, b g )); and (4) a polarity zone of a (respectively b) is associated to several polarity chrons of b (respectively a, see the correlation of chrons a 5 in Fig. 1). This last configuration will be addressed in section 1.3.2. A cost is associated to each of these first three configurations. The best correlation between the two columns is the one for which the sum of the individual costs of all successive configurations is minimal. The cost of the correlation CP (a i , b j ) between column a from a 1 to a i and column b from b 1 to b j , C(i, j) is recursively computed as follows:

C a i , b j = min   c(a i , b j ) + C a i-1 , b j-1 g(a i ) + C a i-1 , b j g(b i ) + C a i , b j-1   (1)
where c(a i , b j ) is the cost of a correlation (called a match and noted cp l = (a i , b j )) between two units and g(a i ) and g(b i ) are the costs of an extra unit (called a gap) in column a (noted cp l = (a i , ∅)) and column b (noted cp l = (∅, b j )), respectively. Applying Eq. ( 1) from i = 1 to 9 and j = 1 to 13 yields the minimum cost correlation between columns a and b.

Specification to magnetostratigraphy

In the following section, we present how the DTW algorithm has been improved in order to achieve correlations between a given magnetic polarity column and the GPTS. Let a denote the polarity column established from paleomagnetic sampling and analyses, and b the polarity reference scale.

Cost computation

We first assume that the reference scale b is complete, meaning that each polarity change a i recorded in the polarity column a must be correlated to a polarity change in the reference scale b. Therefore, we set the cost of a gap in the studied polarity column to be equal to the infinity: g(a i ) = ∞. The cost c(a i , b j ) of a match between a polarity chron (of the reference scale b) and a polarity zone (of the polarity column a), is considered only if the two considered units have the same polarity. Moreover, the cost of a match is designed to minimise local variations of the accumulation rate.

Let t(a i ) refer to the thickness, in meters, of the i th polarity zone recorded in the studied polarity column, p(a i ) its polarity, d(b j ) the duration of the j th polarity chron of the reference scale and p(b j ) its polarity. When the correlation between a i and b j is evaluated within the DTW algorithm, the correlation CP (a i-1 , b j-1 ) is known. Let u and v be the last two units correlated before a i and b j :

CP (a i-k , b j-l ) = [cp 1 , . . . , cp k ] (u, v) = cp k-l where u, v = ∅ with l as close as possible to 1 (2)
The match cost between a i and b j , is defined by:

if p(a i ) = p(b j ) c(a i , b j ) = ∞ else c(a i , b j ) = r 1 /r 2 -1 if r 1 /r 2 > 1 c(a i , b j ) = r 2 /r 1 -1 if not.
where,

r 1 = t(a i ) /d(b j ) and, r 2 = t(u) /d(v) (
3) The use of a ratio in the cost function ensures that compared values are dimensionless. This cost function can be computed under the assumption that the accumulation rate is locally stationary in the two considered chrons. Due to discrete sampling, the thickness of a polarity zone in the studied section is not exactly known. Moreover, error is proportionally more important for thin polarity zones than for large ones. Uncertainty on the thickness of small polarity zones is managed by multiplying the match function by a correction factor k(a i ):

k(a i ) = 0.8 × log 1 + min t(a i ), t(a i-1 ) t(a) (4)
where t(a) is the mean thickness of polarity zones on the studied section.

The cost of a gap in the reference scale g(b j ) is computed as the ratio between the duration of the considered chron d(b j ) and the mean duration of chrons of the reference scale d(b). It corresponds to the cost of an unsampled, non-recorded (i.e non-deposition) or eroded polarity zone on the polarity column. It is therefore computed as g(b j ) = d(b j ) /2 d(b) and respects the intuitive idea that non-sampled or non-preserved polarity zones are preferentially the smallest ones.

Gap management Most of the time, when a polarity change has not been sampled or preserved in the sediments, the resulting polarity column may display a polarity interval which actually includes two or more chrons of the reference scale (see correlation between chrons a 5 and b 5 , b 6 and b 7 in Fig. 1). [START_REF] Waterman | The match game: New stratigraphic correlation algorithms[END_REF] introduce the "one to many matching" possibility in the DTW algorithm. In this case, one unit of a column may be correlated to several units of the second column. However, assuming that the magnetic polarity reference scale is complete implies that the "one to many matching" can only occur in one way, between one polarity zone of the studied polarity column and several chrons of the reference scale. Let mg a i , [b j-1 , b j , b j+1 ] denote the cost of a "one to many match" between polarity zone a i of the polarity column and chrons b j-1 , b j and b j+1 of the reference scale. If the polarity of a i , b j-1 and b j + 1 are the same, the "one to many" correlation is considered and

mg a i , [b j-1 , b j , b j+1 ] = c(a i , b j-1,j+1 ) + g(b j )
where b j-1,j+1 is a chron with a duration d b j-1,j+1 = d(b j-1 ) + d(b j+1 ).

(5)

Management of borders A given polarity column observed in a sediment section contains only few units with respect to the reference scale. Therefore, in the DTW algorithm many gaps should arise at the polarity column extremities (Fig. 2). Let's consider a polarity column a with its first polarity zone called a 1 and its last called a n to be correlated to the reference scale noted b. As proposed by [START_REF] Smith | New stratigraphic correlation techniques[END_REF] these extremities issues are managed by setting g(b j ) = 0 in Eq.1 if i = 1 or i = n.

Constraining the correlation using independent dating Some additional dating information derived from biostratigraphy (e.g. Charreau et al. (2009a); [START_REF] Sen | Magnetostratigraphy and biostratigraphy of the neogene deposits of Kastellios Hill (Central Crete, Greece)[END_REF]), cosmonucleid isotopes (e.g. [START_REF] Hu | Cosmogenic burial ages reveal sediment reservoir dynamics along the yellow river, china[END_REF]), palinology (e.g. Dupont-Nivet et al. ( 2008)), etc., may enable to independently constrain the ages of a stratigraphic section. However, such independent age constraints will only be valid at a specific stratigraphic depth. Let a i be the i th polarity zone of the polarity column where one has independently dated the sediment between 15 and 18M a (Fig. 2). The matching cost function is thus transformed to c(a i ,b j ) = ∞ if polarity chron b j is older than 18M a or younger than 15M a. This ensures that polarity zone a i is correlated to reference chrons dated between 15 and 18M a.

Best and n-best correlations Within the DTW algorithm, the cost of a match between a given polarity zone and a chrons (a i and b j ) must be calculated independently from the correlation of the other units (CP (a i-1 , b j-1 )). Our matching function (Eq.3) does not respect this condition. Therefore, to ensure the finding of the minimum cost correlation CP (a, b), the DTW algorithm is modified to output the k correlations [CP (a i , b j )] k between a 1 to a i and b 1 to b j built from modified version of equation 1:

C a i , b j k = min k    c(a i , b j ) + C a i-1 , b j-1 k g(a i ) + C a i-1 , b j k g(b i ) + C a i , b j-1 k    (6) where C a i , b j k are costs associated to the k correla- tions C a i , b j k .
Tests carried out on various polarity columns show that if k is large enough, the optimal correlation is found (for example k > 1, 000 ensures finding the minimum cost correlation between a polarity column of 25 polarity zone and a GPTS of 290 chrons). Furthermore, this new version of the DTW algorithm computes not only the minimum cost correlation but also the n best correlations of a polarity column, offering insight into uncertainty in dating and accumulation rates estimation. As for the minimum cost correlation, a number k of computed correlations should be large enough to ensure finding the n minimum cost ones (k > 20, 000 correlations are output to obtain the n = 10, 000 minimum cost one in example presented in section 1.4).

Tests carried out on synthetic data sets (columns of 7 polarity zones are correlated to a 25 polarity chrons reference scale) where all possible correlations can be computed (125,423 in this case) show that the correlation costs follow an exponential distribution for the better half of the correlation. In a basic magnetostratigraphic problem (a 25 units polarity column to correlate to a reference scale covering the past 150My, i.e. 290 polarity chrons) the number of possible correlations reaches 10 30 which is technically impossible to compute. Therefore only a few percent of the possible correlations can be done which we call the n best correlations. As presented above, k > n correlations should be computed to ensure finding the n best correlations. If k is too low, the costs distribution does not exhibit an exponential distribution, indicating that the n output correlations do not represent the real n minimum cost ones. k can then be increased until the cost histogram honours the exponential function for realisation 1 to n.

Statistical analysis of the n best correlations

The relationship between the cost of a correlation and its rank, considering only the n best (n being low when compared to the total number of possible correlation) could be modelled by:

[CP (a, b)] m 1 k 1 ln k 2 × m + k [CP (a,b)] 0 1 (7)
where, k 1 and k 2 are two constants greater than one, and [CP (a, b)] m is the cost of the m th best correlation. The cost-rank function computed for the study of the Xishuigou section (section 2.2h) is presented in Fig. 4. Two groups of correlations can be identified on the rank-cost function: (i) a first group, corresponding to the first hundred correlations in which the correlation cost increases rapidly with the rank; (ii) a second group, called cost plateau, where the cost gently increases with the rank. Tests carried out on short synthetic polarity columns show that the second group actually includes about half of the entire possible correlations. We define "interesting correlations" as the correlations with a cost close to the minimum one and which therefore can be considered as likely. The first group actually corresponds to this definition of "interesting correlation".

Analysing the n best correlations of a polarity column offers an insight into uncertainties on sediment dating and resulting conclusions (accumulation rate, . . . ). A statistical analysis of uncertainty in a magnetostratigraphic study requires a number n of correlations large enough to ensure that all interesting correlations has been output. A qualitative evidence that all "interesting correlations" have been considered is the presence of a plateau for the costrank function (e.g. as observed in Fig. 4h).

Visualization of the n best correlations A visual and intuitive representation of the possible correlations (10,000 for the magnetostratigraphic column we studied) is a useful tool for the understanding of magnetostratigraphic uncertainties. Here, we propose a visualisation method of the n best correlations based on a representation of the dating weighted density. For each polarity zone a i of the studied polarity column, its mean age is computed on each of the n best correlations [CP (a i , b j )] n and plotted on a depth-time diagram. We propose the contribution w i of each dating CP i (a, b) to be inversely proportional to its cost C i (a, b):

w i = 1 - C 1 (a, b) -C i (a, b) C 1 (a, b) -C n (a, b) (8)
A kernel smoother [START_REF] Wand | Kernel smoothing[END_REF] is then applied to this set of weighted points, resulting in a representation of the alternative correlations where their relative probability is represented by the normalised weighted density map (Fig. 3).

Validation of the proposed approach

The proposed approach is tested on two synthetic polarity columns designed to evaluate how non-stationary accumulation rate, erosion, faulting events and incomplete sampling are managed by the algorithm.

Synthetic column description A first synthetic polarity column (Fig. 3 A) is built from the GTPS using Result analysis Correlation results can be analysed via the minimum cost correlation (dash line in Fig. 3) but also using the statistical distribution of the n best correlations.

Minimum cost correlation

The minimum cost correlation computed on the first synthetic column exactly reproduces the expected correlation. These results show that our correlation method manages non-stationary sediment accumulation rates. The correlation computed on the second synthetic polarity column is consistent with the way we built it. The un-sampled polarity chron (C4r.1r ) is well recognised by the algorithm and results in a one to many correlation in the DTW table. The erosion occurring in chrons C5n.2n and C5r.1r partially well identified. Contrary to a one to many correlation (i.e. chons c of the polarity column correlated with C4r.1r, C5n.1r and C5n.2n), the algorithm outputs an erosion of chron C5n.1r and C5n.2n.

N best correlations

The probability density plot of the 10,000 best correlations of the two synthetic polarity columns are presented in Fig. 3. The quality of the resulting solutions is similar for the two columns. Dating of the central part of the polarity column (from 300m to 1600m and from 250m to 1000m for the first (Fig. 3a) and the second polarity column (Fig. 3b)) does not display major variation around the minimum cost correlation and attests a good confidence on the given ages in this part. At the point where the result of a faulting event has been mimicked in the second synthetic polarity column, variations arise around the minimum cost correlation, suggesting a possible missing polarity zone. Two other variations appear at the top of both polarity columns. In a real case, these variations may lead to alternative interpretations in terms of accumulation rates as presented in section 3.

Application to recent magnetostratigraphic analyses in Central Asia

Central Asia a key area for tectonics, climate and resources in Asia

Central Asia is constituted by several large and very active mountain ranges including the Tianshan, the Qilianshan, the Kunlun Shan, etc. All these ranges dominate the topography over thousands of metres with some summits higher than 7km. While originally formed during a long-lived Paleozoic to Mesozoic history of subduction/collision of several continental blocks (e.g. [START_REF] Windley | Paleozoic accretion and cenozoic redeformation of the Chinese Tien Shan range, Central Asia[END_REF]), their high present topography attests to an intense and younger deformation: all these ranges were uplifted and reactivated under the India-Asia collision collision started 55Ma ago [START_REF] Patriat | India-eurasia collision chronology has implications for crustal shortening and driving mechanism of plates[END_REF], this region has played a key role during the deformation of the Asian continent accommodating up to 40% of the whole convergence between India and Asia. Moreover, intense deformation and tectonics in this area has probably strongly impacted the Asian and global climate, yielding a strong aridification of Central Asia [START_REF] Fluteau | Simulating the evolution of the asian and african monsoons during the past 30 myr using an atmospheric general circulation model[END_REF]. The material eroded and shed from these uplifting mountains fills several large intracontinental basins (Tarim, Qaidam, Junggar) where thick and continuous sedimentary deposits have recorded tectonic and climate changes. Thanks to numerous tectonic wedges and large entrenchments, thick and continuous sedimentary sections are now well exposed at the toes of these mountain ranges, enabling the analyses of sediments for tectonics, climate change, natural resources, textitetc. However, due to their continental origin, these sediments are poorly dated via classic biostratigraphic method, which greatly inhibits our ability to stratigraphically correlate different sections or to place temporal bound on tectonic and climate processes. Therefore during the past decades, magnetostratigraphy has been the technique of choice to better constrain the depositional ages of these series and numerous studies have been published in international journals. However, subjective magnetostratigraphic correlations sometimes yield controversial ages because several authors may propose several correlations for the same section. Consequently, the interpretation and conclusions derived from these controversial ages are open to debate. Two examples of conflicting interpretation will be discussed and treated below.

The Xishuigou section

The Xishuigou section is located in NE Tibet in the Subei area. In this section three authors have proposed three different magnetostratigraphic correlations yielding depositional ages that vary by about 20Ma. [START_REF] Gilder | Oligo-Miocene magnetostratigraphy and rock magnetism of the Xishuigou section, Subei (Gansu Province, western China) and implications for shallow inclinations in central Asia[END_REF] first dated the section between 26Ma and 19Ma (Fig. 4) while later [START_REF] Yin | Tectonic history of the Altyn Tagh fault system in northern Tibet inferred from Cenozoic sedimentation[END_REF] correlated the same section from 33.5 to 27, and [START_REF] Wang | Danghe area (western Gansu, China) biostratigraphy and implications for depositional history and tectonics of northern Tibetan Plateau[END_REF] from 20 to 9Ma. These contrasting results are very puzzling as they lead to conflicting interpretations. It is nevertheless very important to better assess the depositional ages in this section as they bear on tectonics and climate change linked to the NE Tibet uplift. We first output the 10,000 minimum cost correlations considering all magnetic polarity changes identified in Gilder et al.'s first column. Output correlations are presented in Fig. 5. The minimum cost correlation dates the Gilder et al.'s column from reference polarity chron C5n.1n to C5Er which is equivalent to the age proposed by [START_REF] Wang | Danghe area (western Gansu, China) biostratigraphy and implications for depositional history and tectonics of northern Tibetan Plateau[END_REF]. This also holds for the 10 best ones. The first solution equivalent to Gilder et al.'s interpretation is the 1171st with a cost of 10.33 and a weight of 0.27, which do not let it appear in the weighted density plot. Based on the correlation rules described above, this interpretation may thus be interpreted as improbable. The solution proposed by [START_REF] Yin | Tectonic history of the Altyn Tagh fault system in northern Tibet inferred from Cenozoic sedimentation[END_REF] does not appear in the 10,000 minimum cost correlations output. The analyses of the density plot show two alternative solutions of the top of the polarity column. The most probable one indicates an accumulation rate decrease at ca. 13 Ma. When looking at differences between our computed correlations and the interpretation proposed by [START_REF] Wang | Danghe area (western Gansu, China) biostratigraphy and implications for depositional history and tectonics of northern Tibetan Plateau[END_REF] it appears that: (i) we found no faulting event as interpreted by [START_REF] Wang | Danghe area (western Gansu, China) biostratigraphy and implications for depositional history and tectonics of northern Tibetan Plateau[END_REF] at polarity zone g; (ii) in our results, 12 polarity chrons are missing which represent about 50% of the whole reference chrons for the considered period. Given the sampling density, this is unlikely. Such a large discrepancy may suggest that unknown events in the reference scale may have been identified in the Gilder et al's column as suggested by these authors [START_REF] Gilder | Oligo-Miocene magnetostratigraphy and rock magnetism of the Xishuigou section, Subei (Gansu Province, western China) and implications for shallow inclinations in central Asia[END_REF]. Therefore, we performed a second test where short polarity zones b and 11 supported by only two paleomagnetic samples are not considered. Results presented in Fig. 4c show that the minimum cost correlation (with a cost of 6.67) is then exactly the same as the interpretation proposed by [START_REF] Gilder | Oligo-Miocene magnetostratigraphy and rock magnetism of the Xishuigou section, Subei (Gansu Province, western China) and implications for shallow inclinations in central Asia[END_REF]. The 10 best correlations are also in the same range of age. The first output correlation comparable to the interpretation proposed by [START_REF] Wang | Danghe area (western Gansu, China) biostratigraphy and implications for depositional history and tectonics of northern Tibetan Plateau[END_REF] is the 248th one with a cost of 7.94998. The weighted age density plot shows that the two alternative correlations proposed by [START_REF] Gilder | Oligo-Miocene magnetostratigraphy and rock magnetism of the Xishuigou section, Subei (Gansu Province, western China) and implications for shallow inclinations in central Asia[END_REF] are the most 2001), [START_REF] Wang | Danghe area (western Gansu, China) biostratigraphy and implications for depositional history and tectonics of northern Tibetan Plateau[END_REF] and minimum cost correlation computed considering all polarity changes sampled by [START_REF] Gilder | Oligo-Miocene magnetostratigraphy and rock magnetism of the Xishuigou section, Subei (Gansu Province, western China) and implications for shallow inclinations in central Asia[END_REF] and only polarity zones supported by more than 2 samples. (b) and (c) weighted density age versus depth plot of the 10,000 bests correlations computed considering all polarity zones (b) and only polarity zones supported by at least two samples (c). (d) and (e), related cost distribution for the 10,000 best correlations. (e) and (g) related weighted accumulation rate distribution at polarity zones c, 3, 9 and j. (h) cost-rank function for the 10,000 best correlations considering all sampled polarity changes. The proposed model is computed thanks to equation 7 with k 1 = ln(14.7) and k2 = 10 9 .

probable. Then two other groups of probable alternative solutions appears: one group corresponding to the interpretation proposed by [START_REF] Wang | Danghe area (western Gansu, China) biostratigraphy and implications for depositional history and tectonics of northern Tibetan Plateau[END_REF] and one between 2.5 Ma et 10 Ma. The solution proposed by [START_REF] Yin | Tectonic history of the Altyn Tagh fault system in northern Tibet inferred from Cenozoic sedimentation[END_REF] is not reproduced in the computed 10,000 best correlations.

From these results we conclude that, when loose polarity zones b and 11 are removed, the most probable correlation dates the section from chrons C6n to C7Ar*. Although the correlation proposed [START_REF] Wang | Danghe area (western Gansu, China) biostratigraphy and implications for depositional history and tectonics of northern Tibetan Plateau[END_REF] remains possible, no faulting event seems to be justified in the absence of other observation.

The Yaha section

In the Yaha section located in southern Tianshan, two authors also have proposed different correlations with the reference scale. Charreau et al. (2006) established the first magnetostratigraphic column which they correlated to the reference scale from 12.5 to 5.2Ma and, later, thanks to complementary sampling above the first section, from 12.5 to 1.7Ma (Charreau et al., 2009b). The sedimentation rates derived from this correlation show a stepwise acceleration at 11Ma which authors relate to enhanced uplift in the Tianshan mountains. [START_REF] Huang | Magnetostratigraphic study of the Kuche Depression, Tarim Basin, and Cenozoic uplift of the Tian Shan Range, Western China[END_REF] challenge this view and correlate Charreau et al's first bottom section from 9 to 2Ma which yields accelerated sedimentation rates at 7Ma, concluding that the Tianshan underwent intense deformation at this time. Despite some clues indicating that Huang's proposition was probably invalid [START_REF] Charreau | Comment on "Magnetostratigraphic study of the Kuche depression, Tarim Basin, and Cenozoic uplift of the Tian Shan Range, Western China[END_REF] and the fact that the upper section has been sampled since, the debate may remain open. The approach presented here enables us to decipher which correlation is the most probable following rules described in section 1.3 (fig. 5). The analysis of the 10,000 best correlations confirms the interpretation proposed by Charreau et al. (2009b). The minimum cost correlation faithfully reproduces Charreau et al's correlation with a small variation on the top on the polarity column and the 10,000 best correlations all show an accumulation rate acceleration between 10 and 11 Ma.

Certainty and uncertainty on sediments age and accumulation rates

Magnetostratigraphic data and correlation confidence

Evaluating the degree of confidence in a given magnetostratigraphic correlation is critical. It depends on the sampling density and the uncertainty degree of the polarity column correlation to the GPTS. The sampling density and its ability to ensure identification of all the polarity zones recorded in the studied section can be evaluated via s jackknife test [START_REF] Tauxe | A jackknife for magnetostratigraphy[END_REF]. Uncertainty about the correlation mainly originates from: (i) the number of recorded magnetic polarity zones in the studied section and (ii) the distortion of the recorded signal. For example, a succession of three polarity zones may correlate to many parts of the GPTS while a column of fifty polarity zones, recorded in sediment deposited under constant accumulation rate, will yield a unique dating solution. In contrast, non-stationary accumulation rate, faulting event, erosion, etc. may lead to complex polarity columns that are difficult to correlate to the GPTS. A significant problem in magnetostratigraphy is to evaluate whether the information content of a section is sufficient to date rocks with enough confidence. From the set of realisations, a possible criterion for evaluating this aspect could be to look at the spread for the age of each polarity zone (e.g., standard deviation, interquartile range) or the number of modes of the age distribution for each polarity zone. Taking the average over all polarity zones of age standard deviations or interquartiles as a criterion seems a viable way to assess the significance of the correlations. For example, 10,000 realisations of a three-polarity zones section would undoubtedly yield a very large average age spread, meaning that we don't have enough information to deduce anything. Data are in this case considered as noninformative. Conversely, for a longer section, this average dispersion should decrease. An argument against such a procedure is that the method forces alternative solutions to be considered, hence artificially increases the spread when the number of realisations is large. However, as pointed out in Section 1.3.2 the correlation cost reaches a plateau when the number of computed correlations is high. Hence the weight (Eq.8) of unlikely realisations tends to zero. The use of the weighted standard deviation allows computing a factor characterising the capability of a column to provide dating with confidence. The weighted standard deviation for a given polarity zone sd c w is defined by:

sd 2 w = N × N i=1 w i × (x i -xw ) 2 (N -1) × N i=1 w i (9)
where N is the number of realisations, N the number of realisations with a non-zero weight w i , x i the age of the polarity zone and xw the mean age of the polarity zone considering the N realisations. The mean value of the weighted standard deviation s d w can be computed considering all polarity zones of the considered column. Following results are obtained for the Xishuigou section and the Yaha section:

• sd c w = 0.03M a for the Yaha section (Fig. 5b).

• sd c w = 0.72M a for the Subei section considering all sampled reversals (Fig. 4b).

• sd c w = 7.82M a for the Subei section considering only the polarity zones supported by at least 2 samples (Fig. 4d).

A high sd c w suggests that the data available for the magnetostratigraphic correlation are non-informative and thus that additional data coming from other dating methods (e.g., palynology, cosmogenics) should be used to better constrain the correlation.

Sampling bias and accumulation rate uncertainty

From detailed accumulation rates, one may quantify subsidence, sedimentary fluxes (e.g. [START_REF] Métivier | Mass accumulation rates in asia during the cenozoic[END_REF] or reconstruct past erosion rates (e.g. [START_REF] Charreau | Paleo-erosion rates in Central Asia since 9Ma: A transient increase at the onset of Quaternary glaciations?[END_REF]). Evaluating the uncertainty on accumulation rates is therefore critical to better manage these calculations. Even in the Yaha section, though good confidence exists on the correlation to the GPTS, significant uncertainties on instantaneous accumulation rates remain. Fig. 5c displays a weighted density plot of instantaneous accumulation rates, computed as the depth-time density plot (see section 1.3.2) from the 10,000 best correlations of the Yaha polarity column to the GPTS. Figs. 5d, e andf presents the histogram of accumulation rate calculated for polarity zones o, u and 21 respectively, calculated for each of the 10,000 correlations and weighted by each correlation cost. From these diagrams, it turns out that the accumulation rate is most uncertain for the bottom part of the polarity column where multimodal distributions are observed (polarity zones 21 to w). For example, in polarity zone u, it ranges between 0.07 mm/yr and 0.35 mm/yr, with a mode of 0.16 mm/yr. Accumulation rates calculated for polarity zone 21 are comprised between 0.24 mm/yr and 0.3 mm/yr. Above, the uncertainty decreases with accumulation rates varying between 0.5 and 0.64 mm/yr for polarity zone o. The observed instantaneous accumulation rate is dependent on the thickness assigned to a polarity zone. Since the thickness of a given polarity zone is derived from the depth variation of the magnetic direction, discrete sampling should also yield the uncertainty on the exact position of the magnetic polarity change. This uncertainty could be quantified. Let d 1 a and d 1 b denote the depth of two successive samples of opposite magnetic orientation between which lies one of the two limits of the considered polarity zone is and d 2 a and d 2 b the depth of the two samples defining the second limit of the polarity zone. We note ∆

= max(|d 1 a -d 1 b |, |d 2 a -d 2 b |) and δ = min(|d 1 a -d 1 b |, |d 2 a -d 2 b |).
The probability density function of the polarity zone thickness t(x) (Fig. 6) may be taken as:

t(x) =       0 if x < |d 1 b -d 2 a | 0 if x > |d 1 a -d 2 b | x -d 2 a + d 1 b /∆ • δ if x ∈ |d 1 b -d 2 a |, |d 1 b -d 2 a + δ| 1 /∆ if if x ∈ |d 1 b -d 2 a + δ|, |d 1 b -d 2 a + ∆| -x -d 1 a -d 2 b /∆ • δ if x ∈ |d 1 b -d 2 a + ∆|, |d 1 a -d 2 b |      
(10) Errors on polarity zone thickness have a particular impact on the instantaneous accumulation rates computed on thin polarity zones. Examples of its impact can be observed on Fig. 5c where abrupt variation of the accumulation rate on thin polarity zone (c, 4, 6 and 8 for instance) are questionable. In those cases, the thickness probability density function should be used to better eval-uate uncertainties on accumulation rates. Moreover, uncertainty on polarity zone thickness could also impact the dating. The corrective cost computed thanks to equation 4 allows us to handle uncertainty on polarity zone thickness. However it does not take into account the probability density function of polarity zone thickness and thus does not allow propagating it down to the accumulation rate. A way to propagate the polarity zone thickness uncertainty down to the accumulation rate would be: (i)build n polarity columns sampling possible thicknesses of each polarity zone; (ii) correlate the n columns to the GPTS using the proposed approach; (iii) compute weighted age and accumulation rate density plots from the correlations of the n columns. However, this approach, based on a simple sampling of the polarity zone thickness, requires generating and correlating a substantial number of polarity columns and is not applicable in practice. Additional research is thus required to develop an optimised methodology. Other sources of uncertainty in sediment dating and accumulation rate estimation come from the GTPS itself. [START_REF] Agrinier | On the age calibration of the geomagnetic polarity timescale[END_REF] proposed a probabilistic representation of chron age in the GPTS. This representation should be used to incorporate GPTS uncertainty on sediment dating. However, we consider that this uncertainty is low as compared to bias induced by the non-stationary accumulation rate, fault, erosion and sampling strategy and, to a lesser extent, sediment compaction rate.

Conclusion

This paper introduces a new method for correlating a given magnetostratigraphic column to the GPTS with systematic and explicit rules (or cost function, see section 1.3.2). Uncertainties on magnetostratigraphic correlations interpreted by geologists come from characteristics of the sedimentary record (accumulation variability, fault or erosion) but also from the subjectivity of such an interpretation. The proposed numerical approach does not claim to be objective, because correlation rules proposed in section 1 are also subjective. However these correlation rules are explicitly and mathematically expressed and can be discussed, whereas a handmade correlation cannot be evaluated without additional external constraints on the age of rocks (such as fossil, pollen fission track . . . ) . The two main improvements of our method are: (i) its ability to handle, in the stratigraphic record, potential gaps either due to non sampling, sedimentary hiatus (nondeposition and/or erosion) or faulting and non-stationary and variability of sediment accumulation; (ii) the possibility to compute and analyse not only the most probable correlation but also the n best correlations of a section. The main flaw of our technique is that we assume all polarity changes recorded in a studied sedimentary section are known in the GPTS. In particular cases, when short term polarity zone or local sediment re-magnetisation occurs, this limitation may represent obstacle to a good correlation. However, such debatable polarity zones are supported only by a few samples and represent thin polarity zone on the sedimentary section. Since our correlation method manages gaps well, debatable polarity zones could be ignored in the correlation computation and thus automatically managed by adding a gap if they are actually identified proposed method studies. A large view of possible correlations of a given polarity column can be computed in a short amount of time (about 2 hours for the Yaha section when the entire GPTS is considered). A statistical study of sediment deposition age and related accumulation rate is thus possible. Moreover, the method can shed new light on hand made deterministic correlations by providing alternative solutions.

Figure 1 :

 1 Figure 1: Correlation between two polarity columns and associated correlation path in the DTW table displaying gaps (chrons b 10 and b 11 ) and "one to many correlation" (chrons a 5 , b 5 , b 6 and b 7 )

Figure 3 :

 3 Figure 3: Correlation test with a polarity column where accumulation rate is not stationary (A) and where some polarity zones are not sampled or eroded (B). Dashed lines correspond to the minimum cost correlation of each case.

Figure 4 :

 4 Figure 4: Correlation results for the Xishuigou section. (a)Comparison between correlation proposed by[START_REF] Gilder | Oligo-Miocene magnetostratigraphy and rock magnetism of the Xishuigou section, Subei (Gansu Province, western China) and implications for shallow inclinations in central Asia[END_REF],[START_REF] Wang | Danghe area (western Gansu, China) biostratigraphy and implications for depositional history and tectonics of northern Tibetan Plateau[END_REF] and minimum cost correlation computed considering all polarity changes sampled by[START_REF] Gilder | Oligo-Miocene magnetostratigraphy and rock magnetism of the Xishuigou section, Subei (Gansu Province, western China) and implications for shallow inclinations in central Asia[END_REF] and only polarity zones supported by more than 2 samples. (b) and (c) weighted density age versus depth plot of the 10,000 bests correlations computed considering all polarity zones (b) and only polarity zones supported by at least two samples (c). (d) and (e), related cost distribution for the 10,000 best correlations. (e) and (g) related weighted accumulation rate distribution at polarity zones c, 3, 9 and j. (h) cost-rank function for the 10,000 best correlations considering all sampled polarity changes. The proposed model is computed thanks to equation 7 with k 1 = ln(14.7) and k2 = 10 9 .

Figure 5 :

 5 Figure 5: Magnetostratigraphic correlation of the Yaha section. (a) Comparison between correlation proposed by Charreau et al. (2009b), Huang et al. (2006) and the minimum cost correlation. (b) Weighted density age versus depth plot of the 10,000 bests correlations. (c) Weighted density plot of instantaneous accumulation rate. (d), (e) and (f) Histogram of accumulation rates weighted by correlation cost for polarity zones o, u and 21. (g) Cost distribution of the 10,000 bests correlations.

  t(x) of polarity zone i computed from depth of samples d 1 a , d 1 b , d 2 a , and d 2 b and distance separating these samples a, b and c.

  Figure2: DTW table for correlation between a polarity column (vertical) and the polarity reference scale (horizontal). Cells on which matching and gap costs are computed to manage extremities of the polarity column are displayed in dark grey. Cells on which matching and gap costs are computed to constrain correlation with dated point of the polarity column are displayed in light grey.

	24 Gap Cost =0 25	23	22	19 Matching cost = 20 21 for ai	8	18	Age (Ma) 16 17	15	14	11 correlation path 12 13 Possible	8 Point dated between 9 10 15 and 18Ma

Acknowledgements

This work was performed in the frame of the GOCAD research project. The companies and universities members of the GOCAD consortium (www.gocad.org) are acknowledged for their support. Tom Blanchard, Stuart Gilder and two anonymous reviewers are acknowledged for their valuable suggestions for improving this paper. This is CRPG contribution number 2236.


Alexia Chauzy

Arthur Fourmy

Hélène Mirfendereski

William Couet

Sandrine Marchand

Comparison between pharmacokinetic-phramcodynamic modeling results of ceftaroline against methicillin-resistant Staphylococcus aureus after Time Kill and Hollow Fiber experiments

Comparison between pharmacokinetic-phramcodynamic modeling results of ceftaroline against methicillin-resistant

Alexia Chauzy 1 , Arthur Fourmy 1 , Hélène Mirfendereski 1,2 , William Couet 1,2 and Sandrine Marchand 1,2

• Ceftaroline is a time-dependent antibiotic but its pharmacokinetic/pharmacodynamic (PK/PD) target (%T>MIC) in critically-ill patients is still debated.

• A PKPD modeling approach should help to clarify the issue.

• The aim of this study was to evaluate whether simple Time-Kill (TK) experiments would be sufficient, or if more complex and costly Hollow-Fiber (HF) experiments would be necessary for that.

Background

Materials and methods

Results

Conclusion

TK experiments were able to predict the relative efficacy of the various ceftaroline dosing regimen strategies observed during HF experiments, although HF observations did not match perfectly with TK simulations.

1) Static TK curves

• CFU versus time data demonstrated concentration dependent decay and no regrowth (Fig. 1).

Figure 1. Static time-kill curves of ceftaroline against S. aureus Baa1026

Circles correspond to experimental data and lines to predictions from the model.

3) Comparison of TK-derived simulations with HF results

• HF experiments showed:

-Limited transient CFU decay followed by regrowth after once a day dosing (1200mg q24h) (Fig. 3B) -2-log 10 decay after approximately 24h and no regrowth after twice daily dosing (600mg q12h) (Fig. 3C) -2-log 10 decay in less than 6h with no regrowth after an initial 600mg bolus dose followed by continuous infusion of 1200mg over 24h (Fig. 3D)

2) PKPD model

• TK data were best described by a one compartment model for drug-susceptible growing bacteria (Fig. 2).

• Ceftaroline bactericidal effect was modeled as concentration dependent increase in the killing rate according to an Emax model. • These results were relatively, but not fully consistent with TK-derived simulations. Yet bacterial growth in the absence of antibiotic also differed between TK and HF settings (Fig. 3A).

• Accordingly, a single Emax model successfully described both data sets, but parameters values derived from TK and HF experiments were somewhat different (growth rate constant and maximum kill rate constant 65% higher with HF than with TK).

1) Static TK curves

• Methicillin-resistant Staphylococcus aureus strain: Baa 1026

2) PKPD model using NONMEM • Monte-Carlo simulations after 3 various dosing regimens corresponding respectively to 40%, 60% and 100% T>MIC

3) HF infection model (HFIM)

• PK profiles corresponding to Monte-Carlo simulations reproduced in a dynamic HFIM over 4 days • Comparison of TK-derived simulations with HF experimental results

Poster P1337

Figure 2 .

 2 Figure 2. Schematic representation of the final PKPD model used to simulate ceftaroline bactericidal effect on S. aureus Baa1026 Note that a two-compartments model was previously shown to adequately characterize ceftaroline PK in ICU hyperfiltering patients (ECCMID 2020, Abstract 02024) and was combined to the PD model based on TK data. With Cl ceftaroline , ceftaroline clearance; Q, distribution clearance between central and peripheral compartments; k net , apparent growth rate constant of bacteria B; C, ceftaroline concentration; E max , maximum kill rate constant due to ceftaroline; EC 50 , ceftaroline concentration that results in 50% of E max

Figure 3 .

 3 Figure 3. Bacterial counts profiles for various dosing regimens of ceftaroline, corresponding to 0% (A), 40% (B), 60% (C) and 100% (D) T>MIC in a typical patient with creatinine clearance of 200 mL/min, against a MRSA strain. Symbols represent observed bacterial counts in a dynamic HFIM. Straight lines and shaded areas correspond to simulated bacterial counts from the PKPD model based on TK data.