Bérénice Froment
email: fromentb@obs.ujf-grenoble.fr

Michel Campillo

Philippe Roux

Reconstructing the Green's function through iteration of correlations Reconstruction de la fonction de Green par itération de corrélations

Keywords: ambient-noise data, cross-correlations, C 3 method, passive tomography Mots-clés : données de bruit sismique, corrélations, méthode C 3, tomographie passive

Correlations of ambient seismic noise are now widely used to retrieve the Earth response between two points. In this study, we reconstruct the surface-wave Green's function by iterating the correlation process over the tail of the noise-based correlation function. It has been demonstrated that the so-called C 3 function shows the surface-wave part of the Green's function. Using data from 150 continuously recording stations in Europe, the C 3 results help in the extraction of the travel-times from noise-based measurements, especially through the suppression of effects caused by non-isotropic source distributions. We present the results of the next iterative step (i.e. C 5 ), which shows that some coherent signal is still present in the coda of the C 3 function, and we investigate the evolution of the reconstruction of the Green's function throughout the iteration process. Finally, we discuss the interest for combining information from the different correlation functions to improve noise-based tomography analysis.

Résumé

Les corrélations de bruit sont aujourd'hui largement utilisées pour reconstruire la réponse impulsionnelle de la Terre entre deux points. Nous reconstruisons ici les ondes de surface de la fonction de Green en itérant le processus de corrélation sur la partie tardive des corrélations de bruit. Il a été montré précédemment que le mode fondamental des ondes de surface était clairement visible dans la fonction de corrélation C 3 ainsi obtenue. Les résultats obtenus sur les données de 150 stations en Europe montrent que cette méthode peut aider à améliorer les mesures de temps dans les corrélations, en particulier en supprimant les effets liés à une distribution de sources non isotrope. Nous présentons également les résultats de l'itération suivante (appelée C 5 ) montrant encore du signal cohérent et nous étudions l'évolution de la reconstruction de la fonction de Green lors de l'itération. Enfin, nous discutons l'intérêt de combiner l'information de ces différentes fonctions de corrélation.

Introduction

It has been shown that correlating long time series of seismic ambient noise recorded at two stations makes it possible to retrieve the Earth response (i.e. the Green's function) between these two stations (e.g. [START_REF] Campillo | Phase and correlation in 'random' seismic fields and the reconstruction of the Green function[END_REF][START_REF] Gouédard | Cross-correlation of random fields: mathematical approach and applications[END_REF][START_REF] Larose | Correlation of random wavefields: an interdisciplinary review[END_REF]Sabra et al., 2005a;[START_REF] Shapiro | Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise[END_REF]. This technique has provided useful results for studies of tomographic imaging (e.g. Sabra et al., 2005b;[START_REF] Shapiro | High-resolution surfacewave tomography from ambient seismic noise[END_REF][START_REF] Yang | Ambient noise Rayleigh wave tomography across Europe[END_REF][START_REF] Yao | Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis -I. Phase velocity maps[END_REF][START_REF] Yao | Surface wave array tomography in SE Tibet from ambient seismic noise and two-station analysis -II. Crustal and upper-mantle structure[END_REF] and monitoring (Brenguier et al., 2008a, b;Sens-Schönfelder and Wegler, 2006;Wegler and Sens-Schönfelder, 2007).

In practice, since the ambient seismic noise is dominated by the surface-wave contributions, only the direct arrivals of surface waves are clearly visible in the noisecorrelation signals. However, theory says that in the case of a perfectly isotropic noise field, the noise correlation provides the complete Green's function , including all types of waves that propagate between the two stations in x A and x B (e.g. [START_REF] Gouédard | Cross-correlation of random fields: mathematical approach and applications[END_REF][START_REF] Lobkis | On the emergence of the Green's function in the correlations of a diffuse field[END_REF][START_REF] Roux | Ambient noise cross correlation in free space: Theoretical approach[END_REF][START_REF] Sanchez-Sesma | Elastodynamic 2D Green function retrieval from cross-correlation: Canonical inclusion problem[END_REF][START_REF] Snieder | Extracting the Green's function of attenuating heterogeneous acoustic media from uncorrelated waves[END_REF][START_REF] Wapenaar | Retrieving the elastodynamic Green's function of an arbitrary inhomogeneous medium by cross correlation[END_REF]. If we let u(x A ,t) and u(x B ,t) denote time-dependent random wavefields recorded at two sensors, A and B, we can write:

(1) with , i.e. the correlation over a time interval T. We introduce as the correlation when . A theoretical guess for the reconstruction of the Earth response is:

(2)

where, G + and G -represent the causal and anticausal Green's function, respectively. Note that in practical applications, most studies choose to deal with a normalized correlation function:

(3),

for which the denominator corresponds to the square root of the total energy of the signals recorded in x A and x B . In this case, the maximum of the normalized correlation function corresponds to the coherence level between the two stations.

In the ideal case, the later part of the noise correlation should contain the coda part of the Green's function, and it might be possible to correlate these waves as we correlate standard earthquake codas, to reconstruct the Green's function [START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF]. This is the concept of the so-called C 3 method, i.e. re-correlating the coda of the noisecorrelation functions to reconstruct the Green's function [START_REF] Stehly | Reconstructing Green's function by correlation of the coda of the correlation (C 3 ) of ambient seismic noise[END_REF].

For microseism records dominated by surface noise, the required spatial isotropy of the recorded wavefield is mainly produced by the distribution of sources and by the scattering of seismic waves in the medium. However, usually the distribution of the noise sources does not provide a perfect random noise field. In this case, the theoretical requirements are not completely fulfilled, and the noise correlation does not perfectly reconstruct the surface-wave Green's function, as some fluctuations remain in the correlation results. This, in practice, prevents the identification of other amplitude arrivals that are smaller than the predominant fundamental mode of the surface waves. This is the main limitation in the implementation of the C 3 method, i.e. does the later part of the noise correlation contain the coda part of the surface-wave Green's function? Note that this is also a central issue for noise-based temporal monitoring of seismic velocity, for which measurements are performed on noise correlation codas.

In their recent study, [START_REF] Stehly | Reconstructing Green's function by correlation of the coda of the correlation (C 3 ) of ambient seismic noise[END_REF] demonstrated successful results for the C 3 function, clearly showing the direct surface-wave part of the Green's function. The success of this technique is an indication that the later part of the ambient noise correlations is meaningful and contains multiply scattered waves with properties of regular earthquake codas. [START_REF] Stehly | Reconstructing Green's function by correlation of the coda of the correlation (C 3 ) of ambient seismic noise[END_REF] showed that the C 3 method helps to improve the time symmetry of correlation functions, thereby increasing the quality of the data used in noise-based tomographic studies. Indeed, since the noise correlation functions should converge to the causal and anticausal Green's function (Eq. ( 2)), time symmetry in noise correlations may constitute a way of controlling the quality of the Green's function reconstruction from the correlation process. In particular, travel-time estimations must be the same in both causal and anticausal parts. That is why some tomographic studies use only time-symmetric noise correlations with good signal-to-noise ratios. Since [START_REF] Stehly | Reconstructing Green's function by correlation of the coda of the correlation (C 3 ) of ambient seismic noise[END_REF], some studies have proposed a formalization of the C 3 method [START_REF] De Ridder | Kinematics of iterative interferometry in a passive seismic experiment[END_REF][START_REF] Garnier | Passive sensor imaging using cross-correlations of noisy signals in a scattering medium[END_REF]. In particular, [START_REF] Garnier | Passive sensor imaging using cross-correlations of noisy signals in a scattering medium[END_REF] proved the validity of this method based upon the stationary-phase analysis of the C 3 function leading terms, and they confirmed that this method can enhance the quality of travel-time estimates in the case of anisotropic source distributions.

As indicated above, the noise-correlation technique has provided useful information for seismological applications, and this in spite of the imperfect reconstruction of the Green's function in correlation functions. In particular, an anisotropic source distribution leads to biased travel-times extracted from correlation functions [START_REF] Froment | Estimation of the effect of nonisotropically distributed energy on the apparent arrival time in correlations[END_REF][START_REF] Weaver | On the correlation of non-isotropically distributed ballistic scalar diffuse waves[END_REF]. One of the main issues in the improvement of noise-based measurements is thus to understand and eliminate these source effects in correlation results. In this study, we show the improvement of noise-based results in seismological applications through the use of the C 3 method.

In section 1, we present the data and the processing used in this study to compute the C 3 function. In section 2, we examine the time symmetry of the C 3 function and its relationship with the station distribution in the network. Then, we discuss the interest of combining the information from the different correlation functions. Finally, we show the correlation function that is obtained by iterating the correlation process on the C 3 coda. In particular, we investigate the evolution of the reconstruction of the Green's function throughout the iteration process.

Data and C 3 computation

Throughout this study, we have used pairs of stations that are located in the Alps as part of a regional network of 150 broad-band stations. We consider the vertical motions of seismic ambient noise recorded continuously over one year. The noise processing is the same as in [START_REF] Stehly | Reconstructing Green's function by correlation of the coda of the correlation (C 3 ) of ambient seismic noise[END_REF].

For each station pair, we correlate the noise records, which are pre-filtered in the 0.1-0.2-Hz frequency band, and the noise correlation (C 1 ) is obtained by stacking all of the dayby-day correlations. We then apply the C 3 method between two stations A and B, according to the following successive steps:

1. We take a third station, S, in the network located at x S , and we consider the result of the A-S and B-S noise correlations (C 1 ) as the signal recorded at A and B, respectively, if a source was present at x S . The receiver in S thus has the role of a virtual source.

2. We then select a 1200-s-long time window (T coda ) in the later part of the noise correlation (C 1 ), beginning at twice the Rayleigh wave travel-time (i.e., in the same way as [START_REF] Stehly | Reconstructing Green's function by correlation of the coda of the correlation (C 3 ) of ambient seismic noise[END_REF]). Note that the coda in both the positive and negative parts of the noise correlation (C 1+

and C 1-) are selected.

3. We cross-correlate the codas in the A-S and B-S noise correlations. Note that we whiten the codas between 0.1 and 0.2 Hz before correlation; the motivation for this processing is discussed later. We compute the correlation functions between the positive-time and negativetime codas to finally obtain two correlations that we denote as C 3++ S and C 3-- S , respectively.

4. We average over the N stations of the seismic network to obtain C 3++ and C 3--:

(4)

(5) 5. We stack C 3++ and C 3--to obtain the actual C 3 correlation function between A and B:

(6)

The steps described above highlight the difference in the correlation technique between noise correlations (C 1 ) and the C 3 method, that points out advantages of this iterative method. For instance, since the C 3 method correlates signal that has been extracted from the C 1 computation, it is expected to be more coherent than seismic ambient noise. Furthermore, we can distinguish coda and ballistic contributions in C 1 functions (which is not possible in seismic noise); this provides scattering averaging in the correlation process when only coda waves are selected. Finally, we point out that, by construction, the virtual sources S are uncorrelated in the C 3 method, which may not be the case for noise sources in C 1 computation.

Figure 1 illustrates the C 1 and C 3 correlation functions for one station pair. The good agreement between the pulse shapes in the C 3 and C 1 correlation functions shows that the C 3 correlation process has extracted from the coda of the C 1 correlation the necessary phase information to reconstruct the direct surface wave part of the Green's function.

Network stations as virtual sources

Source signature in correlation functions

In the case of noise sources that are homogeneously distributed, the positive and negative parts of the noise-correlation functions are expected to be symmetric, as they represent the causal and anticausal parts of the Green's function, respectively. In practice, these two parts may not have the same amplitude, which reflects a difference in the energy flux that propagates between the two stations in both directions [START_REF] Paul | Empirical synthesis of time-asymmetrical Green functions from the correlation of coda waves[END_REF][START_REF] Van Tiggelen | Green function retrieval and time reversal in a disordered world[END_REF].

Thus, if the source density (or virtual source power) is larger for one direction than for the other, a time asymmetry will be seen for the correlation function (see Figure 2). This means that the level of symmetry in the correlations, up to first order, is an indication of the source distribution.

What makes C 3 appealing is that, in theory, the "noise" sources that contribute to C 1 and C 3 are different. Considering the two stations A and B, the noise correlation (C 1 ) is directly affected by the noise-source distribution. In contrast, the C 3 method is based on the use of the network stations as virtual sources [START_REF] Stehly | Reconstructing Green's function by correlation of the coda of the correlation (C 3 ) of ambient seismic noise[END_REF]. In this case, the symmetry in the C 3 function should no longer reflect the noise-source distribution, but the distribution of the network stations around the path A-B.

Time symmetry of C 3 and C 1

To verify these theoretical expectations, we compare the time-symmetry of both correlation functions (C 1 and C 3 ) for the particular case of the station pair HAU-BOURR (Figure 3). In the C 3 correlation process, we chose to select only the 55 network stations south of HAU-BOURR (see map in Figure 3a). This means that the theoretical C 3 virtual sources are localized south of the station pair of interest and opposite to the incident-noise direction, as the main source of noise is located in the north Atlantic ocean or on the northern European coast (e.g. [START_REF] Friedrich | Ocean-generated microseismic noise located with the Grafenberg array[END_REF][START_REF] Kedar | The origin of deep ocean microseisms in the north Atlantic Ocean[END_REF]Landès et al., 2010;[START_REF] Stehly | A study of the seismic noise from its longrange correlation properties[END_REF].

Figure 3b shows the results of both the C 1 and C 3 correlation functions for this geometry. The time symmetry of both signals shows two opposite peaks of maximum amplitude, which reflect a main energy flux propagating in opposite directions: from North to South for C 1 , and from South to North for C 3 . Note that the same convention is used for all of the correlation functions, i.e. the signal in negative (resp. positive) correlation times corresponds to waves propagating from North to South (resp. South to North). This observation indicates that the time symmetry of the two correlation functions is consistent with the location of their expected contributing source distributions.

However, the C 3 function has a better time symmetry than C 1 . This indicates that the coda waves used to compute the C 3 function constitute a more isotropic field than the ambient seismic noise. This can be explained by the nature of the coda waves, which correspond to waves scattered from the heterogeneities in the Earth subsurface (e.g. [START_REF] Aki | Origin of coda waves -source, attenuation, and scattering effects[END_REF].

This scattering process is expected to make the coda wavefield more isotropic, depending on the distribution of the scatterers in the medium and the lapse-time considered in the coda [START_REF] Paul | Empirical synthesis of time-asymmetrical Green functions from the correlation of coda waves[END_REF]. To confirm the influence of the scattering in our observations, we compare the C 3 results for the same station pair (HAU-BOURR) computed for two different time windows T, in the C 1 correlation (Figure 4). In the first case, we consider a 1200-s time window in C 1 , including the direct surface waves (see Figure 4a). This time window begins 15 s before the Rayleigh wave travel-time (i.e., twice the dominant period of the signals). We denote it as T all , as it takes into account both the ballistic and the scattered waves. Due to the strong amplitude of the direct waves in the time window, we expect that their contribution will dominate in the C 3 process. The resulting C 3 all function is shown in Figure 4b. In the second case, we use C 3 coda for the C 3 function calculated (as in section 1) from a 1200-s time window T coda of the correlation located entirely after the direct arrivals (Figure 4c).

In view of the nature of the dominant waves involved in both of these computations, we can infer the origin of the energy flux observed in C 3 all and C 3 coda . The energy flux is expected to be totally controlled by the distribution of the virtual sources in the first case, whereas it must also be influenced by scattering (and the distribution of the scatterers) in the second case. The station distribution considered in this case helps in the separation of the different contributions to the time symmetry in the correlation process.

On the one hand, Figure 4b confirms the role of network stations as virtual sources in the C 3 method, as C 3 all shows complete time asymmetry, which reflects a directive energy flux in the opposite direction to the seismic noise direction that is consistent with the distribution of the network stations. On the other hand, the time symmetry is preserved in 4c), which indicates the link between the scattering events and the wavefield isotropy in the correlation process. Finally, we notice the lower coherence observed in C 3 coda compared to C 3 all . This likely reflects the higher level of fluctuations in the coda of C 1 than in its ballistic part. Note that the coherence in C 3 all is about twice as high as that in C 1 in this geometry particularly favorable to C 3 all (i.e., a dense distribution of virtual sources in the alignment of receivers).

C 3 coda (Figure
This analysis shows that the time symmetry of the C 3 function does not depend any more on the distribution of the noise sources, but is instead essentially controlled by the distribution of the stations in the network when direct waves are selected, and is also influenced by the contribution of the scatterers when the C 3 process is performed from the coda. This means that the sources in the C 3 method are controllable to a certain extent, which is a noticeable advantage compared to the C 1 process, which usually suffers from uneven noise-source distribution.

Combining information from the different correlation functions

The C 3 method is an alternative method for the reconstruction from seismic noise of the Green's function between two points, and we have shown that it can improve noise-based measurements, especially by making correlation functions independent of the noise source.

But this iterative method can also be viewed as a way to obtain several estimates of the Green's function between two points that can be used to complement each other for tomography purposes. Indeed, noise-based tomography requires unbiased travel-time measurements between station pairs. In the case of a directive noise, it has been shown that the noise-correlation functions carry the footprint of the noise-source direction, which can lead to biased travel-time estimations. In general, when both the causal and anticausal parts of the noise-based correlation functions are reconstructed, the travel-time estimation is said to be unbiased. Further bias might come from a poor signal-to-noise ratio in the correlation functions, which also prevents satisfactory measurements of the travel-time.

We show in Figure 4 for a given station pair, that the combination of the C 1 and C 3 results may alleviate these bias. In this case, we have seen that both the C 3 coda and C 3 all functions provide information with respect to C 1 since they exhibit also the time-reversed counterpart of the ballistic signal reconstructed in C 1 . Note that in the general case, this information is provided by C 3 coda , since the time symmetry of the C 3 all function is strongly controlled by the station distribution. There are therefore configurations for which the information in C 1 and C 3 all is redundant. On the other hand, in terms of coherence, C 1 and C 3 all can provide an advantage with respect to C 3 coda for better-quality measurements to be obtained. Besides, Figure 4 shows a higher level of coherence in C 3 all than in C 1 .

Several aspects have roles in the improvement of the correlation results. In the present study, we have mentioned in particular the noise-source effects, time symmetry and coherence. The C 3 functions complement C 1 by contributing to one or several of these aspects. The different correlation functions can thus be viewed as complementary information that can be considered together to improve noise-based tomography analysis.

Iterating the correlation process

The results presented in the previous sections have shown that it is possible to reconstruct (at least partially) the Green's function by iterating the correlation process from C1 to C 3 . These results emphasize that the late arrivals in the C 1 function still carry information about the propagation medium since the correlation of the C 1 coda averaged over the station networks allows the reconstruction of the direct arrivals.

We can then ask about the next iterative step, i.e. is it possible to use the coda of the C 3 function to reconstruct part of the Green's function?

For Figure 5, we computed the C 5 function 1 for a particular station pair from the coda of the C 3 functions averaged over the network station pairs. Note that in Figure 5 we compare correlation functions instead of Green's function. The relation between the Green's function and the correlation function (Eq. ( 2)) refers implicitly to a time derivative that is written as in the frequency domain. When dealing with correlation of correlations, an extra -term is introduced into the counterpart of Eq. ( 2), which becomes .

This implies a difference in the spectral content of the Green's function and the correlation function that becomes more pronounced during the iteration process. However, the spectral whitening used in our processing maintains the spectrum as almost constant in a narrow band during the iteration. This allows us to consider that the spectrum variations remain negligible throughout the iteration process.

The good agreement between the direct arrivals of the C 5 and the C 1 functions shows that some coherent signal is still present in the coda of the C 3 function. This result supports the idea that the C 5 function could also be used in noise-based tomography or for monitoring techniques. However, we observe in Figure 5 that the signal-to-noise ratio in the C 5 function is smaller than that in the C 3 function, which was also smaller than that in the C 1 function (note the decrease in the coherence level in Figures 5b, c andd). This might limit the practical interest of this iterative correlation process in geophysical inversion techniques. Sabra et al. (2005c) proposed a theoretical expression for the level of fluctuations in the C 1 function that helps in the understanding of the Green's function reconstruction during the iterative correlation process. They considered a white-noise model and a finite bandwidth B. In the case of normalized correlations (Eq. ( 3)), the expected level of fluctuations, n, is given by: (7).

The stack over N virtual sources, as performed in the C 3 method, leads to an expected fluctuation level in the C 3 function of (e.g. [START_REF] Larose | Fluctuations of correlations and Green's function reconstruction: role of scattering[END_REF][START_REF] Stehly | Reconstructing Green's function by correlation of the coda of the correlation (C 3 ) of ambient seismic noise[END_REF]:

. ( 8)

Measurements of the level of fluctuations carried out for the iterative correlation functions have shown good agreement with Sabra et al. (2005c). This means that the fluctuation levels observed (see Figure 5) are indeed controlled by the evolution of B, T and N during the iterative correlation process. In the following, we examine the evolution of these different parameters in this iteration process:

1. The frequency bandwidth B: In the frequency domain, the correlation consists of multiplying the spectra. If the frequency content is peaked, the correlation process heightens the contrast between the frequencies, which, in practice, amounts to having a narrower spectral content (B is decreasing). This corresponds to a loss of information, which in Eq. ( 7) leads to higher fluctuation levels in the correlations. In our case, the noise spectrum between 0.1 Hz and 0.2 Hz is peaked near 0.14 Hz, and this peak will be accentuated throughout the successive iterative steps. To prevent this loss of spectral information, we whiten the selected coda between 0.1 Hz and 0.2 Hz, so as to correlate a flat spectrum signal in the same frequency band at each iteration step.

2. The duration T: It is well established that correlation results are highly dependent on the duration of the ambient-noise time series. In the present study, C 1 has been calculated from one year of seismic ambient noise, and C 3 results from the average of about 200 correlations (C 3++ and C 3--averaged over a hundred stations) performed on a time window of T = 1200 s extracted from the C 1 coda. Considering Eq.( 8), this results in larger fluctuations in the C 3 function than in C 1 . As the fluctuation level is higher in C 3 than in C 1 , the C 5 function was computed from the selection of a shorter time window in the C 3 coda, which results in a stronger fluctuation level in the C 5 function. This explains why the level of fluctuations will always increase through the correlation iterations after the C 3 processing.

3. The number of virtual sources N: Eq. ( 8) gives a level of fluctuations that is decreasing as . This means that the stacking over many virtual sources improves the C 3 processing.

However, N remains the same after the C 3 function, and therefore this parameter does not influence the fluctuation level in the higher-order correlation functions.

This analysis shows that the level of fluctuations can be improved in the C 3 correlation (with respect to the C 1 correlation) by considering large networks (and thus more virtual sources), although it is limited by the decrease in the time interval T after the C 3 iteration.

Conclusions

The present study shows that the later part of the noise correlation contains coherent signals that can allow the reconstruction of the Green's function by the so-called C 3 method. Indeed, the C 3 correlation function clearly shows the direct surface-wave part of the Green's function, and the success of the next iterative step, which is C 5 , shows that coherent signals are also present in the coda of the C 3 function.

This method can be useful to improve noise-based measurements, especially by suppressing source effects that are caused by non-isotropic source distributions. Two points should be emphasized. First, the C 3 method is computed from scattered waves that have already been extracted in the C 1 function. Secondly, we have shown that the C 3 correlation function loses its dependency on the noise-source spatial distribution, and instead highly depends on the network stations that have the role of virtual sources. Thus, both scattering and an even spatial distribution of stations help in converging toward a two-sided Green's function with this technique. The C 3 method can therefore help in the resolution of problems that arise from uneven source distributions in the primary noise correlations (C 1 ).

However, fluctuations in C 3 coda are higher than in C 1 , and this might mitigate the practical interest in this higher-order correlation function for tomography purposes. We can adopt a somewhat different approach by combining information from the different correlation functions. Indeed, the C 1 and C 3 correlation functions provide several travel-time estimations that can be viewed as complementary information for noise-based tomography inversions.

Figure 1 :

 1 Figure 1: (a) Map of Europe with orange stars corresponding to about 100 stations of the network used as virtual sources to construct the C 3 function between stations EMV and

Figure 1 :

 1 Figure 1 : (a) carte de l'Europe avec les étoiles oranges représentant les stations du réseau (environ 100) utilisées comme sources virtuelles pour construire la fonction C 3 entre les

Figure 2 :

 2 Figure 2: Schematic illustration of the relationship between the distribution of sources and

Figure 2 :

 2 Figure 2 : représentation schématique de la relation entre la distribution des sources et la

Figure 3 :

 3 Figure 3: (a) Map of Europe with orange stars corresponding to the 55 stations of the network used as virtual sources to construct the C 3 function between stations HAU and

Figure 3 :

 3 Figure 3 : (a) carte de l'Europe avec les étoiles oranges représentant les 55 stations du réseau utilisées comme sources virtuelles pour construire la fonction C 3 entre les stations HAU et

Figure 4 :

 4 Figure 4: (a) Noise correlation C 1 between stations HAU and BOURR. The start of the two 1200-s time windows used in the C 3 processing are indicated: T all and T coda . For the sake of

Figure 4 :

 4 Figure 4 : (a) Corrélation de bruit C 1 entre les stations HAU et BOURR. Le début des deux fenêtres temporelles de 1200 s utilisées dans le calcul de C 3 (T all et T coda ) sont indiqués. Par

Figure 5 :

 5 Figure 5: (a) As for Figure 1a, here for the station pair PLONS-ZUR. (b) Noise correlation C 1 between stations PLONS and ZUR. (c, d) Iterated correlation functions C 3 (c) and C 5 (d)

Figure 5 :

 5 Figure 5 : (a) même que Figure 1a pour la paire de stations PLONS-ZUR. (b) Corrélation de bruit C 1 entre les stations PLONS et ZUR. (c, d) fonctions de corrélation itérées C 3 (c) et C 5

As for C 3 , the exponent indicates the number of C in the correlation-function description. The C 5 function corresponds to the Correlation of Coda of the C 3 function.

Acknowledgements

All the seismic data used in this study have been obtained from the IRIS DMC (http://www.iris.edu/), the ORFEUS database (http://www.orfeus-eu.org/), the ETH Zürich, the Commissariat à l'Energie Atomique (CEA, France), and the Laboratoire de Géophysique Interne et Tectonophysique, Grenoble (France). We thank N. Shapiro, U. Wegler and D. Draganov for their helpful comments and suggestions. One of the authors (B.F.) acknowledges the support of Shell Research. This research was supported by European Research Council (Advanced grant Whisper).


Sheng Di
email: sheng.di@inria.fr

Derrick Kondo
email: derrick.kondo@inria.fr

Walfredo Cirne
email: walfredo@google.com

Characterization and Comparison of Google Cloud Load versus Grids

A new era of Cloud Computing has emerged, but the characteristics of Cloud load in data centers is not perfectly clear. Yet this characterization is critical for the design of novel Cloud job and resource management systems. In this paper, we comprehensively characterize the job/task load and host load in a real-world production data center at Google Inc. We use a detailed trace of over 25 million tasks across over 12,500 hosts. We study the differences between a Google data center and other Grid/HPC systems, from the perspective of both work load (w.r.t. jobs and tasks) and host load (w.r.t. machines). In particular, we study the job length, job submission frequency, and the resource utilization of jobs in the different systems, and also investigate valuable statistics of machine's maximum load, queue state and relative usage levels, with different job priorities and resource attributes. We find that the Google data center exhibits finer resource allocation with respect to CPU and memory than that of Grid/HPC systems. Google jobs are always submitted with much higher frequency and they are much shorter than Grid jobs. As such, Google host load exhibits higher variance and noise.

I. INTRODUCTION Miron Livny once said, "I've been doing research in

Clouds before it was called Grids". This quote highlights fundamental questions: what is novel in Clouds, and what are the real differences between Clouds and Grids? Potential differences can arise in many areas, including but not limited to software configuration, resource and job management, and monetary charging models. High-level differences between Clouds and Grids are easy to identify. For instance, in Clouds, software configuration and management is eased via virtual machine images. Resource management is potentially more dynamic as virtual machines can be booted up and down over time. Charging models are pay-as-you-go.

While high-level differences may be easy to identify, the devil is in the details; substantiating such differences quantitatively and statistically, specifically for work load and host load, is critical for the design and implementation of Cloud systems. It is challenging because it requires detailed traces of production and often proprietary Clouds. Moreover, such traces are often enormous, and one requires statistical techniques for summarizing and understanding the data.

At a high-level, one could hypothesize fundamental differences in the time dynamics of load due to differences in the types of users and applications. In Grids and HPC systems, users often submit large scientific applications or experimental research applications. In Cloud data centers, users and applications are more commercial and interactive.

Typical jobs come from web services (for meeting search queries, and online document editing and translations, for instance), and map-reduce applications (for building inverted indices of web documents or email, for instance).

Given those hypotheses, we characterize quantitatively and statistically host and work load in commercial data centers, and compare and contrast properties of Cloud versus Grid load. By work load, we refer to load due to incoming jobs and their corresponding tasks submitted by users. By host load, we refer to load on a particular machine due to executing tasks. We use the sole term load to refer generically to both work load and host load.

A detailed characterization of a system's load is critical for optimizing its performance. In the context of data centers, host and work load characterization is essential for job and resource management including but not limited to capacity planning, virtual machine consolidation, load balancing, and task scheduling. For instance, by characterizing common modes of host load within a data center, a job scheduler can use this information for task allocation and improve utilization. Alternatively, the resource management system can proactively shift and consolidate load via (VM) migration to improve host utilization, using fewer machines and shutting off unneeded hosts.

So, we investigate the following questions. First, what are the characteristics and statistical properties of work and host load in Cloud data centers? Second, what are the key similarities or differences of work and host load between Clouds and Grids, and the implications for job/resource management?

To answer those questions, we study a workload trace of a production data center at Google [START_REF]Google cluster-usage traces[END_REF]. Google offers several Cloud services, such as web search, interactive documents, email, where users (indirectly) submit jobs without having to specify the details of their execution. As such, we believe this trace is representative of real-world Cloud workloads. The trace contains detailed measurements reported every 5 minutes of the states and events concerning users, jobs, tasks, and hosts. In total, the trace contains details about over 25 million tasks executed over 12,500 hosts during 1-month time period.

For comparison, we also characterize the work load and host load from different Grid systems, including AuverGrid [START_REF] Auvergrid | [END_REF], NorduGrid [3], sharcnet [START_REF]Sharcnet grid project[END_REF], as well as other HPC clusters, from Argonne National Laboratory (ANL) [START_REF]Anl cluster[END_REF], RIKEN Integrated Cluster of Clusters (RICC) [START_REF]Ricc cluster[END_REF], METACENTRUM cluster [7] and Lawrence Livermore National Laboratory (LLNL) [8].

We characterize the workload for Google's jobs, based on Google's trace data with over 670,000 jobs. By comparing to the real-world trace of many other Grid systems and HPC systems, we find that Google jobs usually have much shorter length are submitted at higher frequency. Each Google job takes relatively a small share of resources, leading to a much finer resource allocation granularity. Moreover, each Google job usually consists of only a single task; a typical user job, such as keyword search, is relatively compact and self-contained. In contrast, in Grids, each job can consist of several processes running simultaneoulsy over multiple cores.

With respect to host load, we give qualitative and quantitative descriptions of static metrics, such as machine's maximum host load. We also describe dynamic metrics, such as queue state, and relative usage levels compared to capacities. Interestingly, we find that host load in data centers has higher variance than in Grids. This is due to the fact that the Google workload consists of a large number of short tasks that complete within a few minutes.

For the remainder of the paper, we use the terms data center, cluster, and cloud interchangeably. In Section II, we describe our modeling approaches in processing the Google trace data. In Section III, we comprehensively analyze the characteristics of the workload for Google jobs and tasks as compared to those of Grid jobs. We characterize the host load of Google data center in Section IV, by taking into account the variance of priorities, resources, and time periods. Finally, we discuss related work in Section V and conclude with future work in Section VI.

II. SYSTEM MODEL BASED ON GOOGLE'S TRACE DATA

In this section, we first briefly describe Google cluster's task scheduling mechanism, and data contained in the trace.

A Google cluster consists of many machines that are connected via a high-speed network. One or more schedulers receive and process a large number of user requests (a.k.a. jobs), each of which is comprised of one or more tasks. For instance, a map-reduce program will be handled as a job with multiple reducer tasks and mapper tasks. Each task (actually represented as a Linux program possibly consisting of multiple processes) is always generated with a set of usercustomized requirements (such as the minimum CPU rate and memory size). Each job corresponds to one user, and it may contain one or more tasks, which are the smallest units of resource consumption. In addition, different tasks have different scheduling priorities, and there are currently 12 priorities in total.

According to Google's usage trace format [START_REF]Google cluster-usage traces[END_REF], each task can only be one of the following four states, unsubmitted, pending, running and dead, as shown in Figure 1 (excerpted from [START_REF]Google cluster-usage traces[END_REF]). Any newly submitted task will be first put in the pending queue, waiting for the resource allocation (the step (1) as shown in the figure). As long as some available resources meet the task's constraints, it will be scheduled onto the qualified resource for execution (step (2)). Users are allowed to tune their tasks' constraints at runtime, adapting to users' needs and the environment (step (3)). For example, a user can expand the upper bound of CPU rate and memory size for the task in the course of its execution. Any task may be evicted by system, killed or lost during its execution. After its execution, it will enter into the dead state (step (4) and ( 5)), and it can be resubmitted by users later (step ( 6)) to wait for further resource allocation. Basically, the Google scheduler processes high-priority tasks before low-priority ones, and performs the first-comefirst-serve (FCFS) policy on the ones with the same priority. The "best" resources will be used first, in order to optimally balance the resource demands across machines and minimize peak demands within the power distribution infrastructure, reaching approximately optimal resource utilization. Such a model will use the resources according to their abilities, leading to an approximate load balancing situation.

Based on the above task processing model, Google traced over 25 million tasks that were scheduled/executed across over 12500 heterogeneous machines within one month. All the tasks are submitted with a set of customized constraints and priorities, which will be discussed later. More than 10 metrics are collected during the one month of taskevent monitoring, including CPU usage, assigned memory, observed real memory usage, page-cache memory usage, disk I/O time, and disk space.

When releasing the trace, Google normalized almost all floating-point values by its theoretical maximum value. These values were transformed in a linear manner. So, relative information about host load, for example, is preserved.

III. ANALYSIS OF GOOGLE WORKLOAD

In this section, we characterize the workload of jobs on a Google data center, by comparing it to that on other Grid or HPC systems. Based on our characterization, we find that Google jobs behave quite differently from the traditional Grid jobs, especially with respect to job length, frequency of submission, and resource utilization.

Based on Google's trace, we first group the all 25 million tasks in terms of their job IDs, and then compute the statistics (such as mean CPU usage and memory usage) for each job. The Grid/HPC jobs to be compared are also from real-world trace data, available from the Grid Workload Archive (GWA) [START_REF]Grid workloads archive[END_REF] and Parallel Workload Archive (PWA) [START_REF]Parallel workload archive[END_REF] respectively.

1) Job/Task Priority: In a Google data center, each task is submitted with a particular priority, which is selected from 12 levels. Any tasks with high priorities are able to preempt other tasks with lower priorities. All of the tasks that belong to the same job have the same priority. According to the histogram shown in Figure 2, there are three clusters of priorities, low priority (1∼4), middle priority(5∼8), and high priority (9∼12), for both jobs and tasks. Also, we can observe that most of the jobs/tasks are assigned with relatively low priorities (from 1 to 5). So we study later the usage load level or idleness state of the system from the perspective of different task priorities. For example, if a machine's resource utilization is very full but over 90% of execution time is attributed to the tasks with low priorities, the machine can still be considered quite idle w.r.t the tasks that have relatively high priorities (e.g., the values are greater than 4).

2) Job/Task Length: The length of a Google job or a Grid job is defined as the duration between its submission time and its completion time. Its value could be affected by many factors, such as the work to be processed, the priority of the job, the system state (idle or busy) when it is submitted, and so on. We compare the cumulative distribution function (CDF) of the job length between Google and other Grid systems in Figure 3. We observe that Google jobs are quite shorter than Grid jobs: over 80% Google jobs' lengths are shorter than 1000 seconds, while most of Grid jobs are longer than 2000 seconds. This difference is mainly due to the fact that Grid jobs are usually based on complex scientific problems, while Google jobs, such as keyword search, are often real-time.

We also compare Google task lengths to AuverGrid's job lengths. We compare the mass-count disparity of the two systems in Figure 4. (We use the term task length and task execution time interchangeably.) Mass-count [START_REF] Feitelson | Workload Modeling for Computer Systems Performance Evaluation[END_REF] is a very important metric used to summarize the key features (such as heavy tails) for specific distributions. It is made up of the "count" distribution and the "mass" distribution. The "count" distribution simply refers to the cumulative distribution function (CDF) as it counts how many items are smaller than certain size. The "mass" distribution weights each item, specifying the probability that a unit of mass belongs to an item. Specifically, their values are calculated based on Formula (1) and ( 2), where f (t) is referred to as the probability dense function.

F c (x) = Pr(X < x) (1)
F m (x) = x 0 t • f (t) dt ∞ 0 t • f (t) dt (2)
By comparing the two curves, we can determine whether the distribution follows Pareto principle [START_REF] Koch | The 80/20 principle: the secret of achieving more with less[END_REF], heavy tails, or other statistical features. In the analysis, joint ratio (a kind of Gini coefficient [START_REF] Feitelson | Workload Modeling for Computer Systems Performance Evaluation[END_REF]) is a critical measure index, defined as X/Y, meaning that X% of the items account for Y% of the mass and Y% of the items account for X% of the mass. The mm-distance (abbreviated as mmdis.) shown in the figure is defined as the horizontal distance of the two points that are right in the middle of the CDF of the Count curve and Mass curve. Longer distance means a larger number of period lengths each with long duration.

Task lengths are 1.29 times longer on average in Auver-Grid than in Google. However, AuverGrid's maximum task length is 1.61 times smaller than that in the Google cluster. Statistics shows that the average and maximum values of task's execution time among AuverGrid's 340,000 submitted tasks are 7.2 hours and 18 days respectively, while those in Google cluster are 5.6 hours and 29 days respectively. Through our analysis of task length's mass-count disparity (see Figure 4), the distribution of Google's task lengths exhibits the Pareto principle much more than that of Grid's. Specifically, about 94% of tasks' execution times in Google's data center are less than 3 hours. In contrast, only 70% of tasks in AuverGrid are smaller than 12 hours. Such a difference in task length significantly impacts the fluctuation of host load, which we discuss later.

3) Job Submission Frequency: Job submission frequency is evaluated via the submission interval length between two consecutive job submissions. This reflects the interactivity between users/administrators and the systems. Figure 5 presents the CDF of the submission interval length. We can clearly observe that the submission interval length of Google jobs is much shorter than that of Grid jobs, which means that the frequency of Google job submission is much higher than that of Grid jobs. This can also be confirmed by Table I, which shows the minimum/mean/maximum number of jobs submitted per hour in Google and Grid/HPC systems. It is observed that the frequency of Google's job submission is much higher than that of other Grid systems.

We use the fairness index [START_REF] Jain | The Art of Computer Systems Performance Analysis: Techniques for Experimental Design, Measurement, Simulation and Modelling[END_REF] to demonstrate the stability of the job submission frequency. The fairness index is defined in Formula (3), where x i refers to the number of job submissions within one hour in different periods. The higher value of the fairness index, the higher stability of the job submissions. Table I shows that Google jobs are submitted with higher and more stable frequencies than Grid jobs, while Grid job submissions exhibit significantly low fairness because of their strong diurnal periodicity.

f (x) = ( n i=1 x i ) 2 n n i=1 x 2 i (3)

4) Job Resource Utilization:

The last insight about the jobs is on their resource utilization, such as the CPU usage and the memory usage per job. We evaluate the CPU usage by the ratio of the cumulative execution time on one or more processors and its wall-clock time, as shown in Formula (4). For comparison, we use the mean memory size used by the job to evaluate its memory utilization.

CP U U sage = # of CP U • ExeT imeP erCP U W allClockT ime (4)
We present the CPU usage and memory usage in Figure 6 (a) and (b) respectively. In Figure 6 (a), we can observe that the CPU used by Google jobs is always smaller than that of other Grid systems. Specifically, a large majority of Google jobs just need one processor per job at any time, though each job may contain multiple tasks submitted in a sequential order. In contrast, the jobs in AuverGrid and DAS-2 are likely parallel programs simultineously running on multiple execution nodes.

As for memory usage, Google trace does not expose the exact memory size used by jobs but their scaled values compared to the maximum memory capacity of each node. By assuming the maximum memory capacity to be equal to 32GB and 64GB respectively, we could also compare Google jobs and Grid jobs with respect to their memory utilization. In Figure 6 (b), it is observed that the memory size used by Google jobs is always relatively small as users' short-term, interactive jobs are dominant and not resource intensive.

IV. ANALYSIS OF GOOGLE'S HOST LOAD

In this section, we first focus on the statistics of Google clusters' static metrics, such as machine's maximum host load on different attributes. After that, we will analyze two types of dynamic metrics (queuing state and usage levels), which change over time. Our analysis is based on the observation of the Googles' thousands of machines.

A. Analysis of Static Metrics

In the Google cluster, different machines have different capacities for various resource types (CPU, memory, page cache). So it is necessary to characterize the heterogeneity as well as the maximum loads on different machines. In the trace, the normalized hosts' capacities are provided by dividing by the maximum value w.r.t. different attributes. According to statistics, we further infer the upper bound using the maximum resource usage value over the lifetime of the trace for each host. That is, from the task usage on a given machine, we calculate the machine resource usage, and use the maximum resource usage to indicate the machine's valid ability. According to our Google co-author, this is a valuable estimate, since the capacity of usable resource in user space is often less than the full capacity due to system overheads (from the kernel, for instance.)

According to the statistics on 12500 machines, Figure 7 shows the probability distribution of the maximum consumed host load during the whole lifetime of the trace, w.r.t. the four significant attributes, CPU usage, memory consumed, memory assigned, and page cache used. CPU usage is measured based on CPU core seconds per second, that is, the more cores a machine owns, the higher computation ability it would have. The consumed memory indicates the practical value of the memory size consumed by applications, while the assigned memory means the memory size allocated to the applications. The page cache metric records the total amount of Linux page cache (i.e., file-backed memory) including both mapped and unmapped pages. Recall that all the usage values in the Google trace are normalized based on the corresponding maximum node capacity, we can only show the relative heterogeneity instead of the exact amounts. The black dotted line in the figure indicates the normalized heterogeneous capacities, e.g., for all machines, CPU's capacities are 0.25, 0.5 and 1; memory's are 0.25, 0.5, 0.75 and 1. Figure 7 (a) shows that most of machines' maximum host loads are very close to their capacities. Specifically, for the machines with relatively low and middle CPU capacities, and over 80% and 70% of these hosts' maximum loads are equal to their capacities respectively. From Figure 7 (b) and (c), we clearly see that the memory capacity among all machines can be split into four groups, and the expected values of the maximum memory size consumed on machines are kept around 80% of their corresponding capacities, which implies that the memory overhead of each machine must be non-negligable. In comparison, the summed assigned memory size is around 90% of the capacity with high probability. The page cache capacities on machines are the same as each other, while the maximum consumed values shows a clear bimodal distribution.

B. Analysis of Dynamic Metrics

We investigate the characteristics of the dynamic features of the resource metrics appearing in the Google cluster trace data. We also compare the load changes between Grid platform and Cloud platforms, with respect CPU usage and memory usage.

1) Machine's Queuing State: A machine's queuing states refer to the number of tasks that are kept in different states (pending, running, finish, or abnormal) respectively. The distribution of states changes over time as new tasks are submitted/scheduled or old tasks are finished normally/abnormally. Task's abnormal termination is due to one of the following four task events, evicted by higher priority tasks, failed because of task failure, killed by its user, and lost because of its missing source data.

Different task events occur with different frequencies, but different machines exhibit similar distributions on the different types of events. We present the task event trace of a particular machine, in Figure 8 (a). We use the black lines (between start time and end time) to indicate the running period for each task. Based on the different states the tasks are treated over time, we can imagine them being kept in different queues (pending-queue, running-queue, dead-queue) on the machines. For example, we define the running queue state to be the number of running tasks at some time point. Figure 8 (b) also presents the queuing states over time. Specifically, the running-queue state increases gradually from zero to of the trace. The pending-queue state (i.e., the number of pending tasks) is always 0 (except for the bootstrap period of the system), which means that each task is able to be immediately scheduled as it is submitted. Moreover, the number of finished tasks increases linearly and many of them belong to the abnormal-completion state. According to statistics, for the totally 44 million task-completion events , about 59.2% are abnormal ones, among which most of them belong to the fail state (50%) or the kill state (30.7%).

We show more statistics about running tasks on nodes, since the system's dynamic state is mainly attributed to the running tasks instead of the pending or finished ones. We split the running-queue state into 6 intervals based on different number of running tasks, [0,9], [START_REF]Parallel workload archive[END_REF][START_REF] Kondo | Characterizing resource availability in enterprise desktop grids[END_REF], [START_REF] Calder | The entropia virtual machine for desktop grids[END_REF]29], [30,39], [40,49], [50,• • •]. Then, we analyze the mass-count disparity of the period lengths in which the running-queue state is unchanged. We observe the distribution of these lengths over all machines (see Figure 9) exhibits the Pareto principle especially from the perspective of the joint ratio and mmdistance metrics. That is, the majority of the unchangedstate lengths are quite short. In Figure 9, we just show four major intervals for simplicity. It is obvious that the first three intervals follows the 10/90 rule, while the last one ([40,49]) follows the 15/85 rule. The distribution is thus skewed as about 90% of tasks have very short continuous durations and contain 10% mass. Moreover, the first three intervals also show a similar mm-distance, while the last one presents a much smaller mm-distance. This implies that the [40,49] running-queue state changes much more frequently than the other ones. By studying the number of running tasks, we realize that the duration in which the running-queue state is constant on each machine can be several hours on average. This means that the running-queue state is relatively stable.

The relative usage of a machine's resource is most relevant for load prediction. We define a machine's resource usage to be the ratio of the usage of some attribute to its capacity. The range of resource usage is thus in [0, 1]. We will also use resource utilization, usage level, and load level interchangeably. We will characterize in detail CPU and memory resource usage in the following sections.

2) Machine's Usage Level: Based on the scaled CPU capacity information provided by the Google trace, we can compute the relative usage level (or load level) for the CPU and memory relative to their capacity on each machine. As follows, we present a snapshot of the relative CPU/memory usage on 100 randomly sampled machines, before discussing aggregate statistics and distributions.

In terms of the distribution of the number of tasks with different priorities (as shown in Figure 2), we split all jobs/tasks into three categories, with low-priority (1-4), middle-priority (5-8) and high-priority (9-12) respectively. For a particular attribute with a specific priority category (e.g., CPU usage of the mid+high priority tasks), we further divide the whole usage range into five equal intervals, and plot a trace of dynamically changing load levels over time for 50 randomly sampled machines. This depicts the system's idleness level from the perspective of the tasks with different priorities.

From Figure 10 (a)∼(d), we observe that CPU and memory usage differ. Specifically, Figure 10 (a) shows that most of the machines are relatively idle compared to their capacities in most of time (from the beginning to the 21th day and from 25th day till the end). Such an idle situation is attributed to the intention by reserving a certain of CPU resources to meet service level objectives (for instance, a threshold on web request latency) in case of unexpected load spikes. From the perspective of high-priority tasks, however, the percentage of host load is not that heavy, as shown in Figure 10 (b). This means that most of the CPU resources from the Google computer cluster are actually consumed by low-priority tasks, especially in the busy duration (from 21th day to 25th day). Hence, it is relatively easy to allocate idle CPU resources for the tasks with high priorities.

Figure 10 (c) and (d) present the snapshot of the load changes about memory usage. We can see that the majority of machines' memory usage is high compared to their capacities. Moreover, their fluctuations are quite unlike that of CPU usage. Specifically, some machines' memory usage are always relatively lightly-loaded (the green lines shown in the figure); some are often heavily loaded (the black and red lines); some memory usage always alternate between two load levels; and there also exist a few machines that demonstrate completely irregular changing memory usage.

We use Table II and Table III to show the statistics the durations where CPU or memory usage are constant. The statistics are computed using all the tasks (including all priorities) running on all of machines. It is obvious that the Moreover, the frequency of the CPU load changes is very high, in that the average duration is only about 6 minutes and the mm-distance is also small (18∼49 minutes).

Due to the limited space, we cannot show the complete statistics from the perspective of middle/high-priority tasks, but just present a few as follows: As for CPU usage, the uninterrupted load duration is 7∼8 minutes on average, and the joint ratio is always about 30/70, in the situation where only middle/high-priority tasks or high-priority tasks are taken into account; as for memory usage, the joint ratios are about 20/80 and 15/85 for the middle/high-priority tasks and high-priority tasks respectively. In particular for high-priority tasks, CPU's mm-distance mainly appears in 14∼40, except for [0,0.2] whose mm-distance is 371, while memory's is in [560,1700], which clearly indicates that the CPU usage changes much more frequently than that of memory usage over time. Now, we focus on the distribution of the mass-count disparity of the relative resource utilization whose values range from 0% to 100%. Figure 11 and Figure 12 present such statistics about CPU usage and memory usage respectively. By comparing the two figures, we can further confirm that the CPU usage is much lower than memory usage, relatively. Specifically, the percentage load of CPU is about 35% w.r.t. all the tasks and about 20% for the high-priority tasks, while memory's are about 60% and 50% respectively. In addition, we can also know that the distribution of usage is relatively uniform because of the small value of mm-distance and large value of joint ratio. Finally, we explore the difference of load changes between Cloud and Grid, by comparing the resource usage of Google's compute clusters and those of AuverGrid [START_REF] Auvergrid | [END_REF] and SHARCNET [START_REF]Sharcnet grid project[END_REF], whose trace data can be downloaded from the Grid Workload Archive (GWA) [START_REF]Grid workloads archive[END_REF].

In Figure 13, we present the CPU usage load of three machines, selected from the Google cluster, AuverGrid and SHARCNET, within one month, 5 days and 1 day respectively. What is the most interesting is that Grid's CPU usage is always higher than memory usage. This because most Grid jobs are computation-intensive. By contrast, in the Google cluster, a machine's CPU usage is usually lower than memory usage. This implies that Google tasks are not compute-intensive programs and they use other resources such as memory more intensively.

The second key difference observed is that Google cluster's CPU load has higher noise than the other two Grid systems. We measure noise by processing the trace with a mean filter [START_REF] Diniz | Adaptive Filtering: Algorithms and Practical Implementation[END_REF], and then computing statistics on the transformed trace. The minimum/mean/maximum noise of AuverGrid's CPU load are 0.00008, 0.0011, 0.0026 respec-tively, while those for the Google cluster's CPU load are 0.00024, 0.028, 0.081 respectively. It is clear that the noise of Google cluster's usage load is about 20 times as large as that of Grid's on average. (i) sharcnet [START_REF]Parallel workload archive[END_REF][START_REF] Feitelson | Workload Modeling for Computer Systems Performance Evaluation[END_REF] Figure 13. Host load Comparison between Google cluster & Grid Systems Lastly, Figure 13 also shows that the host load of Au-verGrid and HARCNET is more stable during relatively long periods (e.g., a few hours), while Google cluster's load changes rather frequently even in quite short periods, tens of minutes in length. Over all hosts, the mean autocorrelation of CPU usage load is about -8 × 10 -6 , which is much smaller than AuverGrid's CPU autocorrelation (=1.006×10 -8 ). This is mainly due to the fact that Google cluster has to process millions of tasks, a large majority of which are smaller than in Grids, as shown in Figure 8 (a) and Figure 4. The fact that majority of Google's tasks are smaller will definitely lead to much finer resource usage and introduce more instability in terms of host load fluctuation. Moreover, as described in [START_REF] Sharma | Modeling and synthesizing task placement constraints in google compute clusters[END_REF], Cloud tasks' placement constraints may also be tuned by users frequently over time, which may further impact the resource utilization significantly. All in all, it is more challenging to predict Google cluster's host load because of its higher noise and and more unstable state.

V. RELATED WORK

There is much related work on the characterization of workload and host load for Grid systems. H. Li [START_REF] Li | Workload dynamics on clusters and grids[END_REF] studied the workload dynamics on clusters and Grids, based on the real-world trace provided by Grid Workload Archive (GWA) [START_REF]Grid workloads archive[END_REF]. His work revealed many features about Grid load, which can be leveraged to improve the system performance. Specifically, the Grid host load exhibits clear periodic or diurnal patterns, which can be used for host load prediction for the task scheduling [START_REF] Li | Analysis and synchesis of pseudo-periodic job arrivals in grids: A matching pursuit approach[END_REF]. E. Afgan et al. [START_REF] Afgan | Exploiting performance characterization of BLAST in the grid[END_REF] exploited the performance characteristics of the NIBC BLAST application on a cluster environment. In our research, we comprehensively compare the workload between Cloud (Google) and other Grid systems, showing much shorter job length and much higher submission frequency in Cloud systems. D. Kondo [START_REF] Kondo | Characterizing resource availability in enterprise desktop grids[END_REF] analyzed the resource availability of Desktop Grid systems based on the traces collected from different Grid softwares, including the Entropia DCGrid desktop grid [START_REF] Calder | The entropia virtual machine for desktop grids[END_REF], XtremWeb desktop grid [START_REF] Germain | Xtremweb: Building an experimental platform for global computing[END_REF] and so on. Our key contribution is to exploit the key insights about the workload and host load between Cloud and Grid systems, serving as a fundamental basis for the system performance improvement over Cloud systems.

Despite some works [START_REF] Sharma | Modeling and synthesizing task placement constraints in google compute clusters[END_REF], [START_REF] Zhang | Characterizing task usage shapes in google compute clusters[END_REF], [START_REF] Barnes | A regression-based approach to scalability prediction[END_REF] also characterizing workload features for Cloud systems, they mainly focus on modeling running tasks such as tasks' placement constraints [START_REF] Sharma | Modeling and synthesizing task placement constraints in google compute clusters[END_REF] or usage shapes [START_REF] Zhang | Characterizing task usage shapes in google compute clusters[END_REF]. Specifically, B. Sharma et al. [START_REF] Sharma | Modeling and synthesizing task placement constraints in google compute clusters[END_REF] carefully studied the performance impact of task placement constraints based on the resource utilization from the view of tasks. Q. Zhang et al. [START_REF] Zhang | Characterizing task usage shapes in google compute clusters[END_REF] designed a model that can characterize task usage shapes in Google's compute clusters. In comparison to these two works, we intensively compare the load characteristics between Cloud and Grid systems, and summarize many new insights about Cloud load.

A. Khan et al. [START_REF] Khan | Workload characterization and prediction in the cloud: A multiple time series approach[END_REF] designed a model to capture the CPU load auto-correlations in the Cloud data centers by leveraging the Hidden Markov Model (HMM). B.J. Barnes et al. [START_REF] Barnes | A regression-based approach to scalability prediction[END_REF] introduced a regression-based approach for predicting parallel application's workload, and the prediction errors can be limited between 6.2% and 17.3%. In comparison, by using Google cluster's large-scale one-month trace data, we comprehensively studied the host load changes about multiple resource attributes (including CPU, memory and page cache) in Cloud data centers. Specifically, we characterize the machine capacity, queuing size (the number of running tasks), the impact of task priority and task events, usage levels and so on.

VI. CONCLUSION AND FUTURE WORK

We characterize the work load and host load in a Google data center and compare them to those of Grid systems. We use a detailed workload trace of a Google data center with over 25 million tasks executed over 12,000 hosts in a 1month time frame, and the trace data from many Grid/HPC systems. Our key findings are summarized below:

• With respect to Work load of Job or Task:

-Job/Task length: Most jobs finish in tens of minutes and tasks finish quickly on the order of a few minutes. Statistics show that about 55% of tasks finish within 10 minutes and about 90% of tasks' lengths are shorter than 1 hour. A handful of tasks last for several days or weeks and likely correspond to long-running services. Compared to Grid workloads, most Cloud task lengths are shorter, and at the same time the longest task lengths are longer. We believe this difference is due to the differences in Cloud users and applications, which include commercial applications such as web services. This indicates that load can be sporadic and highly variable. -Job Priority: We find that task priorities can be clustered into 3 groups, namely, low, medium, and high. As high priority task can preempt low priority task, load prediction thus can and should be tailored and evaluated for each of these groups. -Job Submission Frequency: Google jobs are submitted with much higher and more stable frequency than that of Grid jobs. The average number of jobs submitted per hour and its fairness index in Google are 552 and 0.94 respectively, compared to 8.4∼126 and 0.04∼0.51 in Grids. -Job Resource Utilization: Google jobs usually have lower resource demand for CPU and memory than Grid jobs, because Google jobs (such as keyword search) are more interactive and real-time in constrast to scientific batch jobs.

• With respect to Host load: -Maximum load: We find that Google host's maximum CPU load is often close to the CPU capacity, and the maximum memory usage is about 80% of the memory capacity. The maximum load is actually controlled in the Google system for guaranteeing the service level of requests in case of unexpected load spikes. In contrast, Grid resources can be highly utilized without having a high risk of losing users or customers. -Machine usage level: CPU and memory usage changes every 6 minutes, indicating again the volatility of load. CPUs are often idle, but memory usage is relatively high. CPU usage in Grids is higher and more stable. Noise of CPU load in the Google cluster is 20 times as high as that in Grids. In the future, we will try to exploit the best-fit load prediction method based on our characterization work, and analyze and improve the job scheduling in Google data centers.

Figure 1 .

 1 Figure 1. State Transitions of Google Tasks

Figure 2 .

 2 Figure 2. Statistics based on Different Priorities

Figure 3 .

 3 Figure 3. The CDF of Job Length of Google and Grid Systems

Figure 4 .

 4 Figure 4. Mass-Count Disparity of Task Lengths (Google v.s AuverGrid)

Figure 5 .

 5 Figure 5. The CDF of Submission Interval of Google and Grid Systems

Figure 6 .

 6 Figure 6. CPU & Memory Usage of Jobs

Figure 7 .

 7 Figure 7. Distribution of Maximum Load

Figure 8 .

 8 Figure 8. Task Events and Queuing State on a Particular Host

Figure 9 .

 9 Figure 9. Mass-Count of Duration in Unchanged Queuing State

Figure 10 .

 10 Figure 10. Snapshot of Resource Usage Load

Table I THE

 I NUMBER OF JOBS SUBMITTED PER HOUR

		Google AG	NG	SN	ANL	RICC MT	LLNL
	max. #	1421	818	2175	22334	132	4919	2315	240
	avg. #	552	45	27	126	10	121	24	8.4
	min. #	36	0	0	0	0	0	0	0
	fairness	0.94	0.35	0.11	0.04	0.51	0.14	0.04	0.23

[0,0.2] [0.2,0.4] [0.4,0.6] [0.6,0.8] [0.8,1] avg value (minute)

			Table II
	CONTINUOUS DURATION OF UNCHANGED CPU USAGE LEVEL
		6	6	6	6	5
	max value (minute)	41269	40302	3590	1452	2575
	joint ratio	26/74	28/72	30/70	30/70	27/73
	mm-distance (minute)	49	25	18	19	24
			Table III
	CONTINUOUS DURATION OF UNCHANGED MEMORY USAGE LEVEL
		[0,0.2]	[0.2,0.4]	[0.4,0.6]	[0.6,0.8]	[0.8,1]
	avg value (minute)	6	9	10	10	10
	max value (minute)	31556	40302	10996	16736	16826
	joint ratio	20/80	23/77	26/74	23/77	18/82
	mm-distance (minute)	119	83	63	95	351

ACKNOWLEDGMENTS

We thank Google Inc, in particular Charles Reiss and John Wilkes, for making their invaluable trace data available. This work was made possible by a Google Research Award and by the ANR project Clouds@home (ANR-09-JCJC-0056-01).