D Clarke

L Zaccarelli

N M Shapiro

F Brenguier

Assessment of resolution and accuracy of the Moving Window Cross Spectral technique for monitoring crustal temporal variations using ambient seismic noise

Keywords:

Temporal variations in the elastic behavior of the Earth's crust can be monitored through the analysis of the Earth's seismic response and its evolution with time. This kind of analysis is particularly interesting when combined with the reconstruction of seismic Green's functions from the cross-correlation of ambient seismic noise, which circumvents the limitations imposed by a dependence on the occurrence of seismic events. In fact, because seismic noise is recorded continuously and does not depend on earthquake sources, these cross-correlation functions can be considered analogously to records from continuously repeating doublet sources placed at each station, and can be used to extract observations of variations in seismic velocities. These variations, however, are typically very small: of the order of 0.1%. Such accuracy can be only achieved through the analysis of the full reconstructed waveforms, including later scattered arrivals. We focus on the method known as Moving-Window Cross-Spectral analysis that has the advantage of operating in the frequency domain, where the bandwidth of coherent signal in the correlation function can be clearly defined. We investigate the sensitivity of this method by applying it to microseismic noise cross-correlations which have been perturbed by small synthetic velocity variations and which have been randomly contaminated. We propose threshold signal to noise ratios above which these perturbations can be reliably observed. Such values are a proxy for cross-correlation convergence, and so can be used as a guideline when determining the length of microseismic noise records that are required before they can be used for monitoring with the moving-window cross-spectral technique.

INTRODUCTION

Stress field variations in time modify the elastic behavior of the Earth's crust, hence they can be recovered through the analysis of the Earth's seismic response and its temporal evolution. This is particularly true when earthquake codas, microtremors or microseismic noise are considered, as these are very sensitive to the effects of the often small perturbations in the Earth's elastic properties as they sample it both randomly and repeatedly [START_REF] Aki | Space and time spectra of stationary stochastic waves, with special reference to microtremors[END_REF][START_REF] Sato | Seismic wave propagation and scattering in the heterogeneous Earth[END_REF]. Much effort has been devoted to the study of waveform variations in space and time for the purpose of understanding the dynamic behaviour of the crust. Of particular interest are tectonically and volcanically active regions in which stress changes are frequent and may precede earthquakes and volcanic eruptions. Initially, almost all studies focused on the spatio-temporal behavior of coda waves, where the observation of variations in their amplitude found a possible application in the forecasting of volcanic activity [START_REF] Aki | Seismic monitoring and modeling of an active volcano for prediction[END_REF]. The inclusion of phase information to the analysis [START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doublets: an application to the Calaveras faults, California[END_REF] gave rise to a new approach which led to the detection of relative variations in seismic velocity between earthquake doublets and multiplets. In the same way, the seismic coda wave interferometry technique developed by [START_REF] Snieder | Coda wave interfermotry for estimating nonlinear behavior in seismic velocity[END_REF] has confirmed the existence of detectable precursory crustal changes [START_REF] Grêt | Monitoring rapid temporal change in a volcano with coda wave interferometry[END_REF][START_REF] Wegler | Increase of shear wave velocity before the 1998 eruption of Merapi volcano (Indonesia)[END_REF]), but is only practicable in cases where records of highly similar earthquake doublets are available.

More recently, seismic noise has become an increasingly popular and promising area of study, as it circumvents the limitations imposed by a dependence on the occurrence of seismic events. This is due to the possibility of retrieving seismic Green's functions from the cross-correlation of records of a random seismic wave field taken at various locations within a region of interest ⋆ (Weaver & Lobkis 2001;Lobkis & Weaver 2001;[START_REF] Campillo | Long-Range Correlations in the Diffuse Seismic Coda[END_REF][START_REF] Shapiro | Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise[END_REF][START_REF] Shapiro | High-Resolution Surface-Wave Tomography from Ambient Seismic Noise[END_REF][START_REF] Sabra | Extracting time-domain Green's function estimates from ambient seismic noise[END_REF]. Indeed, the use of ambient noise cross-correlations for monitoring has been shown to be robust even when conditions prevent the full reconstruction of the seismic Green's function (C. [START_REF] Hadziioannou | Stability of monitoring weak changes in multiply scattering media with ambient noise correlation: laboratory experiments[END_REF][START_REF] Weaver | On the correlation of non-isotropically distributed ballistic scalar diffuse waves[END_REF].

Because seismic noise is recorded continuously and does not depend on earthquake sources, these cross-correlation (cc) functions can be considered analogously to records from continuously repeating doublet sources placed at each station, and can be similarly used to extract observations of variations in seismic velocities.

The main idea for monitoring the evolution of seismic velocities over time using seismic noise is to compare "current" cross-correlation functions that represent the situation at a given time period to "reference" functions that represent an average background state of the studied media. We can distinguish between two different approaches that are used for the extraction of seismic velocity variations from cross-correlations and operate in the time and frequency domains, respectively.

The first method, known as Coda Wave Interferometry, was described by [START_REF] Snieder | The theory of coda wave interferometry[END_REF], and later evolved to Passive Image Interferometry [START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at Merapi volcano, (Indonesia)[END_REF]Wegler et al. 2009) for noise sequence cross-correlations. The second method has been named Moving-Window Cross-Spectral Analysis (MWCS) by [START_REF] Ratdomopurbo | Monitoring a temporal change of seismic velocity in a volcano: application to the 1992 eruption of Mt. Merapi (Indonesia)[END_REF] and is the focus of the present study. In fact, although approaches in both the time and frequency domains have found interesting applications showing similar sensitivities (Wegler et al. 2009), the MWCS technique has the advantage of operating in the frequency domain, where the bandwidth of coherent signal in the correlation function can be clearly defined.

The main goal of this paper is to assess the accuracy of the velocity variations measured from noise cross-correlations with the MWCS technique and, in particular, how this accuracy depends on the quality (i.e., signal to noise ratio, SNR) of the reconstructed cc functions. We start by briefly introducing the main concepts of the MWCS method with most of the technical details described in Appendix A. Then, we use a set of noise cross-correlations computed from records of the seismic stations of the Piton de la Fournaise volcano (La Réunion) monitoring network shown in Figure 1.

First, we study the convergence of these cc functions and their fluctuations in the time and frequency domains. Then, we construct a set of synthetic reference and current cc functions by stretching the observed cc functions to mimic velocity variations within the media and by adding random noise with spectral properties similar to observed random variations. We apply the MWCS measure to these synthetic cc functions and compare the inferred velocity variations with known a-priori introduced values. We finally propose threshold values of SNR above which small velocity variations can be reliably retrieved and subsequently interpreted.

MOVING-WINDOW CROSS-SPECTRAL ANALYSIS

The Moving-Window Cross-Spectral technique was first introduced by [START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doublets: an application to the Calaveras faults, California[END_REF] for the retrieval of relative velocity variations between earthquake doublets. More recently, Brenguier et al. (2008a,b) exploited this technique by applying it to seismic noise records, taking advantage of the possibility of treating noise cross-correlations in analogy with doublets. Here we describe only the general purpose of the technique, leaving all computational details to Appendices A and C. This analysis is applied to time series which are computed by cross-correlating the noise sequences recorded at two different seismic stations, for all possible station pairs. The preliminary step for the analysis is to build up at least one reference and several current cross-correlations.

Since, for computational purposes, the continuous noise records are cut into short sequences (for example, one for each day or hour), it is necessary to stack a certain number of single cc's. In this framework, the reference and current functions are defined by the number of summed cc's: N ref and N cur , respectively. The only requirement is that N ref >> N cur to ensure the reference cc is representative of a background value, while the current cc contains information on the actual state of the crust.

For any couple of reference, cc ref , and current, cc cur , functions, the technique combines two steps. The first step consists in the computation of the time-delay between the two signals within a series of overlapping windows. The second step is the evaluation of the relative velocity variation associated to the current function with respect to the reference. In this second step, it is assumed, for simplicity, that the seismic wave propagation velocity is perturbed homogeneously within the studied media.

It is important to note that the first operation is executed in the spectral domain, through the study of the phase of the cross-spectrum, allowing for precise selection of the frequency band on the basis of the coherency between the two windowed cc's (see Appendix A1). Each computed delay-time corresponds to a cross-correlation lag-time, which is taken as the central point of the window. Therefore, the second step involves the evaluation of the trend, δt/t, of the delay-time estimates over the whole length of the signals (see Figure A1). The slope of their linear regression indicates, to a first approximation, the relative homogeneous velocity perturbation of the current cc with respect to the reference cc.

Critical points in the MWCS analysis are the choices of N ref , N cur , the length and overlap of each window and the total number of windows which are used. These are all required for the first step. The choice of these parameters will depend on the characteristics of the cc functions such as their length, frequency content and how fast the signal decays below the noise level. The aim of our work is to test the reliance of both the resolution and accuracy of the measurements on the quality of the cc functions, which can be quantified in terms of their signal to noise ratio.

RANDOM FLUCTUATIONS AND CONVERGENCE OF OBSERVED NOISE

CROSS-CORRELATIONS

Measuring the signal to noise ratio (SNR) of a stacked set of cross-correlation functions is needed to distinguish between stacks from which reliable delay-times can be measured, and those from which they cannot. Furthermore, the simulated SNR of the cross-correlations we use in our tests must be compatible with this measure. We employ the method described by [START_REF] Larose | Reconstruction of Rayleigh-Lamb dispersion spectrum based on noise obtained from an air-jet forcing[END_REF], which is summarized below.

First, to estimate the level of noise, σ(N, t), in a stack, we measure the variation between each constituent cross-correlation function, cc(t), at each lag-time, t, as follows

σ(N, t) = cc(t) 2 -cc(t) 2 N -1 (1)
Here, • denotes the average over N single functions. We then measure the level of signal, s(N, t), in the stacked cross-correlation by taking its Hilbert envelope

s(N, t) = | cc(t) + iH( cc(t) )| (2)
where H(•) denotes the Hilbert transform of the stacked function cc(t) and i is the imaginary unit. After we smooth s(N, t) and σ(N, t) with a ten-second-wide sliding cosine window, the SNR of the stacked cross-correlation function can be estimated

SNR(N, t) = s(N, t) σ(N, t) (3)
Figure 2 demonstrates the measurement of SNR using this method. The plotted cross-correlations are from stations DSR and TCR near the summit of Piton de la Fournaise volcano, La Réunion (Figure 1), during the thirty days preceding an eruption on its northern flank [START_REF] Peltier | Suivi, modélisation et évolution des processus d'injections magmatiques au Piton de La Fournaise (Réunion), à partir d'une analyse croisée des données de déformation[END_REF]. Neither the variation between the daily functions nor the estimated signal are constant with t. The resulting signal to noise ratio, however, is less variable. For the purpose of our tests, we simulate a signal to noise ratio which is constant for all t when we add noise to our cross-correlation functions. Figure 4 shows the relationship between the error estimated by the MWCS technique (e b defined by Equation A.12 in Appendix A2) and the SNR calculated above. Here, current-and reference functions are formed by grouping daily cross-correlations from stations DSR-TCR into 30-day and 300-day stacks, respectively. For each current and reference function pair between the years 1999 and 2006, we measure delay-times within 6 second-wide lag-time windows. We then attribute the mean current function SNR within each window to its corresponding delay-time δt.

Finally, we calculate δt/t along with its accompanying error as explained in Appendix A. This error, expressed relative to δt/t, is plotted against the mean of the attributed SNR values. For each of the plotted ranges of |δt/t|, these errors appear to be anti-correlated with SNR.

SYNTHETIC CROSS-CORRELATION FUNCTIONS

To test the sensitivity of the MWCS technique, we construct a synthetic data set as follows: first, we take a reference function from real cross-correlations of seismic noise, then stretch it to simulate a series of homogeneous seismic velocity changes. This stretch is achieved by resampling the crosscorrelations with a Fourier-transform based interpolation. Effectively, this involves zero-padding the cross-correlation in the Fourier-domain, then taking the inverse transform. When the original sampling interval is applied, the interpolated cross-correlation becomes a stretched version of the original. We then add random noise to each stretched cross-correlation function to simulate a set of possible signal to noise ratios. Finally, treating the original function as the reference, and each stretched, noise-added function as the current function, we attempt to recover the applied stretch using the MWCS method and to see how the resulting errors depend on the level of the added noise. For these tests it is important to use synthetic noise with properties close to the real random fluctuations of the observed cc functions. Therefore, we first characterize the spectra of these observed fluctuations and then propose a procedure to simulate a random noise series with defined spectral properties.

Spectrum of the observed random fluctuations

A simple way to view a pair of stacked current and reference cross-correlation functions is to treat the current function, cc cur (t), as a contaminated version of the relatively noiseless reference function, cc ref (t).

cc cur (t) = cc ref (t) + n(t) (4)
To observe one realization of the impinging noise, n(t), we simply subtract the current function from the reference function. The difference between the reference and current stack is taken as an observation of noise. Small lag-times are omitted from the calculation as they are not used for delay-time estimation. c and d: Amplitude spectra of the negative-lag and positive-lag segments, respectively, of the observed noise (solid curves) and the reference stack (dashed curves).

to be relatively low compared to that of the reference function, examination of its spectrum reveals that the amplitude of the noise is at least comparable to that of the reference function at certain frequencies.

For every station pair, we calculate n(t) for all available current and reference function pairs between the years 1999 and 2006. Then, by averaging their squared-amplitude spectra, we ob-serve a typical noise spectrum which we can use to contaminate our synthetically stretched crosscorrelation functions (e.g., Figure 6c and6d).

Simulating random noise with pre-defined spectra

To contaminate our synthetic cross-correlation functions, we randomly perturb each value by an amount drawn from a Gaussian distribution. The standard deviation of this distribution is chosen as follows

σ synth (t) = s(t) SNR(t) (5)
where s(t) is the signal envelope of the synthetic function, and SNR(t) is the desired signal to noise ratio. We discuss whether or not the use of a Gaussian distribution to simulate noise in this way is appropriate for our cross-correlations in Appendix D.

To ensure that the noise exists in the appropriate frequency band, we apply the method described by [START_REF] Percival | Simulating Gaussian random processes with specified spectra[END_REF] to produce a random time-series which shares the same spectrum as that observed in real data. After normalizing to unit standard deviation, we scale the resulting noise by σ synth (t), then add it to our synthetic cross-correlation function.

Figure 7 demonstrates how noise is prepared in this way. The original function is a stack of all available daily cross-correlations computed between stations DSR and TCR in the year 2002. Noise is generated using the spectrum which is observed for this station pair, then scaled to produce a constant signal to noise ratio. This noise targets the frequency range in which real noise is observed and in which delay-times are to be later measured.

cross-correlation

RESULTS OF SYNTHETIC TESTS

For every station pair, we simulate 1000 random synthetic realizations of current cc functions with predefined stretching coefficients mimicking velocity perturbations and predefined signal to noise ratios. Then, we analyze the resulting set of δt/t (stretch) estimates and their accompanying errors. We show results of sensitivity tests for vertical-component records from two pairs of stations on Piton de la Fournaise volcano, La Réunion. For each station pair (BOR-SFR, and DSR-TCR),

we stack every available daily cross-correlation function between the years 1999 and 2006 to construct our reference function, then filter between 0.1 and 1.0 Hz.

To construct current functions, we stretch these reference cross-correlations by a range of values, S = {0.01%, 0.02%, ... , 0.10%}, then add noise to simulate signal to noise ratios of SNR = {1, 2, ... , 10}. Finally, we use six-second-wide lag-time windows which overlap by three seconds to compute delay-times (See Appendix A).

For each pair of simulated stretch and SNR values (S, SNR), we obtain 1000 stretch estimates,

S i , i ∈ [1, 1000],
and their associated least-squares errors, e i (standard deviations) from the MWCS technique. Figure 8 shows the distribution of these estimates for two different values of SNR. In both cases, these estimates form an approximately bell-shaped distribution centered around the true stretch of 0.05%. In the case of low SNR, these estimates form a wide distribution (Figure 8a) due to the high level of noise in the cross-correlations, and we cannot confidently recover δt/t.

When SNR is increased (Figure 8b), the distribution narrows, and δt/t is better resolved.

To quantitatively assess the level of systematic error in each set of estimates, we calculate their relative bias as follows Figure 9 shows this measure for stations DSR-TCR (Figure 9a) and BOR-SFR (Figure 9b). For both station pairs, the relative bias is never more than a few percent, provided the simulated stretch and SNR are large enough (say, above 0.02% and 3, respectively). This suggests the MWCS method introduces very little systematic error.

To assess the total relative error over the 1000 trials for each (S, SNR) pair, we calculate their misfit from the true stretch as follows

e total (S, SNR) = 1 S × (S i -S) 2 1000 -1 (7)
This incorporates both the systematic and the random error in each set of 1000 S i estimates. Figure 10 shows these measures for the two station pairs described above. Here, colors indicate the level of error, expressed as a percentage of the true stretch. As expected, this error decreases as either the applied stretch or the simulated SNR are increased.

We compare this error, e total , evaluated from the synthetic test with errors evaluated from the least-squares fit during the MWCS analysis (Appendix A2). For every synthetic current cc function, we evaluate the least-squares error and then compute its mean value for a given pair of stretching coefficient and SNR e i (S, SNR). Figure 11 shows the ratio of errors estimated from the synthetic test and from the MWCS least-squares fit for station pairs DSR-TCR and BOR-SFR.

For both station pairs, the least-squares error underestimates the total variability of the targeted velocity variations by a factor of around six for most values of SNR and for all applied stretching coefficients (Figure 11). We address the cause of this discrepancy in Appendix B.

Finally, we plot e total against SNR (Figure 12) to see if the estimates we obtain during our tests exhibit the inverse relationship between cross-correlation quality and δt/t error that we see for real data (Figures 4 and5). When viewed on a log-log scale, these results exhibit a clear anticorrelation.

Reassuringly, the similarity between this plot and Figures 4 and5 demonstrates the consistency between the SNR we simulate and the SNR we measure from real data.

DISCUSSION AND CONCLUSIONS

To assess the accuracy of the velocity variations measured from noise cross-correlations (CC) with the MWCS technique we constructed a set of synthetic CC functions corresponding to known media velocity variations (stretching coefficients) and perturbed by random noise with statistical properties similar to those observed at the stations of the Piton de la Fournaise seismic network.

Our analysis resulted in simple relations between the accuracy of the recovered velocity variations and the Signal-to-Noise Ratio (SNR) of the analyzed CC functions (Figure 10). In turn the SNR is on average simply related to the duration of the noise record from which the CC function was computed (Figure 3). These results provide us with a simple guidance on how to choose an optimal stack duration to recover a desired level of media velocity variations. In particular, for the case of the seismic stations on Piton de Fournaise volcano, our analysis indicated that recovering a relative velocity perturbation of 0.1% from a single pair of stations requires an SNR of ∼ 5 that can be obtained by stacking a few tens of days of noise correlations. This implies in particular that the accuracy of measurements presented by Brenguier et al. (2008b) by averaging measurements from many station pairs could be barely achieved from analyzing a single pair of stations.

Another important result of our analysis performed with synthetic CC functions is that the formal error computed from the linear regression within the MWCS technique does not match the true uncertainty of the recovered relative delay-times. In the case of our tests, the true error appears to be around six times greater than that which is estimated. This mismatch is mainly due to the fact that the least-squares error is not uniquely defined, but depends on the parameters used in the the application of the MWCS technique (See Appendix B). For a particular choice of parameters, this error may underestimate the real uncertainty of the recovered velocity variations. A further explanation is that the MWCS technique effectively uses only one realization of the CC function with a relatively short duration (because of the fast decay of the coda part of the recovered Green's functions). This single and short realization is not representative of the full variability of the CC functions. The sampling can be improved by using multiple pairs of stations simultaneously as has been done by Brenguier et al. (2008a,b). Nonetheless, the factor relating the MWCS error with the total error is roughly independent of both SNR and the media velocity variation. Furthermore, we observe the same factor (Figure 11) for both pairs of stations considered in our study, BOR-SFR and DSR-TCR. This means that, in the case of Piton de la Fournaise seismic noise crosscorrelations, and for this particular choice of parameters, we can apply a correction to the MWCS errors by simply multiplying their values by a factor of ∼ 6.

A main conclusion from our study is that before systematically applying noise-based MWCS monitoring of temporal media changes in a particular setting, it is important to investigate the statistical properties of the seismic noise and the convergence of noise correlations. This analysis is necessary to establish the correction factor for the MWCS errors and also the optimal durations of correlated time series. So far, our results indicate that recovering relatively weak velocity changes associated with moderate volcanic activity (Brenguier et al. 2008b), intermediate-size earthquakes (Brenguier et al. 2008a) or with seasonal variations [START_REF] Meier | Detecting seasonal variations in seismic velocities within Los Angeles basin from correlations of ambient seismic noise[END_REF]) requires stacking correlations from a few tens of days and averaging measurements from many pairs of stations. Further improvement of temporal and spatial resolution of the MWCS measurements could be eventually achieved by applying additional steps in the data processing such as data adaptive filtering [START_REF] Baig | Denoising seismic noise cross correlations[END_REF].

C(ν) = X(ν) |F ref (ν)| 2 • |F cur (ν)| 2 (A.3)
Here, the overlines indicate smoothing, which in our case is obtained by applying a sliding raised-cosine function with a half-width of 0.1 Hz to the energy density spectra of the two Fouriertransformed time-series and to the real and imaginary parts of the complex-valued cross-spectrum.

The cross-coherence ranges between zero and one, with maximum values approached at those frequencies where the two spectral densities are highly similar.

The time-delay between the two cross-correlations can be found in the (unwrapped) phase, φ(ν), of the cross-spectrum, which will be linearly proportional to frequency.

φ j = m • ν j , m = 2πδt (A.4)
The time shift, δt i (subscript i for the i th window), between the two signals is estimated from the slope m of a linear regression of the samples, j = l, ..., h, within the frequency range of interest (see Figure A2, panel c). During the regression, a weight w j , which depends on the crosscoherence at each sampled frequency, is assigned to each cross-phase value.

w j = C 2 j 1 -C 2 j • |X j | (A.5)
Unlike [START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doublets: an application to the Calaveras faults, California[END_REF], these weights incorporate both the cross-spectral amplitude and the cross-coherence. This generates more differentiated weights in cases where the crosscoherence is relatively constant but the cross-spectral energy is variable. Figure A2, panel b shows such an example. This choice of weighting is described in more detail in Appendix C.

Using a weighted least-squares inversion, the slope m is estimated as

m = h j=l w j ν j φ j h j=l w j ν 2 j (A.6)
The associated error, e m , is calculated using the rule of propagation of errors

e m = j w j ν j i w i ν 2 i 2 σ 2 φ (A.7)
where σ 2 φ is the squared misfit of the data to the modeled slope and is calculated as

Figure A1. a:

A pair of a real reference (from Piton de la Fournaise) and a synthetic current crosscorrelation, along the time interval between -30 and 30 s. The current function is generated by stretching the reference by 0.1%. The black rectangle delimits one of the windows used in the analysis (6s long, sliding 3s). b: The second step in the MWCS analysis: a linear regression through weighted least-squares over the time delays that have been computed during the first step. Each time-delay is associated with the mean time in its sliding window. The straight red line is the fitted slope and the dotted lines highlight the error margins.

σ 2 φ = j (φ j -mν j ) 2 N -1 (A.8)
Following Equation A.4, the time delay, δt, and its error, e δt , between the two signals are taken by simply dividing m and e m , respectively, by 2π.

Repeating this process for all windows, we obtain N w delay-time estimates between the two cross-correlation functions, each corresponding to the central time, t i (i = 1, ..., N w ), of the window in which it was measured.

It is important to keep in mind that, for a given frequency range, e δt is inversely proportional to the square-root of the number of values that are used in the inversion. This means that if the windowed cross-correlations are zero-padded prior to Fourier transformation, the error estimate will be artificially reduced. Multiplying e δt by N f f t , where N f f t is the number of points in the Fourier-transformed time series, removes this dependence. Coherency calculated for the two windowed signals at all frequencies. c: Linear regression for the phase displacement along the frequencies of interest (0.2 to 0.9 Hz), which have been marked with red asterisks and errorbars (1/w 2 j ).

An important feature of this formulation is that, for a given correlation-time interval, the error of the relative velocity variation, e b , is inversely proportional to the square-root of the number of delay-times that are used in the regression. Consequently, if the number of sliding windows N w is increased by reducing the time-step between consecutive windows, then the error will be artificially reduced. This is similar to the dependence of each delay time error e δt (Equation A.7) on the number of points used in the Fourier transform of the windowed data. Multiplying the estimated error by N w is one way to remove this dependence.

Figure A3

. The dependence of δt/t errors on N f f t and N w . Total errors (a, Equation 7) and average leastsquares errors (b, Equation A.12) are plotted as a function of the time-step between six-second-wide sliding windows. One curve is plotted for each value of N f f t (see legend). Cross-correlation functions are taken from station pair BOR-SFR. In all cases, the MWCS technique is applied 1000 times to cross-correlations that have been stretched by 0.05% and have a simulated signal to noise ratio of 5.

APPENDIX D: DISTRIBUTION OF OBSERVED FLUCTUATIONS

In Section 4.2 we simulate noisy cross-correlation functions by contaminating them with a random time-series whose squared amplitude spectrum mimics that of the fluctuations we observe in real data. This time-series is drawn from a Gaussian distribution with random phase. In this section, we determine whether such a series is representative of the fluctuations that exist in real crosscorrelations.

As described in Section 4.1 we observe the real fluctuations in our cross-correlations by taking the difference between corresponding current and reference functions. The cross-correlations we use in the following examples are from station pair BOR-SFR on Piton de la Fournaise volcano, and were measured during the period between 1999 and 2006.

and those we obtain in our simulated tests (Figure 12) suggests that the methods we use to create synthetic noise and evaluate the level of fluctuation in real cross-correlations are adequate for the purposes of these tests.

  Map of stations (named black squares) used in this study. inset: La Réunion island. The black box outlines the part of Piton de la Fournaise shown in the main panel.

Figure 3

 3 Figure 3 shows how the SNR of a cross-correlation stack depends on N. Here, we stack various numbers of consecutive daily cross-correlations from stations DSR and TCR. Days of missing data and of eruptive activity are skipped, and the plotted values of N are the number of remaining days. This plot shows that SNR grows at a rate which is just less than proportional to √ N. While SNR appears to increase monotonically with N, it may in many instances be affected by drastic changes in the geology of the region under investigation. In this example, we avoid the collapse of Dolomieu crater at the start of 2007 on the summit of Piton de la Fournaise volcano by only using cross-correlations between the years 1999 and 2006.

Figure 2 .

 2 Figure 2. The calculation of SNR. a: A set of thirty single-day cross-correlation functions (grey curves)and their stacked mean (solid black curve). The dotted black curve is the signal envelope of the stack, and is smoothed with a ten-second-wide cosine window. b: The smoothed noise measured from this set of cross-correlations. c: The resulting SNR is the ratio of the signal envelope and the noise.

Figure 3 .

 3 Figure 3. A plot of SNR versus number of stacked days (N). Results are separated into bins of ln(N ) and ln(SN R). Counts are plotted as shades of grey after normalization within each column of the plotted grid. Dark shaded bins have high counts relative to lightly shaded bins. SNR is averaged over |t| > 10 seconds. The dashed curve and displayed slope are from a linear regression of the plotted values.

Figure 4 .

 4 Figure 4. Relative δt/t error estimates versus SNR for station pair DSR-TCR. Each panel pertains to a different range of |δt/t|. Dashed curves and displayed slopes are obtained via linear regressions of the plotted values.

Figure 5 .

 5 Figure 5. Lines of best fit determined from relative δt/t error estimates for (a) stations DSR-TCR and (b) BOR-SFR. One line is plotted for each value of |δt/t| (see legend).

Figure 6 Figure 6 .

 66 Figure6shows an example of a current and a reference function computed for stations DSR and TCR on Piton de la Fournaise volcano, La Réunion. The reference function is a stack of the 300 days preceding an eruption on the volcano's south-eastern flank[START_REF] Peltier | Suivi, modélisation et évolution des processus d'injections magmatiques au Piton de La Fournaise (Réunion), à partir d'une analyse croisée des données de déformation[END_REF], while the current function is a stack of the last thirty of those days. Although the level of the noise appears

Figure 7 .

 7 Figure 7. Addition of random noise to a synthetic cross-correlation function from stations DSR and TCR.The signal level (a, dotted curve) is estimated from the original function (solid curve). Panel b shows the spectrum of the cross-correlation function. Random Gaussian noise (c) is produced with a spectrum (d) that simulates that of the noise which is observed for the station pair. This noise is scaled in time (e and f), then added to the original cross-correlation function (g and h).

Figure 8 .

 8 Figure 8. Histograms for two sets of stretch estimates. a: The applied stretch is 0.05 % and the simulated SNR is 3. b: The applied stretch is unchanged, but the simulated SNR is increased to 8. The reference cross-correlation is from stations DSR and TCR.

Figure 9 .

 9 Figure 9. Relative bias calculated for each (S, SN R) pair using Equation 6. Cross-correlation functions are taken from (a) stations DSR-TCR and (b) stations BOR-SFR. Light colours indicate a small relative bias.

Figure 10 .

 10 Figure 10. Total errors calculated for each pair using equation 7. Cross-correlation functions are taken from (a) stations DSR-TCR and (b) stations BOR-SFR. Cold colours indicate (S, SN R) values for which δt/t is well resolved using the MWCS technique.

Figure 11 .

 11 Figure 11. Ratio between the total error and the error estimated from the MWCS least-squares fit as a function of SN R. Cross-correlation functions are taken from (a) stations DSR-TCR and (b) stations BOR-SFR. One line is plotted for each simulated stretch (see legend).

Figure 12 .

 12 Figure 12. Total error plotted as a function of SNR. Cross-correlation functions are taken from (a) stations DSR-TCR and (b) stations BOR-SFR. One line is plotted for each simulated stretch (see legend in Figure 11) along with its estimated slope (dashed lines).

Figure A2 .

 A2 Figure A2.Example of how a delay-time is measured between two windowed cross-correlations. a: zoomin of the two cross-correlations within a single lag-time window (black rectangle in FigureA1, panel a). b: Coherency calculated for the two windowed signals at all frequencies. c: Linear regression for the phase displacement along the frequencies of interest (0.2 to 0.9 Hz), which have been marked with red asterisks and errorbars (1/w 2 j ).

  Figure A4. a) Relative error of time delay estimations versus SNR for different weights. Each symbol corresponds to one of the three definitions of w j . b) Final results of the MWCS analysis in varying the weights and the SNR level. The symbols match those in panel a), the black horizontal line shows the real value of stretching between the two CCs.

ACKNOWLEDGMENTS

The seismogram data used in this study were provided by the Observatoire Volcanologique du Piton de la Fournaise. We thank Michel Campillo, Berenice Froment and Celine Hadziioannou for helpful discussions and G. Moguilny for maintaining the Cohersis cluster on which all computations were performed . This work was supported by Agence Nationale de la Recherche (France) under contracts ANR-06-CEXC-005 (COHERSIS) and ANR-08-RISK-011 (UNDERVOLC) and by a FP7 ERC Advanced grant 227507 (WHISPER).

 Wegler, U., Nakahara, H., Sens-Schönfelder, C., Korn, M., & Shiomi, K., 2009. Sudden drop of seismic velocity after the 2004 M w 6.6 mid-Niigata earthquake, Japan, observed with Passive Image Interferometry, J Geophys Res, 114(B06305), 1-11.

APPENDIX A: MWCS

In the following section, the method of Moving-Window Cross-Spectral analysis (MWCS) is described in the context of stacked reference and current cross-correlation functions.

A1 Time-delay computation

The first step in the MWCS analysis is the calculation of delay-times, δt, between the two crosscorrelation functions within a series of overlapping lag-time windows.

Each cross-correlation function is divided into N w windows, one for each delay-time measurement. The choice of window length, overlap, and N w will generally depend on the frequency content and the SNR of the cross-correlation functions under consideration. The windowed segments are mean-adjusted and cosine-tapered before being Fourier-transformed into the spectral domain.

In Figure A1 (a), an example of a windowed pair of cross-correlation functions is shown. The cross-spectrum, X(ν), between the two windowed time-series is calculated as follows

where F ref (ν) and F cur (ν) are Fourier-transformed representations of the windowed time-series, ν is frequency in Hz and the asterisk denotes complex conjugation. For our purposes, it is more useful to represent the complex cross-spectrum by its amplitude |X(ν)| and phase φ(ν)

One requirement of cross-spectral time-delay estimation is that, aside from being shifted in time, the two windowed time-series are similar. Such similarity is quantitatively assessed using the cross-coherence C(ν) between their energy densities:

A2 Velocity variation results

To a first order approximation, we can consider a stress field perturbation which acts homogeneously over the region sampled by the cross-correlated seismic noise. Under this assumption, the resulting seismic velocity perturbation δv/v within that region will also be homogeneous, and be manifest as a stretching -δt/t of the current cross-correlation function relative to the reference function. This stretching is constant over t, and is numerically the opposite of the velocity perturbation [START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doublets: an application to the Calaveras faults, California[END_REF].

Consequently, to recover δv/v, we apply a linear regression to the N w delay-time measurements (Figure A1).

where the coefficient a represents a possible instrumental drift [START_REF] Stehly | Traveltime measurements from noise correlation: stability and detection of instrumental time-shifts[END_REF], and b corresponds to the relative time variation δt/t. Again, we can estimate these two parameters through a weighted least-squares inversion. Here, the weights, p i , are determined using the estimated error of each time-delay measurement:

with variance

with variance

where t = p i t i / p i , δt = p i δt i / p i and δt 2 = p i δt 2 i / p i are weighted means of t, δt and t 2 , respectively.

APPENDIX B: DEPENDENCE OF ERRORS ON MWCS PARAMETERS

In Section 5, we observe a discrepancy (Figure 11) between the total errors we obtain from the distribution of each set of 1000 stretch estimates (e total , Equation 7) and the estimated least-squares error defined by Equation A.12. One explanation for this is the dependence of the estimated error on the number of sliding windows, N w , into which our cross-correlations are divided (see Appendix A2). In turn, this value is affected by the delay time errors, e δt (Equation A.7 in Appendix A1), which themselves are dependent on the number of points, N f f t , used to transform the windowed cross-correlations into the Fourier domain.

We observe the behaviour of the total error and the estimated error as these two parameters are varied. To this end, we alter N w by adjusting the time-step between consecutive six-second-wide windows, and N f f t by zero-padding the windowed cross-correlations prior to Fourier transformation. Figure A3 shows the total error (Figure A3a) and the average least-squares error (expressed relative to δt/t, Figure A3b) we obtain when a stretch of 0.05% and a signal to noise ratio of 5 are simulated for station pair BOR-SFR. Each point corresponds to 1000 trials for a given choice of N w and N f f t . These plots demonstrate the inverse proportionality between the estimated error and the square-root of both N w and N f f t . Interestingly, the total error also appears to increase slightly with the time-step, suggesting that a choice of broadly overlapping windows improves the precision of the relative travel time measurements that are obtained. However, the associated error estimates must be calibrated in order to accurately evaluate the true precision of these measurements.

APPENDIX C: TEST ON THE WEIGHTS

In Section A1 we introduce weights w j (Equation A.5) to be associated to each φ j when estimating time shifts between cross-correlation functions. In this section, we test the influence of these weights on the results. In order to search for the most suitable formulation of w j , we compare the accuracy of the yielded estimates for three different weight definitions:

Using these weights, we apply the MWCS analysis to a reference and a synthetic current function which has been perturbed from the reference by stretching it to 0.1 %. Starting from this noiseless current function, we add noise (as described in Section 4.2) to reach final signal to noise ratios of 10, 5, 2 and 1. The resulting estimates are shown in Figure A4 (a and b, respectively) for the relative error on time delay computations, and for the relative velocity variation recovered (named stretch). These measurements are in close agreement with one another, revealing only a slight dependence on the weights that are used. We choose to use w j (3) as it produces differentiated weights in cases of near constant coherence, and performs slightly better than the other schemes in these tests. Furthermore, these findings stress the importance of the noise level on the resolution of the MWCS technique as the errors shown in Figure A4 Our tests show that the fluctuations we observe in real cross-correlations deviate slightly from a simple Gaussian distribution with random phase, mostly during eruptions. Therefore, the analysis we present in this paper relies on the fluctuations being Gaussian. One way to improve this analysis for co-eruptive periods would be to characterize the true noise distribution and randomly draw from it when simulating noisy cross-correlations. Nonetheless, the similarity between the measurement errors we observe when applying the MWCS technique to real data (Figures 4 and 5)


G Fau

N Gascoin

P Gillard

M Bouchez

J Steelant

Fuel Pyrolysis through Porous Media: Coke Formation and Coupled effect on Permeability

The development of hypersonic vehicles (up to Mach 10) leads to an important heating of the whole structure. The fuel is thus used as a coolant. It presents an endothermic decomposition with possible coke formation. Its additional permeation through the porous structure involves internal convection. This implies very complex phenomena (heat and mass transfers with chemistry). In this paper, the n-dodecane pyrolysis is studied through stainless steel porous medium up to 820 K and 35 bar (supercritical state). The longitudinal profiles of chemical compositions inside the porous medium are given thanks to a specific sampling technique with off-line Gas Chromatograph and Mass Spectrometer analysis. By comparison with previous experiments under plug flow reactor, the conversion of dodecane is higher for the present experimental configuration. The pyrolysis produces preferentially light gaseous species, which results in a higher gasification rate for a similar pyrolysis rate. The effects of the residence time and of the contact surface area are demonstrated. The transient changes of Darcy's permeability are related to the coke formation thanks to previous experimental relationship with methane production. A time shift is observed between coke chemistry and

A c c e p t e d M a n u s c r i p t 2 permeability change. This work is quite unique to the author's knowledge because of the complex chemistry of heavy hydrocarbon fuels pyrolysis, particularly in porous medium.

Keywords: Permeation, pyrolysis, porous medium, coke, active cooling, catalytic effect. of the accumulated combustion heat. To challenge this problem, several cooling systems have been developed. Two major axis have been particularly studied: the material withstanding against high pressures and high temperatures on the one hand and a reliable cooling of the combustion chamber on the other hand [START_REF] Steelant | ATLLAS: Aero-Thermal Loaded Material Investigations for High-Speed Vehicles[END_REF][START_REF] Bouchez | Material-Aero-Thermal Interaction Computations in the ATLLAS European Programme[END_REF][START_REF] Steelant | Achievements obtained on Aero-Thermal Loaded Materials for High-Speed Atmospheric Vehicles within ATLLAS[END_REF]. The European ATLLAS project (Aerodynamic and Thermal Load interactions with Lightweight Advanced materials for high Speed)

supported by the European Space Agency is investigating this topic, notably the transpiration through Ceramic Matrix Composite (CMC) structures. This cooling technique presents difficulties, which needs to be investigated. The coolant can react with the material (catalytic effect) and the coke formation can modify the structure's permeability [START_REF] Gascoin | Characterisation of coking activity during supercritical hydrocarbon pyrolysis[END_REF]. As a result, the coupling of heat and mass transfer is directly impacted by the chemistry and such phenomenon need to be more understood. For this reason, the present work presents experiments in order to investigate the fuel pyrolysis through porous media and the associated phenomena (permeation, catalysis and coke production).

I . 1. Permeation through porous media

The permeation process is a complex phenomenon widely studied in various research areas such as nuclear power [START_REF] Hino | Helium gas permeability of SiC/SiC composite used for in-vessel components of nuclear fusion reactor[END_REF] or geology [START_REF] Wang | An improved permeability model of coal for coalbed methane recovery and CO2 geosequestration[END_REF]. This particular phenomenon is treated experimentally [START_REF] Tully | Fluid flow and heat transfer in a liquid rocket fuel injector[END_REF][START_REF] Langener | Experimental Investigation of Transpiration Cooling with Subsonic and Supersonic Flows at Moderate Temperature Levels[END_REF], numerically [START_REF] Murthy | Thermal dispersion effects on non-Darcy convection over a cone[END_REF][START_REF] Riccius | Coupled CFD analysis of hot gas and coolant flow in effusion-cooled combustion chambers[END_REF] or mathematically [START_REF] Kim | Fully discrete mixed finite element approximations for non-Darcy flows in porous media[END_REF]. A first approach to understand this specific phenomenon is to express the well-known Darcy's law:

L . g . P KS Q r D = (1)
This law is a linear part (only available for laminar flow) and should be completed by an additional one (for so-called turbulent flow). Thus, the next equation follows [START_REF] Soller | Experimental Investigation of Cooling Techniques and Materials for Highspeed Flight Propulsion Systems[END_REF]:

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

F D K ² V . K V . µ L P r + = D (2)
This law (known as the Brinkman's equation) can be derived and expressed with a large panel of relationships according to the application fields. Additional terms, like the Stokes term [START_REF] Ligaarden | On the Stokes -Brinkman Equations for Modeling Flow in Carbonate Reservoirs[END_REF], can be added and the pressure loss can be written in several ways [START_REF] Gascoin | Comparison of two permeation test benches and of two determination methods for Darcy's and Forchheimer's permeabilities[END_REF]. Depending on the selected expression, the fluid properties are given either by inlet or mean conditions.

The permeation mainly depends on the material structure which can be considered at different scales (macro, micro and also meso scale) [START_REF] Zhou | Analytic characterization of the permeability of dual-scale fibrous porous media[END_REF]. Due to the complexity of composite materials, it is chosen in this paper to consider metallic samples which are homogeneously made of stainless steel agglomerated particles.

I . 2. Pyrolysis of hydrocarbons and hypersonic conditions

The general scheme of pyrolysis is composed of the following successive and simultaneous reactions: random scission (the chain scission produces the free radicals), H-abstraction and beta scission (production of unsaturated compounds with one or two terminal double bond) [START_REF] Safarik | The thermal decomposition of hydrocarbons. Part 1. n-Alkanes (C 5)[END_REF]. These reactions appear in primary and secondary mechanisms. The degree of decomposition is highly dependent on the operating conditions (temperature, pressure, type of flow, nature of reactor) [START_REF] Billaud | Thermal coupling of methane in a tubular flow reactor: parametric study[END_REF][START_REF] Murphy | Analysis of products of hightemperature pyrolysis of various hydrocarbons[END_REF][START_REF] Chen | The Thermal decomposition of methane. I. Kinetics of the primary decomposition to C2H6 + H2; Rate constant for the homogeneous unimolecular dissociation of methane and its pressure dependence[END_REF][START_REF] Chen | The themal decomposition of methane. II. Secondary reactions, autocatalysis and carbon formation; non-Arrhenius behaviour in the reaction of CH3, with ethane[END_REF][START_REF] Zámostný | Experimental study of hydrocarbon structure effects on the composition of its pyrolysis products[END_REF]. Thermal cracking of hydrocarbons have been widely studied in petrochemical industry [START_REF] Murphy | Analysis of products of hightemperature pyrolysis of various hydrocarbons[END_REF][START_REF] Dunkleman | Pyrolysis of propane in tubular flow reactors constructed of different materials, in: Industrial and Laboratory Pyrolyses[END_REF][START_REF] Kassel | The thermal decomposition of methane[END_REF][START_REF] Chakraborty | High pressure pyrolysis of n-heptane[END_REF][START_REF] Zhou | Thermolytic reactions of dodecane[END_REF][START_REF] Billaud | N-Decane pyrolysis at high-temperature in a flow reactor[END_REF]. It appears than the bigger the molecule, the higher the number of reactions which occur. Considering dodecane pyrolysis, the number of reactions largely overpasses 1000. Light hydrocarbons like methane, acetylene or ethylene are produced such as heavy compounds (benzene) and even more complex ones (notably Polycyclic Aromatics Hydrocarbons -PAH -such as naphthalene and pyrene) which are precursors of solid carbon formation [START_REF] Herbinet | Thermal decomposition of n-dodecane: Experiments and kinetic modeling[END_REF]. This coke may interfer with the permeation process by lowering the porosity and reducing the contact surface between coolant and solid surface. Thus, the comprehension of the reaction pathway is one of the key to monitor and control such a A c c e p t e d M a n u s c r i p t 5 phenomenon. Last element to consider is that in real conditions, the fuel temperature can reach 1500 K and the pressure can be higher than 30 bar. This leads the fuel to be at supercritical state. To the authors' knowledge, no study does exist on fuel pyrolysis under permeation configuration. The challenge of the present work is thus to achieve supercritical fuel pyrolysis in porous media and to monitor the longitudinal chemical composition profile inside the structure thanks to appropriate sampling technique.

I . 3. Surface catalytic effect and coke formation

Some interactions between the molecules constituting the wall and the fluid may happen (surface reactions). This could modify strongly the reaction pathway and it can lead to a potential higher production of coke (which justifies considering specifically this point) [START_REF] Ghaly | Reactor Surface Effects During Propylene Pyrolysis, in: Industrial and Laboratory Pyrolyses[END_REF][START_REF] Gascoin | Etude et mesure de paramètres pertinents dans un écoulement réactif[END_REF][START_REF] Abraham | Etude et développement d'une méthode d'analyse par spectroscopie infrarouge appliquée à la pyrolyse d'hydrocarbures en conditions supercritiques et transitoires[END_REF]. Considering this element, the treatment of the inner wall could have a critical impact on the fuel degradation. It was shown that an oxygen treatment generally enhances the catalytic activity (considering a metallic reactor) and leads to a higher production of coke [START_REF] Ghaly | Reactor Surface Effects During Propylene Pyrolysis, in: Industrial and Laboratory Pyrolyses[END_REF][START_REF] Crynes | Pyrolysis of propane in tubular flow reactors. Kinetics and surface effects[END_REF]. The treatment by hydrogen sulphide permits to inhibit this phenomenon and to act as passivation.

A recent study also showed that coating the inner wall with zeolites or with a ceramic-like binder both limit the coke production, enhance the fuel decomposition and, consequently, increase the heat sink capacity of the material [START_REF] Li | Catalytic cracking of endothermic fuels in coated tube reactor[END_REF]. Hence, the choice of the wall material is a critical point and has to be considered with respect to the surface catalytic effect. In the present work, some investigations will be proposed in order to estimate its impact on the distribution of pyrolysis products.

Finally, the present study aims at investigating fuel pyrolysis in porous media to determine the intrinsic effect of the permeable material compared to existing data on fuel pyrolysis in open tubular reactor. The effects of the contact surface, of the residence time, of the catalytic effect and of heat and mass transfers are of interest. This work is even more important since it is the first frontier work, to the author's knowledge, between pyrolysis, catalysis and filtration. The final aim is to be able predicting the coke formation and its effect on permeability in order to ensure the cooling efficiency of regenerative cooling technique. Scientifically, major progresses are waited on the understanding of coupled heat and mass transfer with chemistry by experimentally observing the heterogeneous profiles of permeation and of pyrolysis inside the porous medium.

II.! Material and method

II . 1. Experimental bench

A high pressure and high temperature experimental bench (1200 K, 60 bar) is used with a novel sampling technique to investigate the pyrolysis products distribution as a function of operating conditions. The setup is derived from the one of the COMPARER project (COntrol and Measure of PArameters in a REacting stReam) conducted by the University of Orléans and MBDA-France [START_REF] Gascoin | Etude et mesure de paramètres pertinents dans un écoulement réactif[END_REF]. The liquid fuel is injected by a pump (80 bar, 0.5 g.s -1 ) inside a chemical reactor, which is placed in a furnace. This reactor is made of a preheated stainless steel tubular reactor (Ø int. = 3.45 mm, length = 550 mm) and of a dedicated cell (Ø int = 30 mm, length = 100 mm) in which the porous medium is inserted (filtration diameter of 16 mm only). A cooling reactor (Ø int. = 3.45 mm, length = 300 mm) is placed at the cell outlet to conduct the fluid out of the furnace before being analyzed transiently by Infra Red spectrometry (FTIR) for the gas phase and before being sampled. These samples under liquid and gas phases at ambient conditions are a posteriori analyzed separately by both Gas Chromatograph (GC) and Mass Spectrometer (MS). Further details on the bench description can be found in previous work [START_REF] Gascoin | High temperature and pressure reactive flows through porous media[END_REF][START_REF] Gascoin | Benchmark of Experimental Determination Methods of Gas[END_REF].

An improvement of the described bench consists on a novel technique of sampling (Figure 1). This allows obtaining the longitudinal distribution of pyrolysis mixture composition inside the permeable medium. First, sampling tubes have been placed as close to

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

the medium as possible. The porous sample is drilled to get samples in its middle. A good contact is ensured to get reliable samples. A complex system of small size cylinders, valves and pump has been designed to ensure the smallest possible residence time between the sampling point and the cylinder, a fast filling of these cylinders, a sufficient sample quantity for the GC/MS analysis, a good synchronization of the sampling during the test and finally the good representativeness of the samples compared to the real process itself. The dynamics of this system is of importance to have a "non intrusive" device. The following sampling locations are noted for this work: upstream (point A, Figure 1), inside (point B, Figure 1) and downstream (point C, Figure 1). The distance between each point is equal, around 1.5 mm.

Considering the test cell itself, further information can be found in previous work [START_REF] Gascoin | High temperature and pressure reactive flows through porous media[END_REF].

Figure 1 should be placed here

II . 2. Test methodology

N-dodecane (VWR, Rectapur 99% of purity) is used for these experiments because it is known as a good kerosene surrogate. It has been extensively studied by the authors in numerous pyrolysis configurations [START_REF] Gascoin | Characterisation of coking activity during supercritical hydrocarbon pyrolysis[END_REF][START_REF] Gascoin | Etude et mesure de paramètres pertinents dans un écoulement réactif[END_REF][START_REF] Abraham | Etude et développement d'une méthode d'analyse par spectroscopie infrarouge appliquée à la pyrolyse d'hydrocarbures en conditions supercritiques et transitoires[END_REF]. Thus, comparison with existent experimental data can be easily made. A stainless steel sample (Ø = 30 mm, filtration diameter = 16 mm thickness of 3 mm) is preferred as a porous medium compared to composite material because it is well characterized. Its implication in coke enhancing was notably demonstrated by PRISME team in tubular reactor configuration [START_REF] Abraham | Etude et développement d'une méthode d'analyse par spectroscopie infrarouge appliquée à la pyrolyse d'hydrocarbures en conditions supercritiques et transitoires[END_REF]. Nevertheless, no data are available for permeation configuration so that these works will provide new information on stainless steel surface effect and be compared with precedent results. Its class of porosity number 3 (the socalled SS 3 , Figure 2) presents an overall porosity of 30 %, pores diameter of 4.1 mm and grain diameter of 14.1 µm [START_REF] Gascoin | High temperature and pressure reactive flows through porous media[END_REF]. SEM (Scanning Electron Microscope) images with EDS (Energy Dispersive Spectroscopy) analysis are conducted after the tests on the porous sample to determine the nature of the deposited carbon. Initial composition is given in Figure 2. ( )

t + t - = gas i liquid i i Y 1 Y Y (3)
The permeabilities of Darcy and Forchheimer are computed thanks to the pressure drop and to the flow rate monitored as a function of time, using the Eq. 2. Eq. 2 is rewritten in conformity with the international norm ISO 4022 (Eq. 4) and the fluid properties are the mean one between the inlet and outlet conditions.

F D K . V . K 1 V . . L P m r + = m D (4)
A linear relationship (Eq. 5) derived from previous experiments with stainless steel tubular reactors [START_REF] Gascoin | Characterisation of coking activity during supercritical hydrocarbon pyrolysis[END_REF] has been used in this work to determine the coke formation as a function of time on the basis of methane production which is monitored by FTIR spectrometer at the process

III.! Results and discussion

III . 1. Distribution of pyrolysis product yields

The pyrolysis tests are conducted under the stationary thermal conditions described previously. The samplings are done for fluid temperatures of 725 K, 765 K and 820 K. For points A and B (upstream and inside the porous media), the samplings were only made for the two last temperatures because negligible pyrolysis was observed for the first temperature.

Table 1 presents the mass composition of the pyrolysis mixtures . . The authors precise that the precision is voluntary set at one hundredth in order to avoid the suppression of species which were detected in very low quantity. Nevertheless, the quantity lower than 0.1 wt.% should be considered for illustrative purpose mainly. First, it can be noted that the decomposition of dodecane varies not only with temperature but also as a function of the location in the porous medium (points A to C). The conversion rate increases along the flow.

This result is in qualitative agreement to that obtained under plug flow reactor on COMPARER project [START_REF] Gascoin | Etude et mesure de paramètres pertinents dans un écoulement réactif[END_REF][START_REF] Abraham | Etude et développement d'une méthode d'analyse par spectroscopie infrarouge appliquée à la pyrolyse d'hydrocarbures en conditions supercritiques et transitoires[END_REF]. This is directly related to the residence time distribution and to the convective heat transfer. Considering the pyrolysis products (Table 1), n-dodecane is preferentially decomposed into light species at the entrance of the cell and at the process outlet (respective predominant species: butane / ethane at 765 K and ethane / propylene at 820 K). Alkanes are found in majority compared to alkenes, contrary to what was observed generally in plug flow reactor during previous experiments [START_REF] Gascoin | Etude et mesure de paramètres pertinents dans un écoulement réactif[END_REF][START_REF] Abraham | Etude et développement d'une méthode d'analyse par spectroscopie infrarouge appliquée à la pyrolyse d'hydrocarbures en conditions supercritiques et transitoires[END_REF]. This is due to the heat transfers, which are enhanced in the porous material due to the high contact surface. The flow behaviour in porous medium is closer to the one in perfectly stirred reactor (homogeneous A c c e p t e d M a n u s c r i p t 10 flow) than to the plug flow reactor (spatial heterogeneity). The H over C atoms mass ratio tends to one sixth for all the alkenes; it is one twelfth for benzene and it starts from one third for methane to around 0.181 for dodecane. It was found in previous work [START_REF] Gascoin | Etude et mesure de paramètres pertinents dans un écoulement réactif[END_REF] that the H over C ratio increases (higher than 0.18) when the heat transfers increase, which indicates a higher alkanes production. The present results (Table 1) thus confirm this observation. Indeed, a higher proportion of ethane than ethylene or propane than propylene are found whatever the temperature or the axial position, with a H over C ration always higher than 0.18. Additional comments will be provided on this point at the end of this section with Figure 3.

Table 1 should be placed here

A different repartition is found (Table 1) inside the porous medium (preponderant midweigh species: hept-1-ene 765 K and hex-1-ene and butane at 820 K). Several assumptions can be made to explain these facts. Second, this could be attributed to the chemical diffusion of the species in the porous media. The lack of light species at point B is due to their higher filtration speed. Nevertheless, under steady-state conditions, this effect should be low. Some numerical computations have been done for another project and only small discrepancies in similar conditions have been found between simulations with and without chemical diffusion. Finally, the differences can be related to the coking activity which modifies the surface effect (heterogeneous reactions) and the heat and mass transfers. In case of deactivation of the surface or of the saturation of the sites, a limitation of the surface catalytic effect could be encountered. Such dynamic was already found considering other fuel [START_REF] Ghaly | Reactor Surface Effects During Propylene Pyrolysis, in: Industrial and Laboratory Pyrolyses[END_REF]. Thus, this assumption must be correlated with coke analysis and to the formation and consumption of some species (acetylene, ethylene, benzene). The benzene, toluene and p-xylene formation increases for a temperature of 765 K at point B, and then decreases at 820 K (Table 1). This can sign the coke formation (in a first time, production of coke precursors, and then an aggregation of those ones to form solid carbon). The heat transfers are reduced and the pyrolysis is slowed down (less formation of light compounds). On the contrary, it can also be assumed that the light species, precursors of coke, begin to react and form the secondary species involved in the formation of carbon deposit (decrease also shown at point B, Table 1).

As a consequence, this analysis remains difficult on the basis of these first results and, again, numerical simulation would be of great interest. However, the amplitudes of the variations are quite small and the accuracy of the measures could be responsible for these changes regarding benzene, toluene and p-xylene for example.

The chemical spatial profiles of group species can be further discussed. It is clear that the part, where chemical reactions are the most intense, is the second half of the porous medium.

Significant fall of light species mass fractions (C 1 to C 3 ) can be seen at the point B for the two temperature plateaus. Regarding mid-weigh species (C 4 -C 6 ) the production is quite similar Due to the very low mass fractions, the interpretation of such results remains very difficult.

Nevertheless it could be assumed that this profile is due to the chemical induction delay which depends on the fluid temperature. But, this can not fully explain such an activity. Spatially, point A and B are very close so that the fluid temperature is identical. The presence of the porous medium may increase the residence time and, consequently, the decomposition rate.

The high difference between the profiles may also be linked to a catalytic activity of the porous medium. To investigate the intrinsic effect of the porous media, a comparison can be made with previous COMPARER data obtained under open tubular reactor configurations [START_REF] Abraham | Etude et développement d'une méthode d'analyse par spectroscopie infrarouge appliquée à la pyrolyse d'hydrocarbures en conditions supercritiques et transitoires[END_REF]. The comparison is conducted for the samples issued from the point C. For the plug flow reactor (pressure = 10 -60 bar, residence time = 50 -100 s and mass flow rate = 50 -100 mg.s -1 ), the dodecane decomposition rate is lower than 10 wt.% for a furnace temperature close to 850 K instead of about 48 wt.% for the porous sample (Table 2). Indeed, even by considering a lower pressure (10 bar), a lower mass flow rate (50 mg.s -1 ) and consequently a theoretical higher decomposition, the pyrolysis rate is more than five times lower for plug flow experiment.

Table 2 should be placed here.

To reach a comparable decomposition, the considered temperature must be far higher (over 1250 K, Table 3). The species repartition is quite different: more light and mid-weighted species are found for plug flow reactor experiments while heavy species and, in particular, coke precursors are more present in the permeation configuration (Table 3). So, the increase of the dodecane decomposition and its different profile compared with plug flow reactor may be explained by the presence of the porous medium (i.e. increase of the residence time, of the contact surface and of the catalytic effect).

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

13 Table 3 should be placed here.

Finally, the gasification and pyrolysis rate have been plotted for the three temperatures and three sampling locations (points A to C) (Figure 3). A qualitative parabolic trend corresponding to the plug flow reactor is given from the COMPARER project results (full black curve, Figure 3). The experimental data points obtained for the porous flow are close to the line. This line (from [0%;0%] to [100%;100%]) means that 100 % of the pyrolysis products are found under gas phase at ambient conditions. This clearly highlights the difference between the two types of flow. The strong gas formation, in case of porous flow, is mainly due to the higher contact surface between the fluid and the overall reactor because of the presence of the porous medium. The residence time also explains the differences but it is not possible to estimate the one in porous media for compressible reactive fluid without numerical simulation support.

Additional computations were done for the same operating conditions than the experiments but using a perfectly stirred reactor configuration with the n-dodecane pyrolysis mechanism furnished by Dahm [START_REF] Dahm | Experimental and modelling investigation of the thermal decomposition of ndodecane[END_REF]. The results show strong similarities between PSR and porous flow compositions while PSR computations are not able to represent the multiphysic coupled phenomena involved in open tubular reactor [START_REF] Gascoin | Etude et mesure de paramètres pertinents dans un écoulement réactif[END_REF]. Consequently, the permeation through porous media seems to present more similarities with PSR like conditions than with plug flow ones.

Figure 3 should be placed here

III . 2. Coking activity

At the end of the pyrolysis tests, the central zone on the porous sample -which corresponds to the fluid flow cross-section-is clearly covered on the side of point A by a black layer of about 0.3 mm (Figure 4a). No carbon deposit is found on the opposite side (point C). This coke formation has been observed by SEM and an accumulation of solid micrometric particles A c c e p t e d M a n u s c r i p t 14 is found. and which mean diameter is of the order of 1 mm (Figure 4b). After surface cleaning, another SEM images clearly shows a clogging of the pores by similar particles with much smaller diameter (of the order of 0.1 mm or less) (Figure 4c). The particles of the permeable sample which are visible on the Figure 4c are covered by a carbon layer which composition is not so different from the one of the solid coke particles. Indeed, the coke particles at the surface are almost constituted of pure carbon (96 wt.% with 4 wt.% of O atoms, Figure 4)

while the stainless steel particles reveal a higher carbon content (about 20 wt.% instead of 7 wt.% initially, Figure 2). These 20 wt.% are found with an increase of about 1 wt.% of O atoms (2.34 wt.% instead of 1.31 wt.% initially). The C/O ratio is thus about 24 for the solid coke particles and around 20 for the coke at the stainless steel particles surface (if this 1 wt.% is strictly attributed to the coke formation, which is questionable). This slight difference may be related to the origin of formation of the coke (catalytic, i.e. inside the porous medium, and pyrolytic, i.e. at the surface) as it has been studied in previous work [START_REF] Gascoin | Characterisation of coking activity during supercritical hydrocarbon pyrolysis[END_REF]. Nevertheless, the differences are too small to be quantitatively interpreted and these results should be completed in a future work to determine quantitatively the effect of stainless steel and of the contact surface area in the process. It can only be mentioned that the presence of carbon layer on the particles of stainless steel demonstrates the implication of the solid material in heterogeneous reactions.

Figure 4 should be placed here

In addition to this macroscopic quantification, the mass of carbon deposit can be approximated as a function of time thanks to a previous experimental relationship as described in section II.B. The methane formation has been monitored by FTIR measures during the test and the estimated coke mass is given with and without its normalisation to the final mass measured at the end of the test (Figure 5). Indeed, directly using the relationship as mentioned in section II.B gives a total amount of coke of 3 g at the end of the experiment A c c e p t e d M a n u s c r i p t while only 0.46 g are detected. This over-prediction can be related to the contact surface area differences between the porous flow and the plug flow reactor. As a consequence, it is preferred to normalise the coke mass evolution to the maximum measured one, that is to say 0.46 g. This carbon deposit formation is also given as a function of temperature (Figure 6).

Both parameters are roughly related except when the temperature reaches a thermal plateau (725 K, 765 K and 820 K). The mass of coke increases during these plateaus because of the time of experiment. Figure 6 should be placed here.

III . 3. Clogging effect of the coke formation

The formation of coke deposit inside the pores of the permeable material and on its surface modifies the filtration of the fluid. During the initial heating of the system from ambient to 700 K, before the n-dodecane starts to pyrolyse, a constant value of 1.96.10 -13 m² is found for the Darcy's permeability. This is in complete agreement with the value found at ambient condition with a nitrogen flow [START_REF] Gascoin | High temperature and pressure reactive flows through porous media[END_REF]. Then, the transient darcian value computed during the experiment is given in Figure 7. The vertical straight lines represent the sampling times at the point A. Some decreases are found up to 2 to 3 orders of magnitude, which signify strong clogging of the porous medium. In addition, several sudden increases are observed at the locations marked by the vertical lines. On the basis of Figure 7, it is assumed that making a sampling upstream the porous medium impacts the flow. When opening the sampling line to the cylinder, the inlet pressure briefly becomes lower than the outlet one so that the flow is reversed (because of a lower pressure in the cylinders than in the process). The porous medium is temporarily cleaned and solid particles are removed. As a result, the Darcy's permeability increases.

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

The effect of the coke formation on the Darcy's permeability is mainly due to the accumulation of carbon particles inside the pores of the stainless steel sample because the permeability of the coke deposit on a surface (as observed on point A) has been determined

previously to be about 10 -8 m² [START_REF] Gascoin | Permeation of inert and supercritical reactive fluids through metallic and composite media[END_REF]. Thus, this high value (in comparison to the original permeability of the medium) should not impact the measure of the K D . This conclusion is reinforced by the determination of mass of carbon deposit which is detected at the surface (8.66 mg) compared to the one inside the sample (460 mg). The coke quantity is 50 times more important in the medium than at its surface. Considering the density of solid carbon, this corresponds to a volume of about 0.230 cm 3 , which can be compared to the volume of pore in case of a porosity of 30 %. (about 0.645 cm 3 ). Thus, the coke limits this volume by about 30 % which may have an important impact on the filtration (Figure 7).

Figure 7 should be placed here

The surface weight information can also be used to determine the coking rate. To provide a relationship between permeability variations and chemical phenomena, the coke formation is a good witness of the process (Figure 8). The coking rate is calculated using the following equation:

r C t S m C = & (6)
Where the residence time is varying during the experiment and it is estimated by:

fluid r m v t & ŕ = (7)
Two qualitative similar behaviors are found for the coking rate and for the rate of change in Darcian permeability (Figure 8). The production of coke limits directly the filtration surface. This has an impact on heat and mass transfer because the residence time decreases and it should lead to a lower decomposition. Nevertheless, such comportment is not clearly understood and further work is requirement to confirm these preliminary results. This relationship between coke formation and its modification of the permeability is of great interest in case of active cooling in order to predict the modification of cooling efficiency. In addition, a time delay may appear between both curves because the arrows are not all vertical and their slope is decreasing, which means that the clogging process slows down but the cleaning effect of the sampling at point A modifies the dynamics of the system and this conclusion remains questionable. ~~0,00 0,00 0,00 0,00 0,00 0,00 0,00 p-Xylene ~~0,00 0,00 0,01 0,00 0,02 0,00 0,03 Non-1-ene ~~0,10 0,19 0,15 0,20 0,07 0,09 0,56 Ethylbenzene ~~0,00 0,00 0,00 0,00 0,00 0,00 0,03 Nonane ~~0,11 0,09 0,04 0,14 0,00 0,15 0,26 Dec-1-ene ~~0,14 0,00 0,00 0,31 0,00 0,00 1,05 Decane ~~0,00 0,07 0,05 0,07 0,10 0,06 0,00 Undec-1-ene ~~0,04 0,07 0,01 0,07 0,00 0,00 0,14 Undecane ~~0,09 0,19 0,13 0,17 1,53 0,15 0,19 Dodec-1-ene ~~0,00 0,00 0,00 0,00 0,00 0,00 0,

C 1 =

 1 numerical constant (2096.3) µ = dynamic viscosity (Pa.s -1 ) C 2 = numerical constant (2193.7) ρ = density (kg.m -3 ) C & = coking rate (g.m -2 .s -1 ) τ = gasification rate (wt. %) g = Earth's gravity (m².s -1 the flow line (m 3 ) solid = superscript for solid phase V = fluid velocity (m.s -1 ) total = superscript for all phases Y i = i specie mass fraction (wt. %) A c c e p t e d M a n u s c r i p t I.! Introduction High speed vehicles are facing large thermal load because of the high velocity and because

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

 6

A c c e p t e d M a n u s c r i p tFigure 2

 2 Figure 2 should be placed here

A c c e p t e d M a n u s c r i p t 9 outlet.

 9 This law will be used later in the paper to analyze and to explain the variations of Darcy's permeability of the material during the coking activity.

First, the light

  compounds are sampled preferentially due to a bad sampling technique and the absence of heavy species at the inlet is the consequence. Nevertheless, this is not the most possible reason because heavy compounds are found for the samples related to point B with the same sampling technique. It is thus supposed that the sampling technique provides reliable samples, which are representative of the test conditions. It must be noted that samples at the outlet are satisfactory because they are obtained with a validated methodology which is used since 2005 (6 years of intensive pyrolysis work).

A c c e p t e d M a n u s c r i p t 11 Third

 11 , the radial heterogeneity of the thermal, hydraulic and chemical parameters implies that the samples are mostly representatives of the specific locations where they are made and not of the corresponding longitudinal coordinate. This is a possibility which can not be easily investigated by experimental test. For this reason, numerical simulation should be conducted in a near future.

A c c e p t e d M a n u s c r i p t 12 even

 12 if the drop between A and B is less important. The profile of heavy precursors of coke is comparable for a temperature of 820 K but is clearly different for the 765 K thermal stage.

Figure 5

 5 Figure 5 should be placed here.

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

Figure 8

 8 Figure 8 should be placed here

u s c r i p t 22 Table 1 .

 221 Pyrolysis products distribution and paramaters (in wt.%) at locations A (upstream), B (inside) and C (downstream) analyzed by GC/MS as a function of the maximum fluid temperature.

Figure 1 .

 1 Figure 1. Schematic of the apparatus to be inserted in the furnace with representation of the

Figure 2 .

 2 Figure 2. a) Stainless steel porous medium: Fe (62 wt.%), Cr (18 wt.%), Ni (11 wt.%)

Figure 3 .

 3 Figure 3. Type flow effect on the relationship between gasification and pyrolysis rate.

Figure 4 .

 4 Figure 4. a) Coked porous medium. b) Visualization by SEM (5000 times magnification) and

Figure 5 .

 5 Figure 5. Time evolution of the mass of coke deposit with the one of methane production and

Figure 6 .

 6 Figure 6. Relationship between coke formation and fluid temperature.

Figure 7 .

 7 Figure 7. Variation of the Darcian term and estimated mass of coke.

Figure 8 .

 8 Figure 8. Relationship between coking activity and Darcian permeability.

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

Table 2 .

 2 Comparison of pyrolysis parameters (in wt.%) with or without porous medium

	considering similar conditions
	Table 3. Comparison of pyrolysis products (in wt.%) with or without porous medium
	considering similar decomposition rates

  Through porous mediumTubular reactor 1328 K, 10 bar, 50 mg 1359 K, 10 bar, 100 mg 850 K, 35 bar, 65 mg 1280 K, 60 bar, 100 mg

	Type of reactor	Through porous medium Tubular reactor Tubular reactor	Tubular reactor
	Experimental Hydrogen Product conditions Methane	725 K 914K, 10 bar, 0,00 896 K, 10 bar, 50 mg/s 100 mg/s 4,91	850 K, 35 bar, 0,03 65 mg/s 2,32	765 K 916 K, 60 bar, 100 mg/s	0,04 2,16	820 K	0,01 1,18
	Dodecane Ethylene	91,61	13,06	99,58	51,57 4,68	96,44		3,84	1,03
	Gasification rate Ethane	7,40		6,90	0,00	39,00 2,60	1,00		6,70	1,79
	Pyrolysis rate Propylene	8,39		14,47	0,42	48,43 4,69	0,04		7,84	2,49
	Propane			1,65			1,68			7,17	2,34
	But-1-ene			0,00			3,79			4,45	1,79
	Butane			4,80			1,07			3,72	1,68
	(E)-Pent-2-ene		0,00			2,62			2,05	2,22
	Pentane			3,00			4,53			1,80	3,33
	Hex-1-ene			0,06			1,43			1,95	2,98
	Dodecane Tridec-1-ene Tridecane Tetradecane Hexane Benzene Hept-1-ene Heptane Oct-1-ene Octane Toluene Non-1-ene Nonane Dec-1-ene Decane Undec-1-ene Undecane Dodecene Dodecane Gasification rate > 99 ~~0,00 > 99 83,84 ~~0,01 ~~0,01 0,04 0,00 0,00 0,01 0,06 0,01 0,00 0,37 0,04 1,42 0,11 0,41 0,00 0,00 48,67 41,00 Pyrololysis rate 51,33	90,25 0,01 0,03 0,02	33 51,57 0,00 0,04 0,02 a n u s c r i p t 71,04 83,31 80,92 0,01 0,00 0,00 0,03 0,00 0,02 0,02 0,00 0,02 M 86,70 0,00 0,02 0,01 M 0,58 1,40 0,08 0,08 1,29 1,74 0,54 0,42 1,05 0,96 0,27 0,39 0,18 0,13 0,56 0,61 0,26 0,41 1,05 0,69 0,14 0,16 0,19 0,20 0,33 0,36 0,00 0,06 67,71 0,28 0,00 0,00 0,02 0,08 0,01 0,00 0,43 0,07 1,24 0,14 0,39 0,14 0,06 51,57 71,26 25,50 39,00 15,00 a n u s c r i p t 32,29 48,43 28,74
	Gasification rate	0,00	0,00	13,50	6,00		10,00	26,00	10,00	15,00	39,00
	Pyrolysis rate H/C ratio A c c e p t e d < 1 < 1 16,16 9,75 13,30 --0,2011 0,2011 0,1846 A c c e p t e d	28,96 0,2050	17,00 0,2025	19,00 0,1916	48,43 0,1982
								1

Acknowledgments

This work was supported by the ESA-ESTEC, Contract no.: 3-12861/09/NL/PA. The authors would like to sincerely thank B Le Naour and J Bertrand from MBDA-France for their help involving this project. The highly valuable work of D. Blanc and D. Courilleau involving the measurement cell of PRISME bench has been greatly appreciated. The authors also would like to thank the help of J.-C. Hargé for the SEM and EDS measurements.