Virginie Hamm
email: v.hamm@brgm.fr

Simon Lopez

IMPACT OF FLUVIAL SEDIMENTARY HETEROGENEITIES ON HEAT TRANSFER AT A GEOTHERMAL DOUBLET SCALE

Deep Triassic aquifers of the Paris sedimentary basin are located in clastic series which are mainly composed of fluvial deposits. They lie between 2000 and 3000 m depth and their temperature can reach 120°C in some areas. There were few attempts (3) at geothermal exploitation in the early 80's but none of them was really successful (exploitation was performed during one year on one site but injection pump failure and well damage lead to abandonment). In the context of the renewal of geothermal exploration and development in the Paris basin several new projects are now considered. As a preliminary study of the specificities of Triassic aquifers, we investigated the impact of sedimentary heterogeneities of such environments on the hydrodynamic and thermal behaviour of a geothermal doublet (1 producer well and 1 injector well). Two orientations of wells were simulated: vertical wells and horizontal wells.The main objectives are i) to study the impact of these heterogeneities and of the type of exploitation (vertical or horizontal wells) on the shape of the cold plume developing around the injector and on the thermal breakthrough time, and ii) to explore a relationship between a heterogeneous model of permeability and a more simplistic model that could potentially describe heat transfer processes at the doublet scale.

For the numerical simulations we used a 3D numerical block obtained with a process-based model of fluvial sedimentary deposits (Flumy model developed by Mines Paristech). This synthetic block is 5 km long, 4 km wide and 50 m thick. It is composed of 200 cells in x and z directions for 160 cells in y direction (6.400.000 cells). For numerical limitation, it was upscaled to 100 cells in x and z directions and 80 cells in y direction (800.000 cells). Doublet simulations were performed considering two cases:the doublet is parallel to the main channels; orthe doublet is perpendicular to the main channels.

The simulation results were then compared to those obtained with a simple hydrogeological model with the heterogeneous model being replaced by a 2 layers model which cumulates sand bodies and clay bodies. For each simulation a passive tracer was added to the cold injected water to compare the geometry of the cooled part of the reservoir to the injected fluid paths.

MODEL DESCRIPTION

In this study we carried out numerical simulations to investigate the impacts of sedimentary heterogeneities on geothermal heat exploitation of fluvial aquifers of the triassic series of the Paris basin using a realistic fluvial reservoir model.

Static reservoir modeling of fluvial deposits

In the oil industry, accurate sedimentary heterogeneities description is key aspect of reservoir characterization as they may have a significant impact on dynamic processes such as water flooding and oil recovery (e.g. [START_REF] Pranter | Analysis and modeling of intermediatescale reservoir heterogeneity based on a fluvial point-bar outcrop analog, Williams Fork Formation, Piceance Basin, Colorado[END_REF]. Static reservoir modeling most often relies on geostatistical techniques to directly reproduce the distribution of the parameters of interest (lithofacies, petrophysical properties…) based on estimations from spatial statistics and on the choice of an appropriate probabilistic model. Even if the classification may seem by way much too simple two family of methods may be distinguished: pixel based models where a continuous function is evaluated on a discrete grid and boolean methods that consist in throwing randomly discrete objects into space. As both types of methods can be used to produce conditional simulation honoring well constraint they are intensively used for practical purposes (e.g. Chilès&Delfiner, 1999).

In this study we used a static 3D numerical block representing the sedimentary architecture produced using a combined stochastic and process-based approach. This alternative method has been developed for modeling fluvial meandering channels at the hydrocarbon-reservoir scale [START_REF] Lopez | Analogue and numerical modeling of sedimentary systems; from understanding to prediction[END_REF][START_REF] Lopez | Modélisation de réservoirs chenalisés méandriformes: approche génétique et stochastique[END_REF]. It is especially suited for exploration phases when data are relatively scarce. In this context, deterministic equations from the hydraulic literature are used to palliate the relative lack of subsurface data. As it mimics sedimentary processes this method produces numerical blocks that are geologically constituent and quite realistic. The different elements have been implemented taking into account physical results and case study reports. The three-dimensional evolution of the channel stems from equations developed in hydraulic studies and proven to generate realistic two dimensional shapes. The model allows for the generation of the different elements of fluvial deposits: point bars, crevasse splays, overbank alluvium, sand and mud plugs [START_REF] Lopez | Analogue and numerical modeling of sedimentary systems; from understanding to prediction[END_REF]. This method has been implemented in the Flumy1 software and has been successfully tested by several oil companies.

Figure1: 3D lithofacies block from the Flumy model. (6 million cells) There are 9 different lithofacies (cf. table 1). Relatively permeable sandy facies are represented in red and yellow (channel lags, sand plugs, point bars) whereas more impervious facies are in shades of green (overbank deposits, alluvial plain) (20x vertical exaggeration)

Description of the 3D heterogeneous block

In this study we used a numerical block produced with the Flumy software which is 5 km long, 4 km large, and 50 m high. The reservoir was generated using a 100 m width channel and simulation parameters were tuned to obtain a "ribbon" type reservoir 2 (named "ribbon" because of permeable sand bodies that are relatively narrow and disconnected). Nine lithofacies are used to represent the sedimentary deposits and fill a regular grid with a 25 m horizontal spacing and a 25 cm vertical spacing which makes 6 million cells. A 3D view of the lithofacies block is shown in figure 1 and table 1 lists the nine lithofacies name.

In order to represent the deep series of the Paris basin petrophysical values (porosity, horizontal and vertical permeabilities) have been assigned to each lithofacies using data from the literature on the Chaunoy oil field [START_REF] Eschard | Combining sequence stratigraphy, geostatistical simulations, and production data for modeling a fluvial reservoir in the Chaunoy field (Triassic, France)[END_REF]. Thermal conductivities have been chosen among usual values for sedimentary deposits and some work could be done to further investigate the impact of the variability of these parameters. As the 6 million cells grid was much too big to perform hydrodynamic and thermal simulation with BRGM sequential simulator Marthe, simple upscaling techniques were used to produce a coarser 800,000 cells grid with a 50 m horizontal spacing and a 0.5 m vertical spacing. The formula used to upscale permeabilities and thermal conductivites in a given direction is the one from Le [START_REF] Loc | Étude de la composition des perméabilités par des méthodes variationnelles[END_REF]cited by Renard et al., 2000[START_REF] Loc | Étude de la composition des perméabilités par des méthodes variationnelles[END_REF] showed that the equivalent permeability for a uniform flow in a given direction has an upper bound which is the harmonic mean of the arithmetic means of the local permeabilities calculated over each slice of cells perpendicular to the given direction and a lower bound which is the arithmetic mean of the harmonic means of the local permeabilities, calculated on each line of cells parallel to the given direction. She proposed to take the geometric average of these bounds as the upscaled value, which writes:

            x h z a y a z a y a x h x K        (1)
Where

x K is the upscaled value in the x direction,

x a

 stands for the arithmetic mean in the x direction and

x h  denotes the harmonic mean in the x direction. Porosity was upscaled using a harmonic mean. Figures 2 and3 show the comparison between the original and upscaled porosity blocks.

MODELING OF A GEOTHERMAL DOUBLET EXPLOITATION

The doublet exploitation consists of a closed loop: fluid is pumped atthe production well and after cooling through a heat exchanger is then fully reinjected in the reservoir. Two critical parameters define the lifetime of ageothermal doublet: the thermal breakthrough time and the extent of the cold plume. The former parameter indicates the time from which the production temperature is decreasing due to a variable mixture with time of the reservoir original fluid and of the cold injected water. The latter parameter indicates the volume of the aquifer which is impacted by the cold injection. These two parameters are crucial for geothermal exploitation, because a premature production temperature cooling will affect the profitability of the exploitation and a good knowledge of the development of cold water bodies in the aquifer is a key aspect of its sustainable exploitation.

Wells location and properties

Two set of wells were considered: vertical wells and horizontal wells. For each set of wells two orientations of the double axis were studied: parallel to main permeable bodies represented by the sand channels (Ox direction) or perpendicular to the channel axis (Oy direction). Wells were simply simulated by assigning a very high permeability value (10 -9 m 2 ) to grid cells along the well path. Vertical wells were located 1 km apart and were considered to cross entirely the simulation block. Horizontal wells were 750 m long and 750 m apart.

They were located on the same horizontal level, in the middle of the block (25 m high in the Oz direction). For every studied case the doublet was centered on the simulation block.

Simulation tool

Dynamic modeling was performed using the MARTHE software package developed by BRGM [START_REF] Thiéry | MARTHE Software.Modeling of Aquifers with a Rectangular grid in Transient state for Hydrodynamic calculations of heads and flows[END_REF]. MARTHE is a finite-difference hydrodynamic package that processes threedimensional flow in both saturated and unsaturated condition and can handle heat transport taken into account the temperature dependence of density and viscosity.

The heat loss by conduction to surrounding impervious formations was not taken into account. This is the reason why temperature production declines are relatively high and would probably be much less for a real case study [START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France[END_REF].

Comparison between wells orientation (vertical or horizontal)

Table 2: Production temperature decline

In this paragraph, we compare the production temperature evolution and the plume extension between vertical and horizontal wells simulations. Two configurations are studied: (1) wells are parallel to the main channels direction and (2) wells are positioned perpendicularly to the main channels direction. The conditions for the numerical simulations are a constant flow rate of 200 m 3 /h and a reinjected fluid temperature of 40 °C. The initial fluid temperature was set to 90°C which is the average temperature of the Trias formation.

Table 2 summarizes the production temperature decline with time after 5, 10, 20 and 30 years of geothermal exploitation (typical duration of a doublet operating permit).

Figure 4shows the temperature evolution in the case of wells positioned in the direction of the channels (x direction) where a temperature decline is observed.

Figure4: Evolution of the production temperature during exploitation for horizontal wells

The comparison between the two modes of exploitation (vertical or horizontal wells) shows a clear difference in the thermal evolution. In the case of wells parallel to the channels, (columns 1&2 of table 2and figure4), the production temperature decline is first observed for horizontal wells. From 6 years of simulation, the temperature drop becomes stronger for vertical wells. The difference between the production temperature curves is growing with time and the production temperature is 6.2 °C colder for vertical wells than for horizontal wells after 30 years of doublet exploitation. The shorter distance between horizontals wells (750 m instead of 1 km for vertical wells) induces shorter path lines between the closest extremities of the horizontal wells but also longer path lines because of the horizontality of the wells (each of 750 m long).Nevertheless, as there are relatively more longer path lines than shorter ones, the production temperature decline will start earlier but will be slower. This shows that horizontal wells are particularly interesting for geothermal purpose because the production temperature evolution is more gradual.In the case of wells perpendicular to the channels,(columns 3&4 of table 2), the numerical simulations for the two wells configurations do not show any perceptible decline after 30 years of exploitation.In the case of vertical wells, the flowrate of 200 m 3 /h couldn't be simulated (numerical divergence) because of too low permeability. Thus for this particular simulation, a flowrate of 100 m 3 /h was used. Those simulations show the importance of fluvial sedimentary heterogeneities on wells position and consequently on the production temperature evolution. Thus, if the heterogeneities are well known, it is a good opportunity (mainly for horizontal wells) to set up wells perpendicular to the principal direction of channels. In this case, the duration of exploitation can be easily extended under the condition that the wells remain in good state enough. For vertical wells, as shown by numerical simulation, there could be a problem of injectivity because the risk is greater not to intersect permeable facieswhen drilling and to be consequently unable to pump at a flow rate compliant with an industrial project economics.

Figures 5to8 show for each wells configuration the geometry of the cold plume around the injector after 10 years of thermal exploitation. As mentionedabove, the extension of the cold plume is not the same in the cases where the wells are perpendicular or parallel to the direction of the main permeable bodies. If we compare the cold plume extension for horizontal wells (same simulation conditions), we observe that unlike the first configuration where the well doublet axis is parallel to Ox (fig. 6), in the second case where the doublet axis is parallel to Oy (fig. 8), the cold plume is developing preferentially in the channel direction, and therefore the cold plume will not reach the production well.

It is worth noting that recently, Ungemachet al.

(2011) stresses the potential benefit of using horizontal wells for the geothermal exploitation of the carbonate Dogger aquifer of the Paris basin.

Comparison with a simulation involving density effect for horizontal wells

In this paragraph, we study the effect of density variation with temperature in the case of horizontal wells parallel to channels. The other simulation conditions remain unchanged.

Table 3, shows the cooling at the production well after a period of 5, 10, 20 and 30 years of geothermal exploitation. Figure 9, shows the two production temperature curves for each simulation.

Table 3:Production temperature decline Figure 9: Comparison of the production temperature for horizontal wells between the simulation with and without density effect.

This simulation shows that density effect induces a delay in the production temperature drop. This trend could be amplified on purpose to delay thermal breakthrough if the location of wells is well designed by setting the injection well deeper than the production well (same depth for this theoretical study).

Comparison between the cold plume geometry and a tracer injection

In this second set of simulations, we compare the cold plume extension to the pathways of an injected tracer.Thus, the 40°C injected fluid is "numerically dyed" using a conservative tracer concentration of 1 g/l. The objective is to assess the impact of the heterogeneities on the displacement of the cold front compared with the displacement of the injected fluid.Figure10, shows the production concentration curves evolution in time for each case (vertical or horizontal wells and parallel or perpendicular to channels)

Figure 10: Production concentration curves

Tracer arrival is observed for every simulation, even for horizontal wells set perpendicular to the channels direction,contrary to what could suggest the production temperature that remained relatively constant in that case. Nevertheless,the tracer breakthrough time is much longer.

Figures11to 14show the isocontours of the invaded areas by both the cold injected water and the injected tracer. The red arrow symbolizes the position of the production well. We observe a clear difference between the volume of aquifer disturbed by the cold injection and the volume invaded by the injected tracer.If we try to compare the respective shapes of the cold plume and the tracer plume, the temperature isocontours are somewhat closeto a distorted cylinder.On the contrary, the tracer isocontours are much more impacted by the heterogeneities with the fluid following preferential paths corresponding to high permeability bodies.

Figures15and 16 show respectively for vertical and horizontal wells, a cross sectionof the reservoir between wells.From top to bottom, we can see the horizontal permeability field (orange: high permeability, blue: low permeability), the temperature and concentration fields after 1 and 10 years of injection (orange: temperature of 90°C or tracer concentration of 1g/l, blue: temperature of 40°C and tracer concentration of 0 g/l).This illustrates that the fluid is preferentially following along high permeability paths. The velocity of the traceris always greater than of cold front velocity [START_REF] De Marsily | Groundwater Hydrology of Engineers[END_REF]. For block porosity ranging from 0,08 to 0,2, the tracer velocity is 3 to 6 times higher than the thermal velocity.Furthermore, as the heat diffusion is 10 3 to 10 5 greater than the tracer diffusion, the thermal front tends to homogenize further away from wells where velocities are smaller, which makes that its shape is much less contrasted than the tracer front.

. In the previous paragraph, we discussed about the fact that heterogeneities can lead to channeling along specific paths and lead to short first tracer arrival times. By contrast, the thermal front seems to be less impacted by channeling because the fluid is heated in contact with the rockas it is moving. Here we propose a simple method to replace the heterogeneous model, in the case of wells set in the channels direction, by an alternative homogeneous model that could describe well enough heat transfer processes at this scale.

Methodology used

The method used consists in replacing the threedimensional heterogeneous block characterized by 100,0.5 m thick subdivisions in the z direction, by a two-layer model with one layer cumulating the high permeability zones (net pay) and the second layer cumulating low permeability zones (interstrata).Figure 17shows the distribution of permeability in each block direction (x, y and z) for the heterogeneous model. It shows that more than 60% of cells have horizontal permeability smaller than 2.10 -15 m 2 and more than 20% of cells have horizontal permeability higher than 2.10 -13 m 2 .

In order to calculate the thickness of the potentially productive reservoir (net pay), we defined a threshold of permeability ranging from 1.10 -14 to 4.10 -13 m 2 and we calculated the net pay and transmissivity field for each cells (i,j) using the following equations:

ℎ 𝑖, 𝑗 = 1 𝑘 𝐾 ℎ 𝑖,𝑗 ,𝑘 ≥𝐾𝑚𝑖𝑛 ∆𝑧 (2) 𝑇 ℎ 𝑖, 𝑗 = 𝐾 ℎ 𝑖, 𝑗, 𝑘 𝑘 𝐾 ℎ 𝑖,𝑗 ,𝑘 ≥𝐾 𝑚𝑖𝑛 . ∆𝑧 (3)
Where:

 ℎ 𝑖, 𝑗 is the productive reservoir thickness or net pay (m)  𝐾 ℎ 𝑖, 𝑗, 𝑘 is the horizontal permeability of cell (i,j,k) (m 2 )  𝐾 𝑚𝑖𝑛 is the permeability threshold (m 2 )  ∆𝑧 is the cell thickness in z direction (m)  𝑇 ℎ 𝑖, 𝑗 is the horizontal transmissivity of cell (i,j) (m 3 ) Figure 19 shows the 3D distribution of cells that are above this threshold. This figure illustrates the intricate patterns of connectivity paths between permeable bodies. Figures 20 and21 show respectively the two-dimensional distribution of the net pay and the transmissivity field obtained by vertically averaging the previous 3D distribution.

Figure19: Cells of the coarse grid that have an horizontal permeability above the 5.10 -14 m 2 threshold (K h ≥ K min =5.10 -14 m 2 )(40x vertical exaggeration) The thickness of the second layer representing low permeable areas,that are barriers to flow but that will transmit heat, is obtained as the difference between the total thickness (50 m) and the productive thickness.To assign a permeability to this layer we have chosen to retain the vertical permeability of the original block considering that the exchanges with the main flow that is located in the productive layer will occur perpendicularly to this layer. Moreover, as the permability values involved are relatively low this choice has no significant impact on the heat transfer between the two layers that is dominated by heat conduction. Consequently we have chosen the term of "vertical transmissivity" to note that the transmissivity of this layer was obtained using the vertical permeability. This writes using the same form as equation ( 3):

𝑇 𝑧 𝑖, 𝑗 = 𝐾𝑧 𝑖, 𝑗, 𝑘 𝑘 𝐾 ℎ 𝑖,𝑗 ,𝑘 <𝐾𝑚𝑖𝑛 . ∆𝑧 (4)

The sum being taken on the cells below the K min permeability threshold (complementary set of sedimentary body depicted in figure 19). The notations are identical to equation (3)with the vertical permeabilityK z replacing the horizontal permeability K h and T z . replacing T h .However, as the permeability involved are one order lower than the productive layer "real transmissivity" they will not have a great impact on the flow behavior. Figure22shows the interstratatransmissivity field.

Comparison between the simplified and heterogeneous models

Temperature and tracer concentration curves

Here we compare the production temperature and the tracer restitution between the simple two-layer model and the heterogeneous block simulations. Table 4summarizes the production fluid temperature decline after 5, 10, 20 or 30 years of exploitation for each set of simulations. The results show that the production temperature simulated by the two-layer model is rather pessimistic predicting a premature cooling at the production well. On the other hand, the first arrival time of the injected tracer is slightly later than in the heterogeneous simulations. These simulations lead to the following observations. By homogenizing the permeability on the vertical scale (cumulation of vertical cells distributed into two lithofacies), we create artificial lateral connectivity around and between wells which makes that the fluid will move through a larger volume of the domain and will reach the production well with some delay with the respect to what occurs with the connectivity of the heterogeneous model. After the tracer breakthrough time, the tracer concentration at the production well is increasing more rapidly. The production flow is consequently composed of a greater proportion of the injected fluid what will also induce premature temperature decline in this case. Our conclusion is that with such a rough homogeneous two-layer model we cannot reproduce well enough the production temperature evolution observed using the heterogeneous model.

Extension of the cold plume

We also compare the form of the cold plume between the two models. Figures 27 and28show the temperature isocontours respectively for a vertical and horizontal cold injection. We observe that in the homogeneous model, the cold front is growing mainly in the horizontal direction unlike the heterogeneous model where the cold front is developing also in the vertical direction. The shape of the cold front is like a "ship hull" with isocontours gradually warmer and closed on the bottom whereas it looks like an open distorted cylinder in the heterogeneous model.

CONCLUSIONS AND PERSPECTIVES

In this paper we analyzed the impact of fluvial heterogeneities on the flow and heat transport induced by geothermal exploitation at the doublet scale. We used a 3D numerical block representing realistic fluvial sedimentary heterogeneities as they are supposed to be found in the Triassic clastic series of the Paris basin. The methodology of our work consisted, firstly,to quantify the impact of these heterogeneities according to the mode of operation (vertical or horizontal wells), thento compare the development of the cold front with a conservative tracer plume, and finally to test if the vertical heterogeneities could be replaced by a simpler model.

The simulations show that in such environment, horizontal wells are particularly well adapted with a more gradual temperature decline at the production well. Moreover, when drilling horizontal wells perpendicular to the high conductive paths, the production fluid temperature is still not affected by the cold injection after 30 years of geothermal exploitation. If these results are confirmed at an experimental stage, it would be a promising breakthrough for further development in geothermal well design.

Furthermore, compared to a tracer injection which "marks" the flow paths, the volume of the reservoir impacted by the cold injection is smaller and less impacted by sedimentary heterogeneities. Following these simulation results, we tested a two-layer model that cumulates high permeable cells in the heterogeneous model in one layer (net pay) and low permeable cells in a second layer (interstrata). The horizontal heterogeneous structure was unchanged. This attempt of simplifying the vertical structure of the heterogeneous 3D block failed in fitting temperature decline as observed in the heterogeneous model. Further study need to be done to better understand the impact of connectivity on heat transport, and to develop a method that could substitute complex sedimentary heterogeneous block by a simpler approach giving similar results.

Figure5:

  Figure5: Cold plume extension after 10 years of geothermal exploitation (vertical wells positioned in parallel to channels, x direction).

Figure 6 :

 6 Figure 6: Cold plume extension after 10 years of geothermal exploitation (horizontal wells positioned in parallel to channels, x direction)

Figure 8 :

 8 Figure 8: Cold plume extension after 10 years of geothermal exploitation (horizontal wells positioned perpendicularly to channels, y direction).

Figure 11 :

 11 Figure 11: Temperature and tracer isocontours after 10 years of exploitation (vertical wells oriented according x direction).

Figure 12 :

 12 Figure 12: Temperature and tracer isocontours after 10 years of exploitation (horizontal wells oriented according x direction).

Figure 13 :

 13 Figure 13: Temperature and tracer isocontours after 10 years of exploitation (vertical wells oriented according y direction).

Figure 14 :

 14 Figure 14: Temperature and tracer isocontours after 10 years of exploitation (horizontal wells oriented according y direction).

Figure 15 :

 15 Figure 15: From top to bottom: horizontal permeability field and position of vertical wells (red lines), temperature and tracer fields after 1 year (left) and 10 years (right)

Figure 17 :

 17 Figure 17: Distribution of permeability (m 2 ) in the heterogeneous block

Figure 20 :

 20 Figure 20: Productive thickness in the homogeneous model (K h ≥ K min =5.10 -14 m 2 )

Figure 22 :

 22 Figure 22: Interstratatransmissivity (m3)

Table 4 :

 4 Figures 25and26show respectively the production temperature and tracer concentration in the case of vertical or horizontal wells and for each model structure.

Figure 25 :

 25 Figure 25: Comparison of the production temperature (solid lines) and tracer concentration (broken lines) between the heterogeneous model and the two-layer homogeneous model in the case of vertical wells

Figure 27 :

 27 Figure 27: Comparison of the cold plume extension between the homogeneous model on the top and the heterogeneous model on the bottom. Case of vertical wells.

Figure 28 :

 28 Figure 28: Comparison of the cold plume extension between the homogeneous model on the top and the heterogeneous model on the bottom. Case of horizontal wells.

Table 1

 1 lists petrophysical properties values that have been assigned to each lithofacies(rocktyping).

Table 1 :

 1 Lithofacies and associated petrophysical values2 A similar numerical block may be downloaded from the Flumywebsite (cf. note 1).

	lithofacies	porosity	horizontal permeability	vertical permeability	thermal conductivity
		(%)	(mD)	(mD)	(W/m/°C)
	channel lag	15	400	400	4.5
	sand plug	10	50	1	4.5
	point bar	8	400	5	4.5
	crevasse splay I	12	80	80	2.5
	crevasse splay II channels	8	50	1	3
	crevasse splay II	20	1	0.1	2.5
	mud plug	20	0.1	0.01	1.5
	overbank alluvium	20	0.2	0.02	1.5
	levee	12	1	0.1	2

Wells orientation Parallel to channels Perpendicular to channels Wells position Wells orientation

		vertical	horizontal vertical horizontal
	ΔT (°C)
	after 5	0.78	1.12	0	0
	years
	ΔT (°C)
	after 10	6.86	4.17	0	0
	years
	ΔT (°C)
	after 20	15.19	9.51	0.01	0
	years
	ΔT (°C)
	after 30	19	12.77	0.16	0.02
	years

Type of model structure heterogeneous model homogeneous model Wells orientation Wells orientation that

  Thus those simulations show vertical heterogeneities have an impact on the development of the cold plume though. The two-layer model is not precise enough to reproduce the cold plume geometry.

		vertical	horizontal vertical horizontal
	ΔT (°C)
	after 5	0.78	1.12	4.03	4.67
	years
	ΔT (°C)
	after 10	6.86	4.17	12.49	9.75
	years
	ΔT (°C)
	after 20	15.19	9.51	16.13	12.88
	years
	ΔT (°C)
	after 30	19	12.77	18.32	14.62
	years

More information about the Flumy software and free evaluation versions may be downloaded at : http://www.geosciences.minesparistech.fr/en/organization/presentation-of-thegroup-2/main-projects/flumy-turmy

ACKNOWLEDGMENTS

This study is part of the CLASTIQ-2 project cofunded by BRGM and the French Environment and Energy Management Agency (ADEME). We thank the ADEME for its financial support that permitted to achieve this work.

The authors are also grateful to Mines Paristech who courteously gave them the 3D numerical block produced with the Flumy software.


Guillaume Ploussard

Alexandre De

La Taille

Younes Bayoud

Xavier Durand

Stéphane Terry

Evanguelos Xylinas

Yves Allory

Francis Vacherot

Claude-Clément Abbou

Laurent Salomon

The risk of upstaged disease increases with body mass index in low-risk prostate cancer patients eligible for active surveillance

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


Monica Dapiaggi
email: monica.dapiaggi@unimi.it

Filippo Maglia

Ilenia Tredici

Beatrice Maroni

Giulio Borghini

Umberto Anselmi Tamburini

Filippo Maglia

Complex thermal evolution of size-stabilized tetragonal zirconia

Keywords: Nanopowders, Zirconia, Phase transition

 A c c e p t e d ma n u s c r i p t 3

Introduction

Zirconium oxide undergoes a transition from its monoclinic to its tetragonal polymorph at around 1175°C, and from its tetragonal to its cubic polymorph at 2370°C. The low temperature monoclinic polymorph, the mineral baddeleyite, has very limited practical use, while the two high-temperature polymorphs are employed in a variety of technical applications such as, for instance, solid electrolytes in solid oxide fuel cells and sensors [START_REF] Stevens | An introduction to zirconia[END_REF]. The stabilization of the tetragonal and cubic phases is traditionally realized by doping with lower valence rare earth oxides. Y 2 O 3 is most frequently used in the zirconia solid solutions of practical relevance; a minimum amount of Y 2 O 3 above 2%mol. and 8%mol. is required for the stabilization of the tetragonal and cubic polymorphs respectively. It has been recently shown that tetragonal ZrO 2 can also be size-stabilized at room temperature, without doping, provided that the grain size is smaller than a critical value of around 10-40 nm depending on the synthesis condition adopted [START_REF] Garvie | [END_REF][3][4][5][6][7][8][9].

Since its discovery by Garvie [START_REF] Garvie | [END_REF][3], size stabilized zirconia has been the object of numerous studies devoted in particular to the investigation of the mechanism responsible for the inversion of phase stability. Among the possible causes, difference in surface energy between the tetragonal and monoclinic polymorphs is the explanation that received the largest consensus, although the presence of water vapor in crystallite growth, of intrinsic defects, and of lattice and intergranular strain are also believed to play a role [9].

The structural and kinetic aspects concerning the stability of the size stabilized tetragonal zirconia and its transition to the monoclinic polymorph has been investigated although in less detail [10][11][12][13][14]. This transition occurs upon annealing at high Many of the investigations concerning this point were performed at room temperature examining samples previously annealed at higher temperatures. The details of the transition mechanism and, in particular, its kinetics were consequently overlooked.

Recent studies of pure zirconia [12][13] and on zirconia doped with low amounts of Sc 2 O 3 [14], revealed as the transition towards the monoclinic phase is kinetically more complex than previously thought. It is completed, in fact, only during the cooling stage, with a transition temperature that depends largely on the thermal history of the sample.

Aim of the present study was to investigate, through in-situ X-ray diffraction, the structural and microstructural evolution of Y 2 O 3 -ZrO 2 nanocrystalline solid solutions, with YO 1.5 ranging between 0 and 1%mol., i.e. below the thermodynamic limit for the stabilization of the tetragonal phase, in order to elucidate the role played by the grain size, doping and strain content on the tetragonal to monoclinic transition.

Experimental

ZrO 2 nanopowders with a dopant (YO 1.5 ) content ranging between 0 and 1% at.Y (i.e. from pure ZrO 2 to Zr 0.99 Y 0.01 O 2-δ ) were obtained by a modified Pechini method [15]. An aqueous solution containing the proper amounts of Zirconyl (IV) nitrate and yttrium nitrate was mixed with a solution of citric acid in 1:1 (total metal vs citric acid) molar ratio. The solution was then stirred on a hot plate at 80°C, increasing its viscosity. After a time period that depends on the amount of solution employed, the viscous mass finally turns into a colorless transparent glass. Further heating provided a white solid, which was subsequently ground in a mortar and calcined in a furnace, at 500°C for 1 h.

A c c e p t e d m a n u s c r i p t 5

In-situ high temperature diffraction (an AHT PAP1600 hot chamber was used) patterns were collected in static air, between room temperature and 1200°C. Data were collected during both heating and cooling stages every 100°C, the heating rate between each measuring temperature was of 20°C per minute.

The hot chamber is attached to a Philips X'Pert θ-θ powder diffractometer; data collections were performed on a limited 2θ range, in order to be able to get the data fast enough to follow the phase transition easily. The range was from 15 to 50° 2θ, with a step size of 0.03° 2θ and a counting time of 1 sec per step. The wavelength used is Cu Kα1 (1.5406 Å), with an incident slit of ½ °, an antiscatter of ½ ° on the diffracted beam, and a receiving slit of 0.6 mm, in order to maximize the peak intensity; resolution, on the other hand, was not of primarily concern in our work due to considerable line broadening related to the nano-size of our powder. The sample holder is made of dense alumina, covered with a thin foil of platinum (Goodfellows, 99.9%), as one of the peaks of α-alumina tend to interfere with the main peak of tetragonal ZrO 2 .

The data were analysed with the Rietveld method to obtain the weight percent of the zirconia polymorphs, their cell parameters, and the profile shape analysis. The sample holder peaks were fitted with the Le Bail method [16], in order to remove their influence on the quantitative analysis. The software used was [START_REF] Larson | General Structure Analysis System (GSAS)[END_REF][START_REF] Toby | [END_REF].

As the small number of peaks collected did not allow any sensible microstructural characterization of the samples during heating, the crystallite size and the RMS microstrain were evaluated from the corresponding Lorentzian components of the peak shape of the modified pseudo-Voigt (function 2 in GSAS), by keeping the Gaussian component (which approximately correspond to the instrumental contribution) fixed at all temperatures. In more details, the crystallite size can be estimated from the LX

A c c e p t e d m a n u s c r i p t

6 parameter (D=LX/cos(θ)) and the RMS microstrain from the LY parameter (S=LY tan(θ)). This is clearly not the ideal way to operate, especially due to the fact that some strain component varies with the diffraction angle in a way which is similar to the Caglioti/Gaussian component.

Results and discussion

Three different samples were analyzed with dopant (YO oxides with crystallite and grain sizes in the nanometric range; the crystallite/grain size of the as-synthesized material is generally dependent on the annealing treatment. The latter is necessary to promote the crystallization of the oxide phase, initially obtained in amorphous form, and to remove the organic residuals. On the basis of DTA-DSC and XRD analyses the starting material has been in this case annealed at 500°C for 1 h.

Structural and microstructural characterization of the starting powders was performed, using the software MAUD [19], using high quality datasets collected at ESRF (ID31)

during the experiment HS3654. These high energy patterns, with a wavelength of about 0.4 Å, allowed indeed a much better accuracy, in particular for what concerns the microstructural parameters evaluation, since the reciprocal lattice could be explored in much higher detail.

The results of this analysis are shown in Table 1. All the samples were free of monoclinic polymorph indicating that the synthetic procedure adopted was successful in the stabilization of the high temperature tetragonal polymorph. The cell volume As the temperature is further increased the amount of the monoclinic phase increases up to above 50% at 1100°C; at 1200°C the tetragonal polymorph returns to be the only phase present as at this temperature the tetragonal phase becomes thermodynamically stable. Similar qualitative trends were observed also for the other two compositions. Upon cooling, as the temperature is reduced below 700-800°C, the monoclinic fraction reaches values above 95%; finally, at room temperature no residual tetragonal ZrO 2 could be observed.

The data collected at each temperature were analyzed using the Rietveld method, in order to obtain a quantitative evaluation of the phase composition during heating and cooling. An example of a Rietveld fit is reported in Figure 2. Figure 3 shows the quantitative analysis at each temperature, in terms of the monoclinic polymorph wt.%, The transition is clearly different during the heating and the cooling stages. Upon cooling the transformation kinetic is completed within a narrower temperature interval (< than 100°C) and the onset temperature for the tetragonal to monoclinic transition decreases going from Zr_0 to Zr_1. This indicates, as expected, a larger stability field for the tetragonal phase increasing the doping content. This behavior, that follows the well known behavior for "conventional" microcrystalline sample, was expected as, at the end of the heating cycle, grain coarsening goes well above the limit for any size stabilizing effect.

The situation is indeed more complex during the heating cycle. Not only, in fact, the kinetics is much more shallow, with only about 50% of the tetragonal phase that transforms into the monoclinic one within a 600-700°C temperature interval, but also the temperature for the transformation onset does not follow the doping content, being the Zr_05 sample the first to initiate the transition at ~400°C and the Zr_1 the last one at a temperature of ~600°C.

During the heating cycle, the transition is more likely to occur under microstructural control. Size stabilization is in fact expected to play the leading role, since the doping contents employed are well below the minimum limit for the stabilization of the tetragonal phase.

A qualitative insight on the evolution of the microstructure during the heating cycle can be gained from Figure 1 (upper part) by observing the peak shape evolution of the tetragonal polymorph. As expected, a shrinkage in the width of the peaks of the tetragonal phase is observed as the temperature is increased. at an almost identical temperature of around 400°C and than proceeds with a slope that is somewhat less pronounced for the Zr_1 sample than for the other two.

In the vast majority of the literature, the width of the main diffraction peak is used for a straight estimation of the crystallite size using the Scherrer formula. The temperature evolution of the phase composition of nanocrystalline pure zirconia has been therefore justified, in most cases, solely on the basis of the grain coarsening.

If applied to our results, where for the first time different low doping contents are taken into account, such interpretative model is clearly not sufficient to justify the experimental evidences. The negligible difference in the initial grain size for the three compositions (Table 1) and the results of figure 4 cannot justify the clear differences in the onset temperature for the phase transition.

To clarify this point an evaluation of the RMS microstrain has been performed for all the samples at the different annealing temperatures (heating cycle). The separation of crystallite size and of the RMS microstrain contribution to the Lorentzian peaks was, differently from the case of the high quality synchrotron data of Table 1, possible in a more approximate way, as we are using a smaller 2theta range for our analyses. For that reason it was chosen to show the RMS microstrain results only as variations with respect to the more precise starting values of Table 1.

The results are shown in Figure 5: the relationship between the RMS microstrain and the onset temperature for the phase transition (see Figure 3) is evident. The Zr_05, characterized by the initial highest value of microstrain, starts to release the microstrain at the lowest temperature which also corresponds to the beginning of the monoclinic A c c e p t e d m a n u s c r i p t 10 phase formation. Subsequently, at progressively higher temperature, the same behavior is followed by the Zr0 and finally the Zr1 sample, suggesting that microstrain release is the process responsible for the initiation of the phase transformation at low temperature.

Conclusions

The thermal behavior of zirconia nanopowders with low yttria contents was investigated by in-situ high temperature XRD. All starting materials were characterized by a doping content well below the thermodynamic limit for the stabilization of the tetragonal phase.

The tetragonal phase was on the other hand stabilized thanks to the reduced dimension of the grain size. Room temperature XRD measurements showed as all sample had a crystallite size around 10 nm while the RMS microstrain content varied for the different doping contents.

Upon heating the size-stabilized tetragonal phase transforms partially into the monoclinic one. The onset temperature for this transformation is strongly dependent on the RMS microstrain content and occurs at the temperature where the material starts releasing its stress content. The transition to the monoclinc phase cannot be completed upon heating before entering, at 1200°C the stability range of the tetragonal form. On the other hand, once the grain size is increased as an effect of annealing, the materials rapidly transforms into the monoclinic phase cooling down from from 1200°C to room temperature. The onset temperature of the transition is in this case dependent on the doping content.

A c c e p t e d m a n u s c r i p t 4 temperature

 4 and is believed to be driven by the grain growth of the tetragonal ZrO 2 .

A c c e p t e d m a n u s c r i p t 7 increases

 7 with the doping content coherently with the larger size of the Y 3+ ion with respect to Zr 4+ . The Zr_0 and Zr_05 showed comparable values of crystallite size (around 12 nm) while Zr_1 showed a slightly smaller value. The grain size is for all the prepared samples below the previously reported limit for the complete stabilization of the tetragonal phase in the pure ZrO 2 . Although there is some scatter in the literature about this limit, values between 10 and 40 nm are generally reported. No data are available in the literature for the Zr_05 and Zr_1 compositions as no previous reports on these compositions are at hand. The RMS microstrain variation was on the other hand more significative; the RMS microstrain value for the Zr_05 sample was approximately twice that of the Zr_1.

Figure 1

 1 Figure1shows the qualitative evolution of the HT-XRD pattern collected during the

A c c e p t e d m a n u s c r i p t 8 during

 8 heating (continuous lines) and cooling (dotted lines), that allows a direct comparison among the three samples.

Figure 4 shows the trend A c c e p t e d m a n u s c r i p t 9 of

 9 the Lorentzian peak width at 29.9°, as a function of the annealing temperature, for the three compositions. For all YO 1.5 contents the shrinkage of the diffraction peak begins

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t 12 6Figure 1 .

 121 Figure 1. In-situ high temperature XRD patterns for the Zr_0 sample. Upper part:

Figure 2 .

 2 Figure 2. Example of a Rietveld fit of the Zr_05 sample at 800°C.

Figure 3 .

 3 Figure 3. Quantitative analyses of monoclinic polymorph as a function of temperature

Figure 4 .

 4 Figure 4. Thermal evolution of the peak shape and of the monoclinic polymorph during

Figure 5 .A c c e p t e d m a n u s c r i p t

 5 Figure 5. Thermal evolution of the RMS during the heating stage.

Figure 1 .

 1 Figures

A c c e p t e d m a n u s c r i p t 15 Figure 2 :

 152 Figure 2: Example of a Rietveld fit of the Zr_05 sample at 800°C.

A c c e p t e d m a n u s c r i p t 16 TemperatureFigure 3 .A c c e p t e d m a n u s c r i p t 17 TemperatureFigure 4 .

 163174 Figure 3. Quantitative analyses of monoclinic polymorph as a function of temperature

A c c e p t e d m a n u s c r i p t 18 TemperatureFigure 5 .

 185 Figure 5. Thermal evolution of the RMS during the heating stage.

A c c e p t e d m a n u s c r i p t

13 Tables Table 1 -

 Tables1 Structural and microstructural parameters of the as-prepared powders.

	Sample	<D>v (nm)	RMS strain (10 -3 ) Cell volume (Ǻ 3 )
	Zr0	12.2 (2)	4.18 (1)	67.006 (4)
	Zr0.5	12.7 (1)	5.97 (4)	67.037 (4)
	Zr1	9.81 (3)	2.91 (6)	67.043 (6)

Acknowledgements

The authors would like to thank ESRF for providing beam time (experiment HS3654), and the staff at ID31, in particular Michela Brunelli and Andy Fitch for their help during data collection. Fondazione Cariplo is thanked for funding.


Student Member, IEEE Miguel Altuve
email: maltuve@usb.ve

Guy Carrault
email: guy.carrault@univ-rennes1.fr

Alain Beuchée
email: alain.beuchee@chu-rennes.fr

Patrick Pladys
email: patrick.pladys@chu-rennes.fr

Alfredo I Hernández
email: alfredo.hernandez@univ-rennes1.fr

On-line apnea-bradycardia detection using hidden semi-Markov models

In this work, we propose a detection method that exploits not only the instantaneous values, but also the intrinsic dynamics of the RR series, for the detection of apneabradycardia episodes in preterm infants. A hidden semi-Markov model is proposed to represent and characterize the temporal evolution of observed RR series and different pre-processing methods of these series are investigated. This approach is quantitatively evaluated through synthetic and real signals, the latter being acquired in neonatal intensive care units (NICU). Compared to two conventional detectors used in NICU our best detector shows an improvement of around 13% in sensitivity and 7% in specificity. Furthermore, a reduced detection delay of approximately 3 seconds is obtained with respect to conventional detectors.

I. INTRODUCTION

Apnea-bradycardia (AB) episodes are common on preterm infants. These episodes are defined as a respiratory pause, accompanied with a fall in heart rate, and are associated with short-term morbi-mortality and neurological impairment during childhood [START_REF] Janvier | Apnea is associated with neurodevelopmental impairment in very low birth weight infants[END_REF], [START_REF] Pillekamp | Factors influencing apnea and bradycardia of prematurityimplications for neurodevelopment[END_REF]. AB episodes may be aggravated in the presence of infection, hypoxia or intracranial pathologies. Preterm infants suffering from frequent AB episodes are monitored in neonatal intensive care units (NICU), where the detection of AB triggers a set of nursing therapeutic actions.

In current NICU monitors, AB episodes are detected by processing a single variable extracted from the electrocardiogram (ECG): the cardiac cycle length (RR) [START_REF] Poets | The relationship between bradycardia, apnea, and hypoxemia in preterm infants[END_REF], [START_REF] Cruz | Algorithm fusion for the early detection of apnea-bradycardia in preterm infants[END_REF], [START_REF] Portet | Evaluation of on-line bradycardia boundary detectors from neonatal clinical data[END_REF]. However, a recent work has shown that other features extracted from the ECG, such as R-wave amplitude and QRS complex duration, also show significant changes before and after the episode [START_REF] Altuve | Multivariate ecg analysis for apnoea-bradycardia detection and characterisation in preterm infants[END_REF]. Typically, AB episodes are detected when the cardiac cycle length is higher than a predefined fixed threshold (usually 600 ms) or when it is higher than a relative threshold (usually 33% of the base rhythm), for more than 4 s [START_REF] Poets | The relationship between bradycardia, apnea, and hypoxemia in preterm infants[END_REF].

The aim of this work is to propose an AB detector that exploits not only the instantaneous values, but also the dynamics of the RR series. An approach based on a hidden semi-Markov model (HSMM) is thus presented, with different pre-processing methods. The rest of this paper is organized as follows: the next section introduces the phases of learning and evaluation of HSMM to perform the event detection, as well as the pre-processing methods applied (quantification and integration of delayed version of the observation) and the evaluation methodology. The evaluation of the proposed detector on simulated and real signals is described in Section III. Discussion of the results obtained are exposed in Section IV. Finally, the conclusions and future works are outlined in the last section.

II. METHODS

A. Hidden Semi-Markov Models

An HSMM is similar to a classic hidden Markov model (HMM), but the main difference is that the unobserved process is semi-Markov in the sense that a change to a future hidden state depends on both the current state and the time spent on this state. These models have already been applied successfully in the bioengineering field [START_REF] Faisan | Hidden semi-Markov event sequence models: Application to brain functional MRI sequence analysis[END_REF], [START_REF] Dumont | Detection of myocardial ischemia with hidden Semi-Markovian models[END_REF].

An HSMM is characterized by a number of M states and the set of parameters λ {a ij , b i , π i , p i }, where a ij is the transition probability between states i and j (a ii = 0), b i is the probability of emission of observations, π i is the probability of the initial state and p i is the probability of duration for state i. In this work, the Bayesian information criterion (BIC) [START_REF] Schwarz | Estimating the dimension of a model[END_REF] is used to calculate M . Additionally, b i and p i are represented by a Gaussian distribution: b i ( µ, Σ) and p i (µ d , σ d ), where µ and Σ correspond to the center and covariance matrix of observations and µ d and σ d are the mean and the standard deviation of the duration of each state. In order to avoid negative durations, Gaussian distributions for p i (µ d , σ d ) are always positive (truncated in zero).

The main interest in the application of an HSMM to biomedical series processing is that it can better represent the temporal evolution of a variable through parameters p i . Indeed, the Gaussian law used in the proposed HSMM for p i provides a better representation of the time spent in a given physiological state than the implicit geometric law of HMM, in particular when the duration on a given state is relatively long. In a previous work, we have shown the advantages of the HSMM approach, with respect to classic HMM, for the detection of cardiac ischemia episodes [START_REF] Dumont | Detection of myocardial ischemia with hidden Semi-Markovian models[END_REF].

In this paper, K HSMM are used to represent K different observation dynamics, associated with distinct physiopathological states or events, to be discriminated. A learning phase is firstly applied to each model, in order to estimate each λ k , ∀ k ∈ {1, 2, . . . , K}, from the sequences of observations, provided in a learning dataset. During this phase, λ k is initialized from the parameters of the equivalent HMM, in which the Viterbi algorithm is used to to estimate p j (µ d , σ d ). Additionally, a ij and π i are initialized with uniform probabilities and b i ( µ, Σ) is initialized by a Gaussian mixture model where the center of each Gaussian, which corresponds to the barycenter of each state in the observation space, have been initialized by the K-means algorithm. The Viterbi algorithm, extended to HSMM, is subsequently applied to obtain the final value of λ k , through an expectation-maximization stage. Learning is achieved when the log-likelihood

L k = log P (O 1:T |λ k ) (1)
converges to a maximum value, where

P (O 1:T |λ k )
is the probability that the observation sequence

O 1:T = O 1 , O 2 , . . . , O T is generated by the model with parameters λ k .
Once the learning phase is completed, the K HSMM are applied in a test phase to a set of observations on a specific dataset, in order to perform event detection. In this case, the log-likelihood for instant t and model k,

L k t = log P (O t-T +1:t |λ k ) (2)
is determined using a sliding window of length T . On-line detection of event α (α ∈ {k}) is finally performed when the following equation is verified:

L α t -L k t > δ α,k ∀ k ∈ {1, 2, . . . , K} -{α} (3)
where δ α,k are fixed thresholds that have to be optimized.

B. Quantization of observations

A quantization phase of the observations is proposed in this work to improve the detection performance of the HSMM approach. This quantization phase is justified, since i) HSMM are constituted of a discrete number of states M representing the dynamic range of the observations and ii) inputs to the proposed detector are the result of an automatic ECG segmentation phase that may produce segmentation errors.

Two quantization methods are applied: uniform, characterized by a constant quantization step ∆ QU and non uniform, characterized by a quantization step ∆ QNU , which depends on the distribution of the signal. Vector ∆ QNU is found by comparing the cumulative sum of values of the normalized histogram with respect to a threshold δ QNU . Different values of ∆ QU and δ QNU have been tested to obtain the value that maximizes detection performance. Variables noted with UQ and with NUQ indicate its uniform and non uniform quantization, as for RR UQ and RR NUQ .

C. Integration of delayed versions of the observation

Here, we explore the possibility of increasing the observability of the system while maintaining the same number of sources. In this sense, an observation matrix is constructed, integrating the original and a delayed version of the observed time-series. This observation matrix may be represented as:

O = O t-T +1:t O t-τ -T +1:t-τ ( 4
)
where τ is the predefined time delay. In the case of the RR series, the application of such a time-delay presents also a physiological argument, since it is known that the RR interval is modulated by the autonomic nervous function, with different dynamics for the sympathetic and parasympathetic systems. Different values of τ have been tested to optimize detection performance. These observation matrices will be represented in bold in this paper, such as RR.

D. Evaluation methodology

The proposed HSMM detection approach was firstly evaluated with simulated signals in order to optimize its structure and evaluate the impact of the different preprocessing methods in a controlled way. It was then evaluated with real signals acquired in NICU, in which AB episodes have been manually annotated.

1) Evaluation on simulated signals: Synthetic time series were generated with the FitzHugh-Nagumo model, described by the following differential equations [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF]:

dv dt = 3(v - 1 3 v 3 + r + I) dr dt = - 1 3 (v -a + 0.8r) (5)
Variables v and r are in a "rest" value until a perturbation is introduced into variable I. This induces an excursion in the phase space on variables v and r, before returning to their resting values. The simulated series have a size of 400 s and a perturbation I = -0.2 is introduced during the interval [300, 305] s. Parameter a was modified according to a uniform distribution in the interval a1 ∼ U(0.58, 0.62) and a2 ∼ U(0.78, 0.82) to generate series presenting two slightly different different dynamics. Moreover, the simulated series were divided by its maximum value, in order to generate signals with similar amplitudes. Gaussian white noise was added to the series to obtain a signal to noise ratio SNR = 5 dB. Responses generated with a1 and a2 are thus difficult to differentiate by only analyzing their instantaneous values (such as with a threshold-based detector).

Three synthetic datasets will be used for the learning phase:

• LSFN1: comprising 80 segments at "rest" (before the perturbation) generated with a1 (40 segments) and a2 (40 segments). • LSFN2: with 40 segments from the beginning of the perturbation and generated with a1. • LSFN3: represented by 40 segments from the beginning of the perturbation and generated with a2. All these segments have a duration of T = 10 s.

Once the HSMM have been trained, they are applied to a test dataset consisting of 400 synthetic signals of 400 s each: 200 of these signals were generated with parameter a1 and 200 with parameter a2. Equations 2 and 3 where applied as the online detection approach, with a sliding window of size T = 10 s.

2) Evaluation on real signals: RR series were extracted from the ECG of 32 preterm infants with frequent AB episodes as described in [START_REF] Altuve | Multivariate ecg analysis for apnoea-bradycardia detection and characterisation in preterm infants[END_REF], [START_REF] Dumont | Improving ecg beats delineation with an evolutionary optimization process[END_REF]. RR series were uniformly resampled at 10 Hz. 233 bradycardia episodes were manually annotated from 148 RR series. Two datasets were constructed for the learning phase:

• LS1: composed of 30 segments taken randomly and including an AB event, taken at the beginning of the bradycardia with duration of T = 7 s. • LS2: consisting of 300 segments at rest, taken randomly from the series, with duration of T = 7 s. The length of these segments (7 s) corresponds to the average time measured from the beginning of the bradycardia to the peak RR value within the bradycardia episode. In order to reduce the variability of the first sample of the series, the mean, determined within 5 s before the start of each segment, was removed for all segments. This is particularly important for the estimation of the first state (π i ) of the models.

The evaluation of the models trained from real signals was performed on a test set containing the totality of the 148 RR series (series duration = 26.25 ± 11.37 minutes) with 233 bradycardia episodes. The online HSMM detector is applied to these series as described in the previous section, with a sliding window of size T = 7 s.

3) Performance evaluation: For both simulated and real signals, true positives (T P ), true negatives (T N ), false positives (F P ) and false negatives (F N ) were determined for each sample by comparing the obtained detections with the available annotations. T P occur when a detection falls within a 20 s window, centered at a given annotation. Even if this window seems particularly large, we justify this choice by the fact that we are also evaluating the detection delay as an important marker for the selection of the optimal detection method.

Detection performance was evaluated by estimating the sensitivity (SEN = T P /(T P + F N )) and specificity (SP C = T P /(T N + F P )) of each detector, for different detection thresholds, and represented by means of ROC curves. The detection delay (dd) is also determined, and defined as the time elapsed between the annotation instant and the detection. Detection results were evaluated using the shortest distance to perfect detection (SDP D) according to Eq. 6. Detection delay results are presented in mean ± SD.

SDP D = min (1 -SEN ) 2 + (1 -SP C) 2 (6)
III. RESULTS

A. Simulated signals

This experience is particularly interesting, since, as described earlier, a simple threshold-based detection would be unable to discriminate the series produced with a1 from those produced with a2. Only the series generated with a1 were analyzed, as those generated with a2 do not provide additional information for our purpose of validation. K = 3 HSMM were employed to model each learning dataset: HSMM with set λ k to model LSFNk, k ∈ {1, 2, 3}. The BIC criterion was applied in each case to estimate the number of states per model: M k . A sliding window of size T = 10 s was used to obtain the observation sequence and calculate the log-likelihoods L k . Two thresholds (δ 2,1 and

δ 2,3
) were used to detect the perturbation generated with a1 as in Eq. 3. Detection results for simulated signals are shown in Table I. Optimal parameters were: τ = 5 s, ∆ UQ = 25 and δ QNU = 0.003.

B. Real signals

K = 2
HSMM models were used to detect AB episodes. Each parameter set λ k , was estimated to reproduce the dynamics of LSk, k ∈ {1, 2}. Only one threshold is thus used, δ 1,2 . M k was determined using the BIC criterion. A T = 7 s sliding window was used to construct O t-T +1:t and to determine L k t . The detection of an AB episode is obtained by applying Eq. 3.

ROC curves are shown in Fig. 1 and its respective mean detection delays are shown in Fig. 2. Optimal parameters were: τ = 0, 67 s, ∆ UQ = 1300 and δ QNU = 0.05. Conventional detection methods are also shown as a reference. In these figures, the point where SDP D occurs is shown with character "x". Table II lists the obtained SEN , SP EC and dd at the SDP D point (represented with the "x" symbol in figures 1 and 2), for the detection of AB episodes.

IV. DISCUSSION

Results on simulated signals show that the detection performance is improved when the quantization step is applied (see table I). The best detection performance and the minimum mean detection delay were obtained when a uniform quantization is applied (v UQ ), followed by the non uniform quantization version (v NUQ ). These results confirm the utility of the proposed HSMM approach, and measure the impact of the quantization step and the integration of a delayed version of the observables, on detection performance. The introduction of a time-delayed version of the signal may improve obsevability by synchronizing different effects that modify the same variable, which are not necessarily in phase. The gain in performance depends on i) the existence of such kind of effects on the available series and ii) the appropriate optimization of τ (which is considered constant). The gain obtained from simulated data is rather low, but results from real signals show the benefit of the approach.

Detection results from real signals (Fig. 1 and 2 and table II) show that modeling the dynamics of RR series by HSMM provides an improved detection performance, when compared to conventional detection methods (fixed or relative threshold detections). In addition, the application of the proposed preprocessing methods (RR UQ , RR NUQ and RR), further improves this performance and reduces the mean detection delay. These results indicate that our approach detects, on average, more than 2 s earlier than conventional detectors.

V. CONCLUSIONS AND FUTURE WORKS

In this paper, the difficult problem of exploiting the dynamics of biomedical time series for the early-detection of pathological events was addressed. A detection approach, based on an HSMM and combining different pre-processing methods, was proposed and evaluated using simulated and real signals. Results obtained from synthetic signals demonstrate the interest of the proposed approach, and the fact that the detection is mainly based on the difference of the temporal dynamics of the events under consideration.

On real data, the performance of the proposed approach for the detection of AB episodes has shown to be higher than classical methods, while minimizing the detection delay. The usefulness of the quantization phase and of the integration of delayed versions of the observation into the HSMM was also evaluated and showed to improve the detection performance. This is an important result, as the minimization of the delay from the beginning of the AB episode to the application of the therapy by the nurse is one of the main goals in this clinical field.

The main advantage of the approach proposed in this paper is its simplicity, opening the possibility of an embedded implementation into a mobile device. Future works are directed towards the multivariate extension of this method, in order to integrate other features extracted from the preterm infant's ECG, such as the R-wave amplitude, the QRS complex duration, the PR interval and P-wave morphology. This approach may also improve detection performance and reduce the detection delay.

Fig. 1 .

 1 Fig. 1. ROC curves for AB event detection.

Fig. 2 .

 2 Fig. 2. Mean detection delay of AB episodes.

This work was partly supported by the ECOS NORD project V09S04 M. Altuve is with