Anni Määttänen
email: anni.maattanen@aerov.jussieu.fr

Thierry Fouchet

O Forni

François Forget

H Savijärvi

Brigitte Gondet

R Melchiorri

Yves Langevin

V Formisano

M Giuranna

J.-P Bibring

A study of the properties of a local dust storm with Mars Express OMEGA and PFS data

Introduction

Dust is a ubiquitous component of the Martian atmosphere, and its source is the surface of the planet. Dust particles have a strong effect on circulation patterns via their impact on atmospheric heating rates [START_REF] Gierasch | The effect of dust on the temperature of the Martian atmosphere[END_REF]. They also affect the global climatic cycles of volatiles (CO 2 and H 2 O) via their role as cloud condensation nuclei [START_REF] Michelangeli | Numerical simulations of the formation and evolution of water ice clouds in the Martian atmosphere[END_REF][START_REF] Rodin | Dynamical properties of Mars water ice clouds and their interactions with atmospheric dust and radiation[END_REF]. Dust distribution itself is affected by the atmospheric circulation, and cloud formation redistributes dust via sedimentation. These feedback mechanisms between dust distribution, radiative balance, atmospheric temperatures and circulation, and cloud formation are very nonlinear in nature. A good summary on the Martian dust cycle can be found in [START_REF] Kahn | The Martian dust cycle[END_REF].

Mars hosts a large variety of dust phenomena, starting from the small-scale dust devils going up to planetary scale dust storm events. Surveys and classifications of storms of all scales and dust devils observed in visible spacecraft images can be found in [START_REF] Briggs | Viking Orbiter imaging observations of dust in the Martian atmosphere[END_REF], [START_REF] Cantor | Martian dust storms: 1999 Mars Orbiter Camera observations[END_REF] and [START_REF] Cantor | Mars Orbiter Camera observations of Martian dust devils and their tracks (September 1997 to January 2006) and evaluation of theoretical vortex models[END_REF]. In the thermal infrared, [START_REF] Smith | Interannual variability in TES atmospheric observations in Mars during 1999-2003[END_REF] have monitored the mean dust optical thickness during three consecutive Martian years, while [START_REF] Smith | The annual cycle of water vapor on Mars as observed by the Thermal Emission Spectrometer[END_REF] reviewed specifically the planet-encircling dust storm of 2001. However, the two visible surveys [START_REF] Briggs | Viking Orbiter imaging observations of dust in the Martian atmosphere[END_REF][START_REF] Cantor | Martian dust storms: 1999 Mars Orbiter Camera observations[END_REF] present discrepancies in their spatial distributions of dust storms that are currently unexplained: the Mars Orbiter Camera survey (MOC, [START_REF] Cantor | Martian dust storms: 1999 Mars Orbiter Camera observations[END_REF] found fewer low-latitude dust storms than the Viking survey [START_REF] Briggs | Viking Orbiter imaging observations of dust in the Martian atmosphere[END_REF], and MOC also found most of the regional dust storms to be located in the Northern Hemisphere in contrast with the historical record. [START_REF] Cantor | Martian dust storms: 1999 Mars Orbiter Camera observations[END_REF] explained this difference as arising from a temporal change in dust One important factor affecting the expansion of storms is the radiative effect of dust on atmospheric circulations: suspended dust warms the atmosphere and the warming enhances the circulations. However, the warming also tends to stabilize the vertical temperature profile of the atmosphere, suppress convection, and smooth out horizontal temperature gradients inside the storm. In the end [START_REF] Cantor | Martian dust storms: 1999 Mars Orbiter Camera observations[END_REF] pointed out that this effect can lead to a decay of the storm via a self-destructive mechanism. However, mechanical turbulence caused by wind shear can sustain dust lifting even in an isothermal situation, if the surface winds are strong. Modeling studies of the Martian dust cycle have been conducted by, e.g., [START_REF] Murphy | Three-dimensional numerical simulation of Martian dust storms[END_REF]; Newman et al. (2002a,b); [START_REF] Basu | Simulation of the Martian dust cycle with the GFDL Mars GCM[END_REF][START_REF] Basu | Simulation of spontaneous and variable global dust storms with the GFDL Mars GCM[END_REF]; [START_REF] Kahre | Simulating the Martian dust cycle with a finite surface dust reservoir[END_REF][START_REF] Kahre | Modeling the Martian dust cycle and surface dust reservoirs with the NASA Ames general circulation model[END_REF], and [START_REF] Kahre | Investigations of the variability of dust particle sizes in the martian atmosphere using the NASA Ames General Circulation Model[END_REF].

The radiative effects of dust are dictated by the size, composition, and optical properties of the dust particles. Several authors using a multitude of observations have analyzed these properties in detail or partly, depending, for example, on the wavelength range of the instrument used, and the type of the used observations [START_REF] Pollack | Properties and effects of dust particles suspended in the Martian atmosphere[END_REF][START_REF] Ockert-Bell | Absorption and scattering properties of the Martian dust in the solar wavelengths[END_REF][START_REF] Forget | Improved optical properties of the Martian atmospheric dust[END_REF]Tomasko et al., 1999;[START_REF] Markiewicz | Optical properties of the martian aerosols as derived from imager for mars pathfinder midday sky brightness data[END_REF][START_REF] Fedorova | Remote sounding of the martian dust from ISO spectroscopy in the 2.7 μm co 2 bands[END_REF][START_REF] Clancy | Mars aerosol studies with the MGS TES emission phase function observations: Optical depths, particles sizes, and ice cloud types versus latitude and solar longitude[END_REF][START_REF] Wolff | Constraints of the size of Martian aerosols from Thermal Emission Spectrometer observations[END_REF][START_REF] Mateshvili | Martian ice cloud distribution obtained from SPICAM nadir UV measurements[END_REF]Vincendon et al., 2007). The frequently used single-scattering properties of dust in near-infrared are those presented by [START_REF] Ockert-Bell | Absorption and scattering properties of the Martian dust in the solar wavelengths[END_REF]. They, however, did not present any values in the 2.5 -2.9 μm CO 2 band, since they were not able to retrieve the surface albedo at these wavelengths. Within this band, the analysis of [START_REF] Fedorova | Remote sounding of the martian dust from ISO spectroscopy in the 2.7 μm co 2 bands[END_REF] revealed a strong absorption feature centered at 2.8 μm. [START_REF] Fedorova | Remote sounding of the martian dust from ISO spectroscopy in the 2.7 μm co 2 bands[END_REF] attributed this absorption to the presence of phyllosilicates in the dust. Indeed, new observations from Mars Express and

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

Mars Odyssey have shown the presence of phyllosilicates also on the surface of the planet [START_REF] Bibring | Mars surface diversity as revealed by the OMEGA/Mars Express observations[END_REF][START_REF] Poulet | Phyllosilicates on Mars and implications for early martian climate[END_REF][START_REF] Chevrier | Early geochemical environment of Mars as determined from thermodynamics of phyllosilicates[END_REF].

In this study we present the analysis of a local dust storm observed by two Mars Express instruments. OMEGA (Observatoire pour la minéralogie, l'eau, les glaces et l'activité) was used to retrieve the dust single-scattering albedo in the range 0.4 -4.0 micron, and to determine the structure of the optical thickness within the dust storm. The Planetary Fourier Spectrometer (PFS) was used to determine the impact of the dust storm on the atmospheric temperatures. We also compared these data with 1-D atmospheric model calculations modeling the effect of dust on the atmospheric temperature field.

Section 2 introduces the instruments and the data. Then, in Section 3, we present the different methods used to identify the dust storm in the OMEGA and PFS data. Section 4 presents the method for retrieving the single-scattering albedo as well as mapping the dust optical thickness, and discusses the results.

Section 5 introduces the methods of retrieval and retrieved temperatures from the PFS data, discusses the temperature structure inside the storm, and compares the results to model runs with a 1-D atmospheric column model for Mars. In Section 6, we summarize the results. In this section, we present the two Mars Express instruments involved in this study, we briefly discuss the calibration of their data, and we describe the data selection process.

Observations

Mars

The instruments

The OMEGA spectro-imaging instrument and calibration of data

OMEGA is a visible and near-IR mapping spectrometer, with coaligned channels working in the 0.38 -1.05 μm visible and near-IR range (VNIR channel) and in the 0.93 -5.1 μm short wavelength IR range (SWIR channel). The data products constitute three-dimensional (x, y, λ) image-cubes, with two spatial and one spectral dimensions. The spectral resolution is 7 nm in the VNIR channel, 14 nm in the 0.93 -2.73 μm and 28 nm in the 2.55 -5.1 μm range.

An extensive description of OMEGA can be found in [START_REF] Bibring | Omega: Observatoire pour la minéralogie, l'eau, les glaces et l'activité[END_REF].

OMEGA's photometric calibration was carried out in a laboratory using a tunable black body with temperature stabilized up to 1200 • C and a cold black body (down to 70 K). Additional in-flight calibration were carried out by comparing previous observations and models of Phobos reflectance with OMEGA observations. The absolute photometric calibration is reliable down to the 15% level, while the relative calibration between two spectra is better than 5%. The typical integration time, defined by the spacecraft ground-track

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
velocity and the spatial sampling chosen, is 2.5 ms or 5 ms. With such integration times, a signal-to-noise ratio of higher than 100 is obtained over the entire infrared spectral range, and about 50 in the VNIR channel.

The PFS spectrometer and calibration of data

The Planetary Fourier Spectrometer is an infrared double-pendulum interferometer working in two different channels. The Short Wavelength Channel (SWC) covers wavenumbers from 1700 to 8200 cm -1 (1.2 -5.9 μm) while the Long Wavelength Channel (LWC) operates in the range 250 -1700 cm -1

(5.9 -40 μm). In the present analysis, we have used data from LWC only.

The LWC is characterized by an unapodized spectral resolution of 1.3 cm -1 , a sampling step of 1.0 cm -1 and a field of view of 50 mrad (FWHM). This corresponds to 12 km on Mars surface when viewed from periapsis. A more detailed description of the PFS instrument can be found in [START_REF] Giuranna | Calibration of the Planetary Fourier Spectrometer long wavelength channel[END_REF].

The PFS-LW channel calibration has been presented in detail by [START_REF] Giuranna | Calibration of the Planetary Fourier Spectrometer long wavelength channel[END_REF]. Within the wavenumber range used to retrieve the temperature profile (500 -800 cm -1 ), the noise equivalent spectral radiance (NESR) of a single spectrum varies between 7 × 10 -8 and 1.5 × 10 -7 W cm -2 sr -1 /cm -1 , corresponding to a signal-to-noise ratio of 50 -100 on a mid-latitude daytime spectrum.

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

Data selection

In the beginning of the mission the Mars Express ground-tracks overlapped each other every eleven orbits, i.e. three Martian days apart (later on the resonance has been changed). This property allowed for searching rapid temporal changes in the Martian spectrum that could be attributed to dust lifting events over a limited geographical area (see Section 3 for the spectral criteria).

Within a few of such events, we chose two orbits, #1201 and #1212 (Ls = 135 • , MY27), which overlap south of the equator (Fig. ??). The overlap was focused on observing a surface feature in the Pollack crater and subsequently provided us with a dataset suitable for studying transient atmospheric features. Such features were found, since the latter orbit, #1212, showed strong evidence, in both OMEGA and PFS spectra, for the presence of an optically thick, local dust storm centered at 3 • S, 24.7 • E, at L s = 135 • and 13.30 local time.

[Figure 1] 3 Identification of the local dust storm Still in the solar-reflected range, the contrast (ratio of the band depth to the continuum level) of the strong CO 2 bands at 2.7 and 2.0 μm decreases from the former orbit to the latter, while the weak CO 2 bands at 1.6 and 1.4 μm completely disappear in the spectrum obtained in orbit #1212. In the thermal infrared, at wavelengths longer than 4.4 μm, the brightness temperature decreases from 275 K (#1201) down to 235 K (#1212). There is also a hint of an atmospheric warming revealed by the increased brightness temperature within the CO 2 band at 4.35 μm.

Spectroscopic recognition of atmospheric dust

[Figure 2] These spectral differences between OMEGA spectra recorded about three Martian days apart can only be attributed to the presence of suspended atmospheric dust at the orbit #1212. The Martian dust is known to be bright, with large single-scattering albedo [START_REF] Ockert-Bell | Absorption and scattering properties of the Martian dust in the solar wavelengths[END_REF]. The presence of atmospheric dust above low-albedo terrains thus increases the observed reflectance. The decrease in contrast of the CO 2 bands can also be attributed to the presence of dust. In fact, the reflectance factor in the core of the CO 2 2.7-μm absorption band constitutes one of the best indicators of the presence of atmospheric dust. The core of this band is saturated for surface pressures greater than 200 -300 Pa (Titov et al., 2000), which means that all radiation is absorbed by the CO 2 before reaching the surface and none escapes to space: thus the reflection equals zero in the band center if the atmosphere is clear of scatterers. Inside a saturated band, the observed reflection is nearly insensitive to surface albedo. If any reflection is seen in the band center, it must come from atmospheric aerosol particles that scatter the radiation before it is absorbed by CO 2 . The same effect explains the decrease in contrast of the

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
CO 2 absorption bands centered at 1.4, 1.6, and 2.0 μm, although these bands are not saturated. Finally, the decrease in the brightness temperature at wavelengths longer than 4.4 μm can be explained by a combination of two effects: first, the surface temperature decreases as the atmospheric dust absorbs the solar flux, and second, the thermal emission observed during an optically thick dust storm comes from atmospheric levels located at high altitudes, where the temperatures are colder than at the surface.

In Fig. ??, we compare the PFS spectra taken during the two orbits at 24.7 • E of longitude and between 4 • S and 2 • S. The brightness temperature of the continuum outside the CO 2 band centered at 668 cm -1 decreases by about 70 K from orbit #1201 to orbit #1212. In contrast, the brightness temperature in the core of the CO 2 band increases by about 10 K. Again, this behavior of the thermal spectrum is typical of the occurrence of a local dust storm which cools the surface and the lower atmosphere while it heats up the midatmosphere. The Martian dust presents an absorption feature at 1100 cm -1 due to a silicate vibration mode [START_REF] Smith | Separation of atmospheric and surface spectral features in Mars Global Surveyor Thermal Emission Spectrometer (TES) spectra[END_REF]. Unfortunately, this absorption feature cannot be detected on orbit #1212 as the PFS spectrum in dusty conditions is too noisy in the range 1000 -1300 cm -1 . Therefore, we are not able to measure the dust optical depth in the solar reflected and thermal range simultaneously.

[Figure 3]

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

Independent Component Analysis

As discussed in the previous section, the spectra obtained during the orbit #1212 from both OMEGA and PFS-LW clearly show the presence of a local dust storm between 4 • S and 2 • S. However, we will be in position to retrieve the dust parameters from the OMEGA spectrum only if the measured reflectance is independent of the surface albedo. In order to assess this possibility, we have performed an Independent Component Analysis (ICA) on the OMEGA cubes centered on the region of interest. ICA is a method that can retrieve the separate sources that, mixed together, constitute the observed signal. The method is described in detail in the Appendix A. Our goals were: i) to separate the dust source from the surface albedo source and to estimate their relative contribution to the OMEGA spectra, ii) to map the pixels were the dust source dominates the observed signal, hence to map the geographical extension of the dust storm. We applied the ICA method to the logarithm of the observed reflectance during orbit #1212 in the spectral range 2.6 -3.0 μm. Our choice of the logarithm of the observed reflectance is justified by the fact that the dust reaches high altitudes, where the growth of the CO 2 transmissivity is linear in the logarithmic scale. The use of the logarithm of the reflectance also enhances the noise, but since the signal in the core of CO 2 band is high due to dust scattering, this does not pose problems in our case. corresponds to the surface albedo. We estimated the contributions of the different components to the observed signal by using the quadratic mean of the mixing coefficients. For the ten brightest pixels in the core of the dust storm, the dust component contributes to 99% of the power. Within the whole dust storm, the dust component contributes to about 85% of the total observed signal. Hence, the ICA method proves that spectra within the core of the dust storm are independent of the second major component, i.e. the surface albedo, and thus our single-scattering analysis will be valid.

[Figures 4]

Optical properties of dust

Radiative transfer code

We have modeled the OMEGA spectrum of the local dust storm in the solarreflected range using a radiative transfer code including the gaseous absorption of CO 2 , CO and H 2 O and taking into account the multiple scattering events. Although CO and H 2 O are minor absorbers, their bands could be misinterpreted as dust mineralogical signatures if they were not included in the calculation. Multiple scattering was handled using a two-stream approximation [START_REF] Meador | Two-stream approximations to radiative transfer in planetary atmospheres: A unified description of existing methods and a new improvement[END_REF]. The gaseous absorptions were first calculated in a line-by-line approach using the spectrometric parameters extracted from the GEISA databank [START_REF] Jacquinet-Husson | The 1997 spectroscopic geisa databank[END_REF]. The line-by-line opacities were then used to calculate correlated-k distributions. The modeled atmosphere was divided in 45 layers equally spaced in log(pressure). For each layer we specified in entry the gaseous abundance, the dust optical thickness at 1 μm, and the single-scattering albedo, the asymmetry parameter and the

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

extinction cross-section as a function of wavelength. The gaseous abundances were set to the standard volume mixing ratios of 0.95 and 8 × 10 -4 for CO 2 and CO, and 1.5 × 10 -4 for H 2 O as derived from PFS observations on orbits 1201 and 1212 by Fouchet et al. (2007).

Retrieval methods

Single-scattering albedo

Our analysis method for the dust optical properties stems from the practically complete decoupling of the surface and the dust reflectance in the case of a very thick dust storm. We first averaged the ten pixels with the highest reflectance in the SWC located within the dust storm (as identified by the ICA method).

Then each spectel (i.e. spectral individual element) of the averaged spectrum was fitted by adjusting only the dust single-scattering albedo. Except for the dust optical thickness, all other parameters remained fixed to the following values:

• The surface albedo was retrieved from the overlapping orbit #1201, assuming a dust-free atmosphere (τ = 0.0).

• The g and Q ext dust parameters were fixed to the [START_REF] Ockert-Bell | Absorption and scattering properties of the Martian dust in the solar wavelengths[END_REF] values interpolated to the OMEGA spectral resolution.

• The temperature profile was retrieved from PFS data orbit #1201 (see section 5.2). We have estimated the sensitivity of our retrieved ω 0 to several parameters:

• The surface albedo was determined assuming a dust-free atmosphere (τ = 0.0). This assumption is supported by Spiga et al. (2007), who analyzed the 2.0-μm CO 2 band to determine a low dust content of τ = 0.1-0.2. However, the Mars Exploration Rover (MER) observations yielded a dust content of τ = 0.6 at the time of the orbit 1201 of Mars Express [START_REF] Wolff | Constraints on dust aerosols from the Mars Exploration Rovers using MGS overflights and Mini-TES[END_REF].

We calculated that the higher dust content would lead to a 5% change in the retrieved surface albedo compared to the dust-free case, leading to a change of less than 0.1% in the determination of ω 0 . This analysis provides quantitative evidence for the decoupling of the reflectance from the surface albedo whenever the dust optical thickness is large, hence validating our method.

• We assumed an asymmetry parameter of g = 0.63, as suggested by most studies [START_REF] Ockert-Bell | Absorption and scattering properties of the Martian dust in the solar wavelengths[END_REF][START_REF] Clancy | Mars aerosol studies with the MGS TES emission phase function observations: Optical depths, particles sizes, and ice cloud types versus latitude and solar longitude[END_REF]. This corresponds to particles of about 1 μm in radius. It may vary within a few percent, but as shown by the right panel of Fig. ?? this would affect our retrieved singlescattering albedo less than 0.5%.

• For the analysis we used a temperature profile determined from the dustfree orbit #1201. This may not correspond to the atmospheric conditions

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
prevailing in the core of the dust storm. However, using warmer profiles, we found that our retrieved single-scattering albedo is insensitive within 0.1% to the assumed temperature profile.

• Our method to take into account the multiple scattering events, the twostream method, is relatively simple and may not be the most accurate. We have compared the reflectance calculated with our method with the simulations made for a standard Mars atmosphere by (Tran and Rannou, 2004, http://tttrung.free.fr/) using 3D spherical Monte-Carlo algorithm, which were intended to produce as realistic results as possible. We found our calculations to lie within 5% of the true reflectance, which leads to an error of +0.004 -0.002 on ω 0 . Moreover, Vincendon et al. ( 2007) also determined the dust single-scattering albedo from the same OMEGA data, but using the 'infinite' dust optical depth approximation with a Monte-Carlo algorithm. Their ω 0 values differ from ours by less than a percent.

• The signal-to-noise ratio for the high reflectance spectra used in this study is larger that 200. Hence, the largest error source comes from the uncertainty on the OMEGA absolute calibration, estimated to lie within ±15%. As shown on the right panel of Fig. ??, this corresponds to an uncertainty on the determination of the dust single-scattering albedo of +0.02 -0.03 in the center of 2.9 μm band and +0.012 -0.006 in the rest of the spectral interval.

Finally, adding quadratically all the errors, our uncertainty on ω 0 amounts to +0.021 -0.031 in the center of 2.9 μm band and to +0.014 -0.008 in the rest of the spectral interval.

A C C E P T E D M

A N U S C R I P T ACCEPTED MANUSCRIPT

Dust optical thickness

Outside the core of the dust storm, the reflectance varies with the pixels. This variation cannot be attributed to a change in the dust properties, but rather to a change in the dust optical thickness: surface and dust contributions to the reflectance can no longer be considered decoupled. Using the surface albedo retrieved from orbit 1201, this situation allows us to map the optical thickness outside the core of the dust storm. We proceeded as follows:

• We used the dust single-scattering albedo retrieved previously.

• The surface pressure was calculated for each pixel using the European Mars Climate Database, EMCD (http://www-mars.lmd.jussieu.fr/, [START_REF] Lewis | A climate database for Mars[END_REF]). Since the spatial grid of the EMCD is much coarser than the resolution of the OMEGA images, the pressure was hydrostatically extrapolated to each pixel using the MOLA topography.

• The surface albedo was retrieved from observations obtained during the orbit 1201. The observed area was divided in three albedo regions: low, intermediate (with a strong gradient), and high. The region-averaged values at 1 μm were 0.16 (low), 0.35 (intermediate), and 0.37 (high). The areas were limited by latitudes 4.7 • S and 0.4 • S.

• Then, finally the optical depth was adjusted to fit the observed reflectance for every wavelength of each pixel.

The errors in the determination of the dust optical thickness can be due to small-scale variations of the surface albedo and to the propagation of uncertainties in the dust single-scattering albedo ω 0 . We have estimated the smallscale variations of the surface albedo by calculating the standard deviation of the observed reflectance in the three albedo regions. The resulting uncertain-

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
ties depends on the dust optical thickness. At the low-end limit, dust opacities lower than τ = 1.0 for the low albedo region, τ = 2.3 for the intermediate, and τ = 3.0 for the high albedo region cannot be distinguished from small-scale surface albedo variations. They hence constitute our lower limits of sensitivity.

For larger opacities, the errors caused by the surface albedo decrease rapidly with τ : ±0.2 for τ = 4, ±0.1 for τ = 6, and are negligible for τ > 7.

The uncertainties on the OMEGA absolute calibration, and the accuracy of the two-stream method propagate into uncertainties in the retrieved singlescattering albedo ω 0 . This also causes an uncertainty in the dust optical thickness, which increases with decreasing τ . We have calculated this error to be +0.5

-0.7 at τ = 8, quite constant, +0.5 -0.9 , for smaller values, and decreasing towards τ = 10.

Summing quadratically the errors give error bars of +0.5 -0.9 for τ , translating to +13.5/-23% for τ = 4, +8.5/-15.2% for τ = 6, and +6.75/-11.5% for τ = 8.

For a summary of the errors in the τ retrieval, see Table 1.

The analysis proved to be unreliable in the longer wavelength channel for two reasons. Due to the presence of thermal emission, it is difficult to retrieve the surface albedo from the orbit 1201 at wavelengths longer than 3.5 μm. In addition, in orbit 1212, thermal emission will also contribute significantly to the spectrum as the optical depth decreases. Around 3 μm, the atmospheric dust and the surface reflectivity are very similar, preventing us from separating surface and atmospheric dust contributions. Thus we are unable to distinguish variation of the surface albedo from variation of the dust optical depth at longer wavelengths, and our dust optical thickness retrieval is not reliable above 3 μm. depend on the assumed size distribution. In contrast, the OMEGA experiment is mostly sensitive to the spectral properties and insensitive to the particle radius.

Centered at 2.9 μm, our retrieved single-scattering albedo shows a strong spec-

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
tral feature, attributed to hydrated minerals. This is not a surprise, since this hydration band is also ubiquitous in the spectrum of Martian soils. However, the shape of our band differs from that retrieved from ISO data by [START_REF] Fedorova | Remote sounding of the martian dust from ISO spectroscopy in the 2.7 μm co 2 bands[END_REF]. Their analysis revealed a narrow absorption feature centered at 2.8 μm, characteristic of phyllosilicates. Phyllosilicates are also observed in OMEGA spectra [START_REF] Bibring | Mars surface diversity as revealed by the OMEGA/Mars Express observations[END_REF][START_REF] Poulet | Phyllosilicates on Mars and implications for early martian climate[END_REF], but on specific geographical units covering altogether a small fraction of the planet's surface.

Retrieving the Martian dust single-scattering albedo in the 2.6 -2.9 μm strong CO 2 band is difficult, since the dust scale height affects the analysis. We consider our determination of the single-scattering albedo in this spectral region more reliable than that of [START_REF] Fedorova | Remote sounding of the martian dust from ISO spectroscopy in the 2.7 μm co 2 bands[END_REF] for the following reasons: Besides the signature of hydrated minerals at 2.9 μm, the Martian dust is devoid of any mineralogical signature. There is a hint of ferric oxide absorption at 0.9 μm, but at the limit of the OMEGA relative calibration. The rest of the spectrum shows no spectral absorption at 0.5%-level in the 1.0 -2.5 μm range, and at the 1%-level in the 2.7 -4.0 μm range.

An interesting result is that altogether the Martian dust single-scattering albedo spectrum looks very similar to that of Saharan desert dust with par-

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
ticle size of approximately 1 μm (Tegen and Lacis, 1996). This particle size is close to the average effective radius of Martian dust, and the asymmetry parameter g = 0.63 we have chosen for the study corresponds to particles of about this size.

Dust optical depth

The analyzed τ in the continuum at the wavelength 1 μm is shown in Fig. ??.

The values of the retrieved τ fall in the full range of 0.0 -9.9 in this image that is zoomed to the storm area. However, it should be kept in mind that in this albedo region our method is reliable only for τ > 2.3 (see Section 4.2.2).

The storm shows a lot of structure in the range of τ analysed. Several optically thick, roundish dust plumes can be seen, and their diameters, when defined with the boundary τ = 5 (color change from blue to green in the image), vary from some kilometers to about 10 km, although the estimates are limited by the width of the OMEGA session. The latitudinal extent of the storm area imaged in Fig. ?? is 4.24 degrees, which corresponds to about 250 km.

The width of the OMEGA image is about 0.36 degrees, corresponding to approximately 20 km, but the storm seems to be extending outside the image, and thus we can not evaluate its exact extent in the longitudinal direction.

[Figure 7]

We have estimated a range for the lifting rate of the dust in the storm assuming that 1) the storm has started right after the previous observation at orbit #1201, i.e. 3 sols earlier (minimum lifting rate), or that the storm started 2 hours before the observation (maximum lifting rate), and 2) that the lifted dust has not been advected away from the source region, i.e., the source region

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
is the area where the optical thickness for the lifting rate estimate is calculated.

For the lifting rate estimation we used the region of the storm seen by the PFS FOV for the temperature retrieval (the oval in Fig. ??) for purposes of later temperature profile comparison. The mass of the lifted dust can be estimated from the optical thickness using the relation M = 4 3 ρrτ Qext , where M is the dust mass in the unit column of air, ρ the dust density, r the particle radius, τ the dust optical thickness and Q ext the extinction coefficient. The values of τ and Q ext should be defined on the same wavelength, here at 1 μm, the ρ was taken as 2500 kg m -3 and r=1.5 μm. The average optical thickness at 1 μm in the area of the PFS FOV turned out to be τ avg =4.23. Using the Q ext at 1 μm from [START_REF] Ockert-Bell | Absorption and scattering properties of the Martian dust in the solar wavelengths[END_REF] the above formula gives (0.25 -9.1) × 10 -7 kg m -2 s -1 for the lifting rate range within the time frame of 2 hours -3 sols. For comparison, some examples of averaged lifting rates calculated from modelling are of the order of 10 -9 kg m -2 s -1 , (globally averaged, [START_REF] Basu | Simulation of the Martian dust cycle with the GFDL Mars GCM[END_REF], 10 -7 -10 -3 kg m -2 s -1 (zonally averaged, [START_REF] Kahre | Investigations of the variability of dust particle sizes in the martian atmosphere using the NASA Ames General Circulation Model[END_REF].

The storm is located on a topographical slope (highland in the south, lowland in the north) at the edge of an area of an abrupt albedo change. OMEGA data was used in retrieval of surface pressure (Forget et al., 2007;Spiga et al., 2007), and in their analysis Spiga et al. (2007) noted some occurrences of atmospheric waves via their signature in the surface pressure field. One occurrence of these waves was observed in orbit 1201, only some degrees away from the location of the dust storm observed 3 sols later. The waves occurred at a large crater centered at 8 • S, whereas the storm is centered at 3 • S.

Even though the dust plumes appear similar in shape as convective clouds, we can not compare the results with convective cell sizes from, for example, LES simulations [START_REF] Michaels | Large-eddy simulations of atmospheric convection on Mars[END_REF]Tyler et al., 2008): convection

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

in such simulations has been traditionally modeled in low-dust situations, and even though such simulations do produce vortices that could lift dust (dust devils), they do not initiate dust storms. Comparison with global modelling [START_REF] Murphy | Three-dimensional numerical simulation of Martian dust storms[END_REF]Newman et al., 2002a,b;[START_REF] Basu | Simulation of the Martian dust cycle with the GFDL Mars GCM[END_REF][START_REF] Basu | Simulation of spontaneous and variable global dust storms with the GFDL Mars GCM[END_REF][START_REF] Kahre | Simulating the Martian dust cycle with a finite surface dust reservoir[END_REF][START_REF] Kahre | Modeling the Martian dust cycle and surface dust reservoirs with the NASA Ames general circulation model[END_REF][START_REF] Kahre | Investigations of the variability of dust particle sizes in the martian atmosphere using the NASA Ames General Circulation Model[END_REF] is difficult because such studies focus on initiation mechanisms, details of the dust lifting parameterizations, and modeled lifting rates, but do not study the smaller scale structure of the storms in more detail, naturally also in lack of resolution. Toigo et al. ( 2002) modeled observed dust lifting events in the southern polar areas, but focused on the aforementioned features as well.

Thus we conclude that to study further the structure of the storm, and unveil the mechanisms responsible for it, a high-resolution mesoscale model study is required, and even though it is out of the scope of this paper, our analysis gives very good constraints for such a study.

The effect of dust on the atmosphere

The temperature structure analysis from PFS observations

The temperature structure retrieval is based on the fitting of the ν 2 band of CO 2 centered at 667 cm -1 as described in detail in Fouchet et al. (2007). The forward radiative transfer modeling is carried out using a line-by-line model including the gaseous absorption by CO 2 , but without taking dust into account since we cannot determine the dust opacity from the 9 μm band. In large optical depth conditions, this will tend to underestimate the temperature in the lower atmosphere. The inversion of the temperature profile is

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
carried out using the regularization and stabilization methods described by Fouchet et al. (2007). Using this technique, vertical temperature profiles can be retrieved from the surface up to a pressure level of 10 Pa (∼ 35 km) with a vertical resolution of about one scale height. The uncertainties on the retrieved temperature for a single PFS spectrum typically amount to ∼ 2 K from 10 to 35 km, but increase in the first scale height above the surface to reach about ∼ 5 K at the deepest pressure level. The thermal fingerprint of the local dust storm is evident at 3 • S in the temperature structure of orbit #1212, with a net cooling close to the surface and net heating in the upper atmosphere, as expected from theoretical considerations [START_REF] Gierasch | The effect of dust on the temperature of the Martian atmosphere[END_REF]. In the core of the dust storm, the temperatures near the surface and at 30 Pa, were ∼ 10 K colder and ∼ 5 K warmer, respectively, than at the surrounding unperturbed latitudes (Fig. ??). Several measurements of temperature profiles have been performed in the millimetric range or in the thermal infrared at the occasion of regional or global dust storms. [START_REF] Moreno | Wind measurements in Mars' middle atmosphere at equinox and solstice: IRAM Plateau de Bure interferometric CO observations[END_REF] and [START_REF] Clancy | Mesospheric winds and temperatures from JCMT sub-millimeter CO line observations during the 2003 and 2005 Mars oppositions[END_REF] out of a regional dust storm at the end of MY27. In more drastic conditions, [START_REF] Hanel | Infrared Spectroscopy Experiment and the Mariner 9 mission: Preliminary results[END_REF] and [START_REF] Smith | Recovery of surface reflectance spectra and evaluation of the optical depth of aerosols in the near-ir using a Monte Carlo approach: Application to the OMEGA observations of high-latitude regions of Mars[END_REF] using spacecraft monitoring of the CO 2 -band, and [START_REF] Moreno | Wind measurements in Mars' middle atmosphere at equinox and solstice: IRAM Plateau de Bure interferometric CO observations[END_REF] in the millimetric range, found a vertical temperature profile close to isothermal during a global dust storm.

Retrieved PFS temperature profiles

Therefore, atmospheric dust tends to homogenize the vertical temperature profile. We observe the same trend in our temperature profiles (Fig. ??), but with a weaker amplitude than for the regional or global dust storms described above.

[Figure 9]

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

Comparison with 1-D atmospheric models

The 1-D version of the general circulation model of the Laboratoire de Météorologie Dynamique (hereinafter LMD 1-D model, [START_REF] Lewis | A climate database for Mars[END_REF] We conducted a model run for the location and season of the dust storm with the average dust optical thickness of τ =4.23, calculated for the area of the PFS FOV, and a well-mixed dust vertical profile with top of the dust at 55 km.

The model results at 13.30 LT exhibit a very warm, nearly isothermal middle atmosphere (the short-dashed line in Fig. ??). There is a clear difference between the retrieved PFS temperature profile (the solid line in Fig. ??) and the profile predicted by the model (short-dashed line in Fig. ??). The warming of middle atmosphere in the PFS profile seems to be only 5 K and the cooling of the surface about 10 K, whereas the model profile is significantly warmer in the middle atmosphere (60 K at 1 Pa), and remains isothermal between 1 Pa The radiative response time of the Martian atmosphere is quite short (shorter than the three sols between the two orbits studied here). Hence, it is unlikely that the difference between the PFS-retrieved and the modeled temperature vertical profile can be attributed to the fact the atmospheric column did not reach thermal equilibrium. This was tested by looking at the results of the first sol of the aforementioned run with dust top at 55 km (see Fig. dust exhibit an inversion of differing strengths near the surface, and the shape of the profile is dictated by strong warming in the dust layer. The general shape of the profiles is, however, less steep than the observed profile, and changing the dust top altitude only moves this shape with altitude. So we can get the overall values of the modeled temperature profile to agree with that observed, in particular in the middle atmosphere, where also the PFS resolution is adequate, but the slope of the profile is not fully satisfactory.

The best overall agreement seems to be found at dust top altitude around 10-20 km.

However, it should be noted that the PFS can not resolve the details of the temperature profile near the surface, so the fact that we do not see the inversion in the PFS profile does not mean it would not exist. The retrieval procedure of the PFS temperature profile suffers from the thick dust, which degrades the reliability of the retrieval, and thus the temperature vertical profile retrieved from PFS data in the core of the dust storm must be considered with caution. As mentioned before, since we cannot retrieve the dust opacity in the thermal infrared from the silicate band at 9 μm due to a too low signalto-noise ratio, we inverted the temperature profile from the CO 2 band for a clear atmosphere. This tends to underestimate the temperature in the lower layers. In addition, the surface temperature and the atmospheric temperature in the first half scale height are highly correlated, which hampers an accurate inversion close to the surface.

The confinement of the dust in the first scale height would be difficult to reconcile with OMEGA data. In the well-mixed infinite dust model, an optical thickness of 10 is required to fit the 2.0-μm CO 2 band in the dusty spectrum of orbit 1212. If the dust is confined to the first 10 km, an optical thickness of

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
20 is required, which seems an extremely high value.

Faced with these contradictions and uncertainties, we are not in position to firmly conclude, as suggested by the 1-D atmospheric model, that the dust is confined close to the surface. Other OMEGA and PFS orbits will be analyzed to investigate this issue in a future work.

Summary and conclusions

We have used the observations of two Mars Express instruments, OMEGA and PFS, to distinguish and analyze a local dust storm near the equator of Mars in late northern summer. The same geographical unit was observed on two orbits only 3 sols apart, and the later orbit revealed a local dust storm, whereas in the first orbit the atmosphere was very clear. These two orbits were exploited to study the optical properties of the dust and the morphology of the dust storm.

Since the dust storm is optically very thick, dust reflection is decoupled from the surface albedo and all scattering comes from the dust contribution. In this case we were not limited to the CO 2 absorption bands in the retrieval, but were able to derive the dust single-scattering albedo reliably throughout the wavelength range examined. The single-scattering albedo is a basic parameter required for radiative transfer codes that are used in many applications, for example in atmospheric models and data inversion. This study has provided a rare dataset of the single-scattering albedo in the visible and near-infrared with good spectral resolution. Comparison revealed that our results agree quite well with previous studies. The results show that the single-scattering albedo of this

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
particular Martian dust is constant and high, 0.97, in the shorter wavelengths (1.0 -2.6 μm) and that there is significant absorption by hydrated minerals between 2.7 and 3.6 μm. However, our results do not show the hydration band of phyllosilicates that was seen by [START_REF] Fedorova | Remote sounding of the martian dust from ISO spectroscopy in the 2.7 μm co 2 bands[END_REF]. The discrepancy between the data of [START_REF] Fedorova | Remote sounding of the martian dust from ISO spectroscopy in the 2.7 μm co 2 bands[END_REF] and this study may be, at least partly, explained by the ISO experiment mapping the full disk of Mars and the season being characterized by fairly low dust optical thickness.

We also retrieved temperature profiles from simultaneous observations of the PFS instrument and compared them with profiles acquired with 1-D atmospheric modeling with a dusty atmosphere. The model profiles and the retrieved dust storm profile from PFS reach an acceptable degree of agreement only if the dust storm is confined in the lowest 1-2 scale heights of the atmosphere. The short duration of the storm (less than 3 sols) can support this: the dust has not had the time to mix vertically higher. However, the OMEGA data do not support this hypothesis. In addition, the details of the temperature profiles can not be studied because of limited vertical reliability of the retrieved PFS profile in this situation of very high optical thickness.

We also present a novel approach for dust storm discovery, the Independent Component Analysis. In this paper we have shown that this method can be successfully used to distinguish suspended dust from the surface features, if the dust optical thickness is high enough.

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

A The Independent Component Analysis method

The Independent Component Analysis is a technique that comes from developments in the Blind Source Separation (BSS) research. The goal of ICA is to estimate h source signals, assumed to be stationary and independent, using n observed signals (n ≥ h) that are independent unknown mixing of the source signals [START_REF] Comon | Analyse en composantes indépendantes et identification aveugle[END_REF][START_REF] Comon | Independent component analysis[END_REF][START_REF] Comon | Independent component analysis, a new concept?[END_REF]). The most common example of this problem is the "Cocktail Party" problem. In spectroscopy, it is considered that the observed spectra are linear combinations of pure sources or poles that we want to estimate. For a complex ensemble of data in which more signals contribute to the observed data, the gaussian hypothesis of the Principal Component Analysis (PCA) will not be valid in all cases. Indeed, the goal of the PCA is to maximize the dispersion, and consequently, to decorrelate the signals assuming that all the sources follow a gaussian probability function. The extraction of pure and statistically independent spectra of a melange cannot be reliably obtained with the PCA. However, ICA like PCA, tries to transform the original representative space by searching for directions in a new space, so that the resulting vectors are independent, and not only uncorrelated. The ICA is then a method of linear transformation in which the representation is the one that minimizes the statistical dependence of the components. This is achieved, using a criterion, related to the information entropy theory, that gives back the statistical independence assuming that the data follow a non-Gaussian distribution.

The use of ICA for multispectral data analysis in planetology has been firstly performed by Forni et al. (2005) on OMEGA data proving its efficiency. It has also been applied on the same data set in combination with the Bayesian Positive Source Separation [START_REF] Schmidt | Classification de la surface de mars par imagerie spectrale omega. suivi spatio-temporel et étude des dépôts saisonniers de co 2 et h 2 o[END_REF]. [START_REF] Gomez | Nfindr method versus independent component analysis for lithological identification in hyperspectral imagery[END_REF] have applied ICA on hyperspectral terrestrial data with the fully constrained least squares spectral mixture analysis (FLCS) to estimate the end member abundances.

A.1 Hypotheses and modeling

In a first approximation, it is assumed that the mixing is linearly instantaneous, that the noise can be neglected and that n = h. Be s t and x t , the source and observation vectors, there is a relation such that

x t = As t (1)
whre A is the so called "mixing" matrix, assumed to be constant and invertible.

The goal is to estimate G, the separation or de-mixing matrix, from the only observation matrix x t , such that GA is a diagonal matrix up to a permutation.

The starting hypothesis to solve this kind of problem is the independence of the sources.

In the frame of the linear model X = AS there exists a very significative link between dependence, correlation and non-gaussianity. Two random variables are said independent when the knowledge of the value of one of the variables does not bring any information about the value of the other. Independence is characterised by the fact that the joint probability distribution is equal to the product of the probability distributions. Hence, independence is a stronger hypothesis than uncorrelation. In other words, two variables that are independent are uncorrelated but the opposite is not necessarily true. This means that in the search of independent component, the decorrelation is half of the job and one may seek the basis in which the data is represented by maximally This is done by maximizing a non-Gaussianity criterion based on higher order statitistics.

A so-called "whitening" of the data is first applied so that the mixing matrix is orthogonal and hence the de-mixing matrix as well. In fact, by multiplication with the inverse of a square root of the covariance matrix of the data one obtains:

xt = Ãs t (2)
with à = P D 1/2 P T A where D is a diagonal matrix containing the eigenvalues of the covariance matrix of x t and P is the orthogonal matrix of the associated eigenvectors. This insures that E[x t xt T ] = I is orthogonal, being I the identity matrix and E[X], the expectation value of the random variable X.

The different components of xt are decorrelated and xt is "white" Assuming that the sources s t have zero mean and unit variance, we obtain:

E[x t xT t ] = ÃE[s t s T t ] ÃT = Ã ÃT = I (3)
The mixing matrix à is orthogonal and since G à = I, G is also orthogonal.

This leads to so-called orthogonal algorithms to retrieve the de-mixing matrix G.

After whithening of the data, one is then seeking a rotation R (which leaves things white) so that ŝt defined by ŝt = Ã-1 R xt (4) is an orthogonal ICA method and, like most mainstream ICA techniques, it exploits higher order statistics related to the fourth order cumulant

k(X) = E[X 4 ] -3E[X 2 ] 2 . ( 5
)
This estimator is relatively sensitive to the noise and the number of samples should be high enough to minimize its influence. From the matrix of sources S, it is possible like in the PCA, to compute C the factorial coordinates or "scores" from the following equation:

C = XS T (S.S T ) -1 (6)
Another way to do this, is to compute the correlation coefficient of the sources S with the signal X. The advantage is that, since the sources are independent, the sum of the correlation coefficient squared is 1. This property allows to have an estimate of the contribution of each independent component to the signal.

More details on JADE can be found in [START_REF] Cardoso | Blind signal separation: statistical principles[END_REF][START_REF] Cardoso | High-order contrasts for independent component analysis[END_REF]; [START_REF] Hyvärinen | Independent Component Analysis[END_REF]. The color scale ranges between 0-9.9.

  atmospheric conditions favorable for dust storms.

  Express orbiter was successfully inserted into Martian orbit on December 25th, 2003 (Ls = 322 • ). The satellite is located on an eccentric quasi-polar orbit (inclination of 86.35 • ). The periapsis altitude is approximately 250 km, while the periapsis latitude has drifted during the mission from close to equator at A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT the beginning of the mission to over the South Pole in July 2004 (Orbit #523, Ls = 48.8 • ), and then northward to reach the North Pole in April 2005 (Orbit #1536, Ls = 183.4 • ). The spatial resolution obtained at a given latitude hence depends on the orbit or season of the observations.

Figure

  Figure ??, panels a and b, compares the OMEGA spectra obtained by averaging of pixels common to both orbits #1201 and #1212 between 2 • S and 4 • S. In the solar-reflected part of the spectra, the reflectance factor in the 1 -2.5 μm range rises from ∼ 0.3 for the orbit #1201 to about 0.40 for the

Figure

  Figure ?? shows a map of the two major components of the ICA decompo-

  Figure ?? shows a map of the two major components of the ICA decomposition of orbit #1212: i) the first component (left panel) corresponds to the dust storm in the middle of the map, ii) the second component (right panel)

  The single-scattering albedo was determined for several different dust optical thicknesses at 1 μm. As shown in Fig. ??, for τ = 5 or τ = 20, the retrieved ω 0 showed an artifact, either positive or negative, of the CO 2 gaseous A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT band at 2.0 μm. This unrealistic behavior for dust spectral properties allowed us to determine the dust optical thickness for the averaged spectrum. A single-scattering albedo devoid of CO 2 band residuals was obtained for τ = 10 (Fig. ??), which corresponds to the highest dust optical thickness in the core of the dust storm. [Figure 5]

  -scattering albedo is shown in left panel of Fig.??. In the 0.9 -2.5 μm range, it agrees within 0.5% to the single-scattering albedo retrieved using the same OMEGA data byVincendon et al. (2007). The small difference between the two studies can be accounted for by our choice of a different radiative transfer scheme than that usedVincendon et al. (2007) (two-stream against Monte-Carlo). However, our study expands the analyzed spectral range to the entire OMEGA spectrum dominated by solar reflection (0.4 -4 μm). In this range our retrieved values are consistent in absolute value and in spectral shape with the single-scattering albedo derived by[START_REF] Ockert-Bell | Absorption and scattering properties of the Martian dust in the solar wavelengths[END_REF]. However, our values are larger than the dust parameters retrieved from Pathfinder visible data(Tomasko et al., 1999;[START_REF] Markiewicz | Optical properties of the martian aerosols as derived from imager for mars pathfinder midday sky brightness data[END_REF][START_REF] Johnson | Dust deposition at the mars pathfinder landing site: observations and modeling of visible/nearinfrared spectra[END_REF]. Particularly around 1 μm, this difference cannot be attributed to the 15%-uncertainty of the OMEGA absolute calibration as shown in Fig.??. We do not have any explanation for this disagreement. We just note that the Pathfinder and the OMEGA dataset are not sensitive to the same parameters. The Pathfinder data are mostly sensitive to the particle mean size, while the spectral properties inferred byTomasko et al. (1999)

  i) ISO observations were averaged over the disk of the planet, hence encompassing regions of different dust opacities, ii) ISO observations were carried out in a period of low dust opacity when the dust should have been confined close to the surface, while our observations pertain to a local dust storm efficiently mixing dust in the first atmospheric scale heights. Indeed Fedorova et al. (2002) needed a higher single-scattering albedo at the wavelength of the strongest CO 2 absorption, which points to the dust being confined close to the surface during their observations.

Figure ? ?

 ? Figure ?? presents the latitude-pressure cross-sections of the temperature structures inverted for orbit #1201 and #1212. Except in the local dust storm area, the two cross-sections are quite consistent with each other at pressures larger than 70 Pa. They are also consistent with the predictions of the LMD General Circulation Model (Forget et al., 1999) retrieved from the European Mars Climate Database, EMCD (http://www-mars.lmd.jussieu.fr/, Lewis et al. (1999)) with atmospheric temperatures of 230 K at tropical latitudes in the lower scale height, and of 180 K at the 70-Pa pressure level. At pressures lower than 70 Pa, the two temperature structures differ in the Southern Hemisphere. At 30 • S, the temperature retrieved was ∼ 170 K at 30 Pa during the orbit#1201 while it was over 180 K during the orbit #1212. In fact, the temperature structure retrieved during orbit #1212 was closer to the predictions of the EMCD, with the warming at high altitudes in the Southern Hemisphere corresponding to the return branch of the Hadley circulation. Possibly the difference between the two orbits indicates a temporal variability of the Hadley circulation cell.

  observed the Southern Hemisphere in its late summer season in 2003 (MY26) and 2005 (MY27) using the CO rotational lines. The temperature profiles they derived in 2003 and 2005 differ with each other, the profile obtained in 2005 showing a weaker vertical temperature gradient than in 2003. This corresponded to the break-

  was used to study the effect of dust on the atmospheric temperatures: the 1-D model version includes all the physical process descriptions as the global model, and an improved, more accurate radiative transfer scheme to compute the absorption and scattering of solar radiation by airborne dust, based on the Toon et al. (1989) scheme. The dust optical properties given as input for the model were those used and retrieved in this work. The model was initialized with an isothermal temperature profile with T = 200 K at all model layers and the results were taken from the sol 10 at 13.30 local time (the time of the observation), if not stated otherwise. The model-produced profiles were convolved with PFS weighting functions to optimize the comparison with the PFS temperature profile retrieved in the storm area (PFS FOV of the temperature retrieval is shown in Fig. ??).

  ??): already the 01 LT profile of the first sol displays the warm, isothermal middle atmosphere that was seen at sol 10 in Fig. ??. So the dust heats the atmosphere very fast and it seems very unlikely that the storm would have started so recently before the observation that the atmosphere would not yet have responded the observed and modelled temperatures could rather indicate that dust vertical distribution in the dust storm observed by OMEGA is not uniform throughout the atmosphere, but that the dust is confined in the lowest levels. A new set of model runs was made with a modified dust distribution of a well-mixed τ = 0.2, on top of which a τ = 4.23 profile limited in vertical extent was embedded. The altitude of the top of the dust storm was varied and the produced profiles were compared with the retrieved PFS profile. Examples of the model profiles are shown in Fig. ??. The dust top altitude has a notable effect on the shape of the temperature profile (see Fig. ??). All the LMD 1-D profiles with vertically heterogeneous

  those bases in which the data is decorrelated.

  A -1 = Ã-1 R are estimations of the sources and of the inverse of the mixing matrix. A.2 JADE As in Forni et al. (2005) and Erard et al. (2008), we have used the Joint Approximate Diagonalization of Eigenmatrices method (JADE) algorithm, which

Fig. 1 .

 1 Fig. 1. Tracks of orbit#1201 session 3 and orbit#1212 session 3 mapped on cylindrical projection of the Martian surface. The dust storm studied in this paper was seen on orbit #1212 between the latitudes shown by the dashed lines. The PFS FOV of orbit #1212 for the temperature profile inversion is shown with the oval. The local time of observations was about 13.30.

Fig. 3 .

 3 Fig. 3. Comparison of the PFS spectra averaged over the range 2 • S -4 • S obtained during orbit#1201 (dotted line) and orbit#1212 (solid line).

Fig. 4 .

 4 Fig. 4. Spatial maps of the two major components of the ICA decomposition of orbit #1212. Left panel: the dust component. Right panel: the surface aldedo component.

Fig. 5 .

 5 Fig. 5. The behavior of the single-scattering albedo ω 0 with 3 different values of dust optical thickness τ = 5, 10, 20. The value τ = 10 with the flattest behavior of ω 0 was used for the retrieval.

Fig. 6 .

 6 Fig. 6. Left: The retrieved single-scattering albedo (solid line) and comparison with previous studies of Fedorova et al. (2002) (dashed line), Ockert-Bell et al. (1997) (stars), and Tomasko et al. (1999) (crosses). Right: Contours of the I/F as a function of single-scattering albedo ω 0 and asymmetry parameter g for the infinite optical depth case. The star represents the ω 0 derived from our OMEGA data at 1.07 μm.The error bar shows the uncertainty on ω 0 due to the 15% uncertainty on OMEGA absolute calibration.

Fig. 7 .

 7 Fig. 7. The dust optical thickness at 1.0 μm in the area of the storm (0.4 -4.7 • S).

Fig. 11 .

 11 Fig. 11. Temperature profiles from the first sol of the model run with τ = 4.23 and dust top altitude 55 km at local times 00, 08, 13 and 20. The model was started at 00 LT and initialised with an isothermal temperature profile.

  Forget, F., Hourdin, F., Fournier, R., Hourdin, C., Talagrand, O., Collins, M.,

	Lewis, S. R., Read, P. L., Huot, J.-P., 1999. Improved general circulation
	models of the martian atmosphere from the surface to above 80 km. J.
	Geophys. Res. 104, 24155-24176.
	Forget, F., Spiga, A., Dolla, B., Vinatier, S., Melchiorri, R., Drossart, P.,
	Gendrin, A., Bibring, J.-P., Langevin, Y., Gondet, B., 2007. Remote
	sensing of surface pressure on Mars with the Mars Express/OMEGA
	spectrometer: 1. Retrieval method. J. Geophys. Res. 112, E08S15,
	doi:10.1029/2006JE002871.

Acknowledgements

We thank the reviewers for comments that significantly improved the paper. We thank Dr. A. Fedorova for providing the data for comparison and A. Spiga for discussions about the surface pressure values and providing the program for plotting nice maps of OMEGA images. AM wishes to thank Association Franco-Finlandaise pour Recherche Scientifique et Technique, the Wihuri foundation, the Enhrnrooth foundation and the Väisälä foundation for their financial support that made possible the visits to Paris (approximately 5 months in total) to do this work during the period Nov 2005 -Dec 2007. Funding from the Kordelin foundation, the Graduate school in astronomy and space physics and the Centre National d'Etudes Spatiales are gratefully acknowledged.

Dust storm seen by MEx/OMEGA and MEx/PFS


F Brissaud

C Smidts

A Barros

C Bérenguer

Dynamic Reliability Modeling of Cooperating Digital-Based Systems

Dynamic reliability explicitly handles the interactions between the stochastic behavior of system components and the deterministic behavior of process variables. However, its industrial level applications are still limited, notably due to the inherent complexity of the theory and the lack of a generic modeling framework. The increased use of digital-based systems has also introduced additional modeling challenges related to the interactions between cooperating digital components. For solving these challenges, the present paper first extends the mathematical framework of dynamic reliability to handle 1) information and data computed and exchanged between digital components; and 2) random parameter deviations. A formalized Petri net approach is then proposed to perform the corresponding reliability analyses, using a finite element method. Finally, the framework's effectiveness is demonstrated on a simplified model of a nuclear reactor case study.

INTRODUCTION

Probabilistic risk (or safety) assessments (PRA) provide a general framework for managing risks linked to engineering systems, notably in nuclear power plants, aerospace, and chemical industries. The main purpose of PRA is to identify the possible accidental scenarios, to rate their consequences, and to assess their likelihoods. An essential challenge is to deal with the system complexity, that is, to treat system interactions both efficiently and realistically.

In the seventies, the event tree/fault tree methodology was introduced for PRA (US NRC 1975), focusing on system components and their "static" relationships. In the late eighties-early nineties, dynamic reliability (or probabilistic dynamics) methods were developed to explicitly handle the influence of time, process dynamics and human operations, on system failures and PRA scenarios. These approaches include dynamic discrete-time event trees (e.g. DY-LAM, DETAM), the GO-FLOW methodology, Event Sequence Diagrams (ESD), the Dynamic Flowgraph Methodology (DFM), Markov models, and Petri Nets. Most of them have been presented, compared and discussed by Aldemir et al. (1994), [START_REF] Siu | Risk assessment for dynamic systems: an overview[END_REF], and [START_REF] Labeau | Dynamic reliability: towards an integrated platform for probabilistic risk assessment[END_REF]. Although there is a broad consensus on the need for dynamic reliability methods [START_REF] Aldemir | Special issue on reliability and safety analysis of dynamics process systems[END_REF], such approaches have not yet penetrated the arena of industrial applications [START_REF] Labeau | Dynamic reliability: towards an integrated platform for probabilistic risk assessment[END_REF]. Given reasons are the theoretical flavor and the lack of a generic platform for performing such analyses.

The formal mathematical framework of dynamic reliability was established under the name of Continuous Event Tree (CET) theory, (Devooght & Smidts 1992a). The latter introduced two sets of variables to define the complete system state:  physical (continuous) variables of the process (e.g. level, pressure, temperature);  state (discrete) variables of the system components (i.e. system configuration). The evolution of the process variables was characterized by a set of first order (non-stochastic) differential equations for each value of state variables. Changes in state variables were assumed stochastic, defined by transition rates, and dependent on process variables. Human operations were also included as additional variables (Devooght & Smidts 1992b).

More recently, the use of digital safety-related systems has introduced new issues for dynamic reliability modeling, especially due to the interactions between system components [START_REF] Aldemir | Current State of Reliability Modeling Methodologies for Digital Systems and Their Acceptance Criteria for Nuclear Power Plant Assessments[END_REF]. Typical examples of such systems include the "intelligent transmitters" which are able to exchange information, to perform internal data processing and advanced functionalities (Brissaud et al., sub.).

The present paper focuses on the dynamic reliability modeling of cooperating digital-based systems. First, the mathematical framework of dynamic reliability is extended in Section 2 to deal with the specific characteristics of these systems. A formalized Petri-net approach is also proposed to perform the related reliability analyses. Finally, a case study involving a nuclear reactor is presented in Section 3.

DYNAMIC RELIABILITY FRAMEWORK FOR

COOPERATING DIGITAL-BASED SYSTEMS

Problem Formulation

In addition to time t, four types of variables are used to describe the complete system state. The process and components state variables remain the same as those defined by Devooght & Smidts (1992a). Data variables are introduced in order to characterize digital systems; and deviation variables extend the possibilities of failure modeling. All these variables are time-dependent. The process, data, and deviation variables are continuous and depicted by vectors of reals, respectively denoted x(t), y(t), and e(t). Components state variables are discrete and depicted by a vector of integers, denoted i(t). (cf. Table 1).

The process variables x(t) represent the physical variables involved in the system dynamics (e.g. pressure, temperature, volume). They evolve deterministically, given components state, and with deviations as parameters (e.g. the level in a tank is determined by the configuration of the valves, and the amount of leakage). The evolution of process variables may then usually be defined by a set of first order differential equations, indexed by the components state:

dt d x(t) = x' i (x, e, t) (1)
The data variables y(t) represent any information or data which are computed, stored, and/or exchanged between system components (e.g. commands, measurement results, diagnostic information). By nature, they do not directly affect process and deviation variables, but are used to change components state. The data variables may usually be expressed as a function of process and deviation variables, given components state (e.g. when a transmitter is in a non-fully operating mode i.e. a component state, its measurement results depend on the quantities to be measured i.e. process variables, and drifts i.e. deviation variables).

The data variables may also depend on their previous values (e.g. stored data, locked signals), denoted -y(t), with -y(t) = y(tε) and ε which tends to 0 + :

y(t) = y i (x, -y, e, t) (2)
The deviation variables e(t) represent continuous errors or deviations in system properties, which evolve stochastically (e.g. system degradations, drifts), depending on process variables and components state (e.g. the leak in a closed valve follows a random distribution influenced by the flow rate). Because the evolution of deviation variables is continuous, it may usually be defined by a set of first order differential equations which include random variables (e.g. the rate of crack growth is a random variable which depends on the current crack level), indexed by the components states:

where E' i (x, e, t) is a (function of) random variable.

The components state variables i(t) represent the structure (configuration) of the system and is a function of the states of its components (operating or failed) and of human operations (e.g. opening or closing a valve). The state of any system component (e.g. operational, degraded, or failed) can be described by integers which are arranged in vector i(t). The components state variables may evolve both deterministically and stochastically, depending on process and data variables (e.g. a valve is controlled by a signal; a transmitter failure rate depends on the temperature), and deviations (e.g. after a certain level of degradation, a component transition occurs from a degraded mode into a fully non-operating mode). The components transition rate from state i k to state i l at time t, given process, data, and deviation variables, is denoted p(i k → i l | x, y, e, t). Then, the total components transition rate from state i k is:

λ ik (x, y, e, t) =   k l i i p(i k → i l | x, y, e, t) (4)
The transitions between components states are assumed instantaneous. When the components state at time t is i(t) = i k , the probability that the components leave this state before time t + τ is therefore:

F ik (τ | x,y,e,t) = 1 -exp        0 λ ik (x,y,e,t+u)•du     (5)

Mathematical Solution Using a Finite Element Method

For the numerical analyses of the complete system evolution according to time, a finite element method is adopted. A time step Δt is used which should be small enough to assume that variables x(t), y(t), e(t), and i(t) are constant in any time interval [t, t + Δt[ without loss of accuracy. These variables at time t + Δt can then be determined by their values at time t. In particular, a components state transition which occurs between time t and time t + Δt is considered to occur exactly at time t + Δt. In the same way, the evolution of the process, data, and deviation variables between time t and time t + Δt are considered to occur as "jumps" exactly at time t + Δt. It is then possible to approximate the values of process and deviation variables at time t + Δt, according to the complete system state at time t, using the finite differences of the derivatives given in Section 2.1:

x(t + Δt) ≈ x(t) + Δt • x' i (x, e, t) (6) e(t + Δt) ≈ e(t) + Δt • E' i (x, e, t) (7)
In addition, the probability that the components remain in their current state at time t, denoted i(t

) = i k , up to time t + Δt, that is i(t + Δt) = i k , is: Pr[i(t + Δt) = i k | i(t) = i k , x, y, e, t] = 1 -F ik (Δt | x, y, e, t) ≈ 1 -Δt • λ ik (x, y, e, t) (8)
And similarly, the probability that the components leave their current state at time t, denoted i(t) = i k , for another specific state at time t + Δt, denoted i(t +Δt) = i l , with i k ≠ i l , can be approximated by:

Pr[i(t + Δt) = i l ≠ i k | i(t) = i k , x, y, e, t] ≈ Δt • p(i k → i l | x, y, e, t) (9)
Note that a deterministic (certain) transition can then be modeled using a rate equal to 1/Δt. The couple (Ï, P i (x, y, e, t)) is defined, with Ï a vector composed of all possible combinations of components states, and P i (x, y, e, t) a vector with all its components equal to 0 except one that is equal to 1 whose components are determined randomly according to Equations ( 8) and ( 9), in such a way that:

i(t + Δt) ≈ P i (x, y, e, t) T • Ï (10)
Once the components state, process and deviation variables are defined at time t + Δt, the data variables at time t + Δt can also be determined, using ε = Δt which implies that -y(t + Δt) = y(t):

y(t + Δt) ≈ y i (x, -y, e, t + Δt) (11)
Equations ( 6)-( 11) show that the complete state of the system at time t + Δt, i.e. (i, x, y, e, t + Δt), can be fully determined according to its state at time t, i.e. (i, x, y, e, t), according to deterministic and stochastic evolutions. The system is therefore a piecewise-deterministic process (PDP), [START_REF] Davis | Markov models and optimization[END_REF]).

Petri Net Formalism for Numerical Analyses

Petri nets and their extensions, including stochastic and colored characteristics, provide natural and effective tools for modeling dynamic systems [START_REF] David | Petri nets for modeling of dynamic systems: A survey[END_REF], notably for risk analysis [START_REF] Vernez | Perspectives in the use of coloured Petri nets for risk analysis and accident modelling[END_REF]. Stochastic Petri nets were also used efficiently in dynamic reliability [START_REF] Dutuit | Dependability modelling and evaluation by using stochastic Petri nets: Application to two test cases[END_REF]). In the present paper, a Petri net formalism, using stochastic and colored properties, is proposed in order to provide a generic framework to:  flexibly model the dynamic reliability of a system with the help of a visual interface easy to handle;  simulate the evolution of the complete system state, using a finite element method. In the proposed approach, each place of the Petri net is associated to one set of variables, and vice-versa. The number of places then increases linearly according to the number of variables, which avoids any combinatorial explosion. According to the nature of the variables (continuous or discrete, stochastic or deterministic), different representations are used for graphical convenience (cf. Figure 1). The values of the variables are given by the (colored) token, with real or integer numbers (according to the places) inside the corresponding places, and are changed by the transitions. Each place therefore always contains one and only one token, and thus, all the transitions are always enabled. Guards are then used for each transition and denoted s j [Δt], which means that the transition is fired at each time instant

s j + k • Δt, with k = 0, 1, 2,…
Each transition of the Petri net is associated to a place, denoted "managed place." This place is linked to the transition by an input arc, which means that the corresponding variables are changed by the transition ("the token is removed from the place"); and linked to the same transition by an output arc, which attributes to the variables their new values ("the token is deposited in the place"). An expression is given on this output arc to specify these new values, which may depend on the previous values of the variables (handled by the input arc) and also on variables from other places. The latter places are then denoted "dependence places" and are linked to the corresponding transition by bi-directed arcs, which means that the values of their variables may be used by the transition, but are not changed.

The new values of variables (specified by the output arcs from the transitions to their "managed places") may also depend on random variables. Contrarily to the "classical" stochastic Petri nets, the stochastic aspects therefore do not relate to time instants, but to values of variables (i.e. the token "color"). In this approach, the time, as a variable, is modeled by a place. Besides, all the transitions are fired at deterministic time instants (specified by the transition guards). This approach can therefore be classified as an "untimed stochastic (and colored) Petri net." A generic Petri net for the dynamic reliability modeling of systems is depicted in Figure 2, using the elements described in Figure 1. The five types of variables defined in Section 2.1 are modeled, and depicted using different varieties of places. In Figure 2, the variables are represented by vectors (except for time t). For more detailed models, it is also possible to split these vectors into subsets with one place for each (subset), and to treat them separately.

Each transition is fired at every time step Δt, changing the variables modeled by the corresponding "managed place," according to the equations given in Section 2.2 and specified on the output arcs. The transition that changes the time variable t is not linked to any "dependence place" and is simply used to increment time t by Δt at each solicitation. On the other hand, the components state variables i(t) are changed by random variables defined by Equation (10), which depend on all the variables. In each time interval [t, t + Δt], all the values of variables at time t + Δt are computed following a specific order defined by the s j (cf. top of Figure 2). (Such specific orders may also be required to deal with dependencies between subsets of variables.) Note that the values of variables at time t + Δt depend on values at time t (cf. Equations ( 6)-( 11)). In particular, x(t + Δt) and e(t + Δt) both depend on x(t) and e(t).

To avoid "losing" the values of x(t) (resp. e(t)) after computing x(t + Δt) (resp. e(t + Δt)), "Metatransitions" are introduced. They are used to first compute all the derivatives at time t (i.e. x' i (x, e, t) and E' i (x, e, t)), storing them as additional variables (cf. Figure 3), and then to change the variables of the complete system state. A "Meta-transition" has therefore a double guard denoted (s j , s k ) [Δt], which means that the derivative is computed at each time instant s j + k • Δt, and the variables of the "managed place" at s k + k • Δt.

CASE STUDY OF A FAST REACTOR

The Europa Fast Reactor

The case study is the primary circuit of the Europa fast reactor. This system has been proposed as a benchmark problem on accident sequences [START_REF] Wider | Comparative analysis of a hypothetical loss-of-flow accident in an irradiated LMFBR core using different computer models for a common benchmark problem[END_REF]). Several dynamic reliability analyses have been also performed on this application [START_REF] Amendola | DYLAM-1: A Software Package for Event Sequence and Consequence Spectrum Methodology[END_REF], Smidts & Devooght 1992[START_REF] Swaminathan | An application of the ESD framework to the probabilistic risk assessment of dynamic systems[END_REF]. A comprehensive description of the system is provided by Smidts & Devooght (1992). In the present paper, the physical variables have been simplified to allow us to focus more specifically on other aspects. In particular, new transmitter features have been introduced (communications and drift corrections), and deviation variables.

The simplified model of the primary circuit of the Europa fast reactor is depicted in Figure 4. This system is comprised of two channels (C1 and C2) where sodium is introduced as a coolant by a pump (PM). The lack of coolant, for example in case of a pump failure, increases the temperature in the channels and may yield hazardous events. The sodium temperatures in the channels (T 1 and T 2 ) are therefore monitored by transmitters (ST1 and ST2) which send their results (T1 s and T2 s ) to a common controller (CT). Similarly, the flow rate (G) is monitored by a third transmitter (SG) which sends its result (G s ) to another controller (CG). CT and CG send their signals (y CT and y CG ) to a central controller (SDL). If a high temperature threshold (T max ) or a low flow rate threshold (G min ) is detected, then SDL sends a signal (y SDL ) which should activate an emergency shutdown (SCM). This safety device, commonly named SCRAM, consists in inserting (under gravity) control rods into the core which quickly stop the nuclear reaction by absorbing neutrons. The system variables are described in Sections 3.2 to 3.5.

Components State Variables

The state of each of the eight system components is represented using a finite integer and constitutes one of the components of vector i(t), i.

e. i(t) = (SPM(t), SSG(t), SCG(t), SST1(t), SST2(t), SCT(t), SSDL(t), SSCM(t)).

In the following, and in accordance with the formalism presented in Section 2, each component of vector i(t) is modeled separately. The components state variables are then defined in Table 2.

The normal (full) operating modes of the components are represented by state variables equal to 1. When a component state variable is equal to 0, the corresponding failure mode is "dangerous," that is, it may directly yield an inability of the system to perform its safety function (i.e. inserting the control rods into the core). On the contrary, when a component state variable is equal to 2, the corresponding failure mode is "safe," that is, it may directly yield a spurious activation of the safety function. Other values of component state variables (3 or 4) correspond, for example, to "degraded" modes of operation.

The state variables of the "mechanical" components, that is, the pump and the SCRAM, directly affect process variables as described in Section 3.3. On the other hand, the state variables of the controllers and transmitters directly determine data variables defined in Section 3.4. Finally, the effects of the "degraded" modes of operations are modeled using deviations variables defined in Section 3.5.

The transition rates between the possible states of each component are given in Table 3. Note that transition rates depend on time t, process variables, deviation variables, and states of other components. A deterministic (certain) transition is also assumed for the SCRAM activation (using a rate equal to 1/Δt).

(t) 1 or 3 0 1•10 -3 • exp(δ M (t) • 5•10 -5 ) 1 3 1•10 -2 SSG(t) 1 0 2•10 -3 1 2 2•10 -4 SCG(t) 1 0 or 2 1•10 -5 SST1(t) 1, 3 or 4 0 4•10 -4 •(1+Ι 1 (SST2(t)=0))•r(T 1 (t)) 1, 3 or 4 2 4•10 -5 •(1+Ι 1 (SST2(t)=2))•r(T 1 (t)) 1 3 or 4 1.5•10 -2 SST2(t) 1, 3 or 4 0 4•10 -4 •(1+Ι 1 (SST1(t)=0))•r(T 2 (t)) 1, 3 or 4 2 4•10 -5 •(1+Ι 1 (SST1(t)=2))•r(T 2 (t)) 1 3 or 4 1.5•10 -2 SCT(t) 1 0 or 2 1•10 -5 SSDL(t) 1 0 or 2 1•10 -6 SSCM(t) 1 0 3.125•10 -8 • t • r((T 1 (t)+T 2 (t))/2) 1 5 (1/Δt) • Ι 1 (y SDL (t)=1) __________________________________________________ * T 1 (t)
and T 2 (t) are process variables defined in Section 3.3, δ M (t) is a deviation variable defined in Section 3.5, Ι 1 (.) and r(.) are functions defined by Equations ( 12) and ( 13).

The following function is used to model the effect of the temperature process variables (T 1 and T 2 ) on transition rates (cf. Table 3):

r(T) = 6.17•10 -2 • exp(5.21•10 -3 • T) (12)
The following function is used to model dependencies with component states (e.g. between the temperature transmitters, cf. Table 3):

Ι 1 (A) = 1 if assertion A is true

Process Variables

The thirteen process variables are x(t) = (ω(t), G(t), T c,1 (t), T c,2 (t), T 1 (t), T 2 (t), P(t), C 1 (t), C 2 (t), C 3 (t), C 4 (t), C 5 (t), C 6 (t)) and defined hereafter. Expressions of the derivatives are given by the following equations, with the initial conditions at time t 0 = 0. The parameters and their values are reported in Table 4. Angular speed of the pump, denoted ω(t) [rad/s]:

                    s rad K C I t K t SPM t C t dt d M M M / 6000 ) 0 ( ) ( 0 ) ( ) ( ) ( 1     (14)
where δ M (t) is the pump torque deviation defined in Section 3.5, and SPM(t) is the pump state variable defined in Section 3.2 (cf. Table 2).

Sodium linear momentum flow rate, denoted G(t) [kg/m²s]:

               s m kg G G t C G t G t G dt d nt ² / 4000 ) 0 ( ) 0 ( ) 0 ( ) ( ) 0 ( ) ( ) ( 2 2 2 2     (15) Fuel temperature in channel i, denoted T c,i (t) [K], with i = 1, 2:                           2 1773 1 1437 ) 0 ( ) 0 ( ) ( ) ( ) ( ) ( , , , , i if K i if K P w R T T t T T t T t P w R t T dt d i e i c i c i e i c i i c  (16)
                                               K G t T C A P w T T T t T t T t G dt d t G t G t T C t T A t P w t T dt d i R i i e i e i i c i i R i c i i i i 873 ) 0 ( ) ( 2 ) 0 ( ) 0 ( ) ( ) ( 1 ) ( ) ( 1 ) ( ) ( ) ( 2 ) ( ) ( , ,   (17)
Power generated by the core, denoted P(t) [W]:

                                   MW P t C t P t SSCM t t P dt d i i i i i SCM 2000 ) 0 ( ) ( ) ( 5 ) ( ) ( 6 1 6 1 1    (18)
where SSCM(t) is the SCRAM state variable defined in Section 3.2 (cf. Table 2).

Precursor concentration, denoted C i (t) [W], with i = 1, 2, 3, 4, 5, 6:

                 ) 0 ( ) 0 ( ) ( ) ( ) ( P C t P t C t C dt d i i i i i i i     (19)
Note that the introduction of function Ι 1 (A) allows the use of only one set of differential equations, which depend on components state variables as parameters. These equations also depend on deviation variables. Finally, Equations ( 14)-( 19) are appropriate for the analyses presented in the present paper, but are not applicable to any other cases.

To illustrate the evolution of the flow rate G(t) and temperatures T 1 (t) and T 2 (t), three scenarios are assumed and depicted in Figures 5 and6:  s1: SPM(t) = 3, and SSCM(t) = 0 for any time t;  s2: SPM(t) = 3 for any time t, and SSCM(t) = 1 up to time t = 53 s, then SSCM(t) = 5;  s3: SPM(t)=3 up to time t = 30 s, then SPM(t)=0, and SSCM(t) = 1 up to time t = 33 s, then SSCM(t) = 5.

In addition, the pump torque deviation δ M (t) has been simulated according to Equation (26), using the same results for all these scenarios. Processing data, denoted α i (t) with i = 1, 2, are computed by the temperature transmitters according to other data variables (using communications between transmitters), in order to deduce characteristics on measurement results (cf. Table 5). These data are based on the derivatives of measurement results, compared to theoretical values (cf. Section 3.3):

i (t) __________________________________________________ T i s (t) = T i (t) & G s (t) = G(t) α i (t) = 0 T i s (t) < T i (t) or G s (t) > G(t) α i (t) < 0 T i s (t) > T i (t) or G s (t) < G(t) α i (t) > 0 __________________________________________________ 3.4
(t) [K], with i = 1, 2:       2 ) ( 0 ) ( ) ( 1 ) ( ) ( ) ( 1 min 1 1               t SSG G t SSG t t G t SSG t G t G s s (20)             2 ) ( 0 ) ( ) ( 4 ) ( ) ( 3 ) ( ) ( 4 , 3 , 1 ) ( ) ( ) ( ) ( 1 max 1 1 1 1                        t SSTi
        e s i i c i s s s s s i R i c i i i s i s i i T t T T t t G t G t t G t G t T C T A P w t t T t t T t                             ) ( ) 0 ( 1 ) ( ) ( ) ( ) ( ) ( ) 0 ( 2 ) 0 ( ) ( ) ( ) ( , ,    (25)
According to α i (t), parameters of drifts correction, denoted T i c (t), with i = 1, 2, are then defined according to the rules given in Table 6, assuming T 1 (t) = T 2 (t) before SCRAM activation (cf. Figure 6).

Figure 1 .

 1 Figure 1. Petri net tool box.

Figure 3 .

 3 Figure 3. Meta-transition for Petri net.

Figure 2 .

 2 Figure 2. Generic Petri net for the dynamic reliability modeling of systems.

Figure 4 .

 4 Figure 4. The Europa Fast Reactor

  subject to negative drifts = 4 results subject to positive drifts temperature SCT(t) = 1 correct signals controller = 0 signals locked to "unsafe" value (CT) = 2 signals locked to "safe" value central SSDL(t) = 1 correct signals controller = 0 signals locked to "unsafe" value (SDL) = 2 signals locked to "safe" value SCRAM SSCM(t) = 1 normal operation (SCM) = 0 full failure = 5 SCRAM activation __________________________________________________

Figure 5 .Figure 6 .

 56 Figure 5. Flow rate evolution according to time t (G(t)), in scenarios s1, s2 (equivalent to s1), and s3, cf. Section 3.3

  .m²/s constant for pump friction I = 10 kg.m² pump moment of inertia υ = -5.00•10 -2 s -1 constant for flow rate evolution C nt2 = 5.56•10 -6 kg/m² constant for flow rate evolution R = 9.521•10 -7 K/W thermal resistance w 1 = 0.41175 proportion of power generated by C1 w 2 = 0.58825 proportion of power generated by C2 T e = 653 K sodium reactor inlet temperature τ 1(T) = 1.1223+1.5215•10 -3 •T-1.0471•10 -6 •T 2 +2.7476•10 -10 •T 3 τ 2 (T) = 1.5714+2.1303•10 -3 •T-1.4661•10 -6 •T 2 +3.8470•10 -10 •T 3 A 1 = 0.518867 m² section of sodium passage in C1 A 2 = 0.741238 m² section of sodium passage in C2 C R (T)= 1629-8.3290•10 -1 •T sodium specific heat γ SCM = 4.40•10 -1 s -1 reactivity induced by control rods β 1 = 8.2100•10 -5 delayed neutron fraction β 2 = 7.4480•10 -4 delayed neutron fraction β 3 = 6.6150•10 -4 delayed neutron fraction β 4 = 1.3277•10 -3 delayed neutron fraction β 5 = 6.1480•10 -4 delayed neutron fraction β 6 = 1.8940•10 -4 delayed neutron fraction Γ = 3.98•10 -2 s mean neutron production lifetime γ 1 = 1.29•10 -2 s -1 precursor decay constant γ 2 = 3.11•10 -2 s -1 precursor decay constant γ 3 = 1.34•10 -1 s -1 precursor decay constant γ 4 = 3.31•10 -1 s -1 precursor decay constant γ channel i, denoted T i (t) [K], with i = 1, 2:

Table 1

 1

	. Nomenclature __________________________________________________
	Variable __________________________________________________ Description
	t	time
	x(t)	vector of process variables at time t
	y(t) -y(t)	vector of data variables at time t vector of previous values of y(t) up to time t
	e(t)	vector of deviation variables at time t
	i(t)	vector of components state variables at time t
	x' i (x, e, t)	expression of derivatives of process variables
		at time t, given i

(t), x(t), and e(t) y i (x, -y, e, t

	)	expression of data variables at time t, given
		i(t

), x(t), -y(t), and e(t) E' i ( x, e, t

	)	vector of random variable which provides the
		rate of change of deviation variables at time t,
	i k p(i k →i l |	given i(t), x(t), and e(t) specific value of vector i(t), indexed by k

x,y,e,t) components transition rate from state i k to state i l at time t, given x(t), y(t), and e(t) λ ik (x, y, e, t

	)	total components transition rate from state i k
		at time t, given x

(t), y(t), and e(t) F ik (τ | x, y, e, t

	) probability of leaving components state i k in
	time interval [t, t + τ[
	P

i (x, y, e, t) T • Ï product of vectors which determine randomly i(t+Δt) according to transition rates (i, x, y, e, t) notation for (i(t), x(t), y(t), e(t)), description

	of the complete system state at time t
	Δt __________________________________________________ time step
	dt d e

(t) = E' i (x, e, t) (

Table 2 .

 2 Components

	state variables __________________________________________________
	System	State	Possible value with description
	component variable __________________________________________________
	pump	SPM(t)	= 1 normal operation
	(PM)

Table 3

 3

	. Transition rates between components states __________________________________________________ State From To Transition rate* [s -1 ]
	variable state __________________________________________________ state
	SPM

Table 5

 5

	. Relationships between processing data α i (t), and measurement results T i s (t) and G s (t), with i = 1, 2 __________________________________________________
	Measurement results	α

  Di (t) are drifts defined in Section 3.5. Note that data variables at time t depend on their values at time t -Δt (cf. Sections 2). Controller signals, denoted y CG (t), y CT (t), y SDL (t):

	T i	s	t					T i	t	T i	c	t						SSTi	t	or
						Di	t		SSTi	t						Di	t	SSTi	t	(21)
					T i	s	t	t				SSTi	t	T
	y	where δ   ( 1   t SCG ) ( 1   G t s ) (  ) 2   G t CG	min	  1   SCG	( t	)		 1	(22)
	y	 CT	1 (  t	 SCT ) 1  	 ( t T ) 1 s	  1 1   SCT ) max  T ( t 2 or T t  ( s ) (  ) t	2   T max		(23)
	y	 SDL	( 1  t	 SSDL  ( ) 1   t y CG )	(  t	2 )	 	1	or	y	CT	( t	)		  1 1   SSDL	( t	)		 1	(24)

Table 6

 6

	. Computation of T i __________________________________________________ c (t) according to α i (t), with i = 1, 2 Criteria* on α i (t) Computation of T i c (t) __________________________________________________ |α 1 (t)| > α ref & |α 2 (t)| < α ref T 1 c (t) = T 2 s (t) -T 1 s (t) & T 2 c (t) = T 2 s (t -Δt) |α 1 (t)| < α ref & |α 2 (t)| > α ref T 2 c (t) = T 1 s (t) -T 2 s (t) & T 1 c (t) = T 1 s (t -Δt) α 1 (t) < -α ref & α 2 (t) > α ref T 1 c (t) = ((T 2 s (t) -T 1 s (t))/2) or α 1 (t) > α ref & α 1 (t) < -α ref & T 2 c (t) = ((T 1 s (t) -T 2 s (t))/2)
	otherwise	T i

c (t) = T i s (t -Δt) with i=1,2 __________________________________________________ * For the basic case: α ref = 1

Deviation Variables

The three deviation variables are e(t) = (δ M (t), δ D1 (t), δ D2 (t)) and are defined by the following equations.

The pump torque deviation, denoted δ M (t) [N.m]:

with E M (µ, σ) a random variable which follows a Log-Normal distribution of parameters µ = 5, σ = 1.

The drifts, denoted δ Di (t) [K], with i = 1, 2:

with E D (δ Di (t), T i (t)) a random variable which follows an Exponential distribution of parameter 3•10 4 /(T i (t) • δ Di (t) + 20), with i = 1, 2.

Results

The threshold values used for the analyses are T max = 923 K and G min = 3345 kg/m²s. At time t 0 = 0, all the system components are in the full operating modes (the components state variables are equal to 1), except for the pump which is in the "degraded" mode of operation (SPM(t 0 ) = 3); the process variables are in the initial conditions defined in Section 3.3; the initial values of data variables are defined by Equations ( 20)-( 25); and deviation variables are nil.

The system has been modeled using the formalism presented in Section 2, and its evolution has been simulated using CPN tools [START_REF] Jensen | Coloured Petri Nets and CPN Tools for modelling and validation of concurrent systems[END_REF]), a Petri net software which provides all the characteristics required for the proposed modeling approach (including the color properties, the specification of random functions, and the definition of "metatransition"). The scenarios of system evolution are then classified according to the temperatures in the reactor channels. If the SCRAM is activated while both T 1 and T 2 are lower than T max -5 = 918 K, then a "spurious trip" is assumed. On the other hand, if either T 1 or T 2 exceeds T max + 10 = 933 K for more than 5 s, then a "hazardous event" is assumed. In the other cases, the system is "under control."

The analyses have been performed by Monte Carlo simulations (cf. Equation ( 10) and Section 2.3), and the results are reported in Table 7. To assess the effects of the "intelligent transmitter" features (the possibility to exchange information, to process data, and to offset drifts), a case without drift correction (i.e. α ref = ∞) is also presented. It is then shown that the use of such "intelligent transmitters" allows a greater percentage of "under control" cases, obtained by a reduction of the number of "spurious trips." However, a drawback of decreasing the probability of "spurious trips" is a slight increase of the percentage of "hazardous events." It is therefore deduced that adjustments of parameter α ref should be used to balance reliability and safety.