Ben Ismail

Dhouha Kbaier Ben

Lazure Pascal

Puillat Ingrid

Dhouha Kbaier
email: dhouha.kbaier@ifremer.fr

Statistical properties and time-frequency analysis of temperature, salinity and turbidity measured by the MAREL Carnot station in the coastal waters of Boulogne-sur-Mer (France)

Keywords: Continuous wavelet transform, cross-correlation, empirical mode decomposition, Hilbert-Huang transform, stationarity, time dependent intrinsic correlation, time series, wavelets Section

In marine sciences, many fields display high variability over a large range of spatial and temporal scales, from seconds to thus ousands of years. The longer and longer recorded time series, with an increasing sampling frequency, in this field are often nonlinear, nonstationary, multiscale and noisy. Their analysis faces new challenges and thus requires the implementation of adequate and specific methods. The objective of this paper is to highlight time series analysis methods already applied in econometrics, signal processing, health,etc. to the environmental marine domain, assess advantages and inconvenients and compare classical techniques with more recent ones. Temperature, turbidity and salinity are important quantities for ecosystem studies. The authors here consider the fluctuations of sea level, salinity, turbidity and temperature recorded from the MAREL Carnot system of Boulogne-sur-Mer (France), which is a moored buoy equipped with physico-chemical measuring devices, working in continuous and autonomous conditions. In order to perform adequate statistical and spectral analyses, it is necessary to know the nature of the considered time series. For this purpose, the stationarity of the series and the occurrence of unit-root are addressed with the Augmented-Dickey Fuller tests. As an example, the harmonic analysis is not relevant for temperature, turbidity and salinity due to the nonstationary condition, except for the nearly stationary sea level datasets. In order to consider the dominant frequencies associated to the dynamics, the large number of data provided by the sensors should enable the estimation of Fourier spectral analysis. Different powerspectra show a complex variability and reveal an influence of environmental factors such as tides. However, the previous classical spectral analysis, namely the Blackman-Tukey method, requires not only linear and stationary data but also evenly-spaced data. Interpolating the time series introduces numerous artifacts to the data. The Lomb-Scargle algorithm is adapted to unevenly-spaced data and is used as an alternative. The limits of the method are also set out. It was found that beyond 50% of missing measures, few significant frequencies are detected, several seasonalities are no more visible, and even a whole range of high frequency disappears progressively. Furhermore, two time-frequency decomposition methods, namely wavelets and Hilbert-Huang Transformation (HHT), are applied for the analysis of the entire dataset. Using the Continuous Wavelet Transform (CWT), some properties of the time series are determined. Then, the inertial wave and several low-frequency tidal waves are identified by the Please note that this is an author-produced

Introduction

The birth of observatories contributes to improve the environmental knowledge. Indeed, long-term research infrastructures such as NEPTUNE Canada and in Europe: EMSO, FIXO3 and JERICO are providing longer and longer time series with high frequency records in open sea and coastal areas. The analysis of such time series requires the implementation of procedures and methods suitable to this increasing data flow towards an automatisation. Consequently, the objective of this paper is to explore adequate analysis methods, assess advantages and inconvenients and compare classical techniques with more recent ones. Coastal ecosystems, especially those of transitional waters, display high variability over a large range of spatial and temporal scales. To study their dynamics and estimate their variations at all scales, high frequency measurements are needed [START_REF] Dickey | Optical and physical variability on timescales from minutes to the seasonal cycle on the new England shelf: July 1996 to June 1997[END_REF][START_REF] Chavez | Moorings and drifters for real-time interdisciplinary oceanography[END_REF][START_REF] Dickey | The emergence of concurrent high resolution physical and bio-optical measurements in the upper ocean and their applications[END_REF][START_REF] Dickey | An experiment in duplex communication with a multi-variable moored system in coastal waters[END_REF] together with appropriate analysis techniques. We present here the data recorded every twenty minutes by an autonomous monitoring device, settled in coastal waters of Boulogne-sur-Mer (eastern English Channel, France), recorded from 2005 to 2009. The fixed buoy station can record various biogeochemical parameters simultaneously. The time series of sea level, temperature, salinity and turbidity are considered. The fluctuations of the four quantities are important for ecosystem studies. The paper is structured in the seven following sections. After the description of the background and objectives in section 1, section 2 is dedicated to the presentation of the data. Then, stationarity issues are addressed in section 3. A preliminary analysis of stationarity focuses on the graphs and the correlograms. The conclusions from the direct observation of these are reinforced with application of the Augmented Dickey-Fuller (ADF) tests [START_REF] Dickey | Distribution of the estimates for autoregressive time series with a unit root[END_REF]. Section 4 is devoted to a spectral analysis of the time series. First, a harmonic analysis is described. Secondly, the Blackman-Tukey spectral analysis is performed and comparisons between different periodograms are made. The problem of missing data is described and the Lomb-Scargle algorithm [START_REF] Lomb | Least-Squares Frequency Analysis of Unevenly Spaced Data[END_REF][START_REF] Scargle | Studies in Astronomical Time Series Analysis. I. Modeling Random Processes in the Time Domain[END_REF] is used. This alternative does not require any interpolation of the data, in contrast to the Blackman-Tukey method. Furthermore, the limits of the Lomb-Scargle are illustrated. Continuous Wavelet Transform (CWT) and Hilbert-Huang Transformation (HHT) approaches are presented in section 5, while section 6 deals with an Empirical Mode Decomposition (EMD) based cross correlation method: Time Dependent Intrinsic Correlation (TDIC). Finally, a conclusion section (section 7) provides also possible developments of this study.

A C C E P T E D

M A N U S C R I P T Table 1 presents some statistical properties of the time series. The normalized third central moment is called the skewness and represents a measure of asymmetry of the tails of a distribution. The mode is the most common value, or the peak in the histogram. Fig. 4 displays different probability density functions (PDFs) of the data sets. We notice for the salinity that mean < median < mode, which corresponds to a negative skewness of 49 . 2  and indicates that the distribution is spread out more to the left of the mean value [START_REF] Panofsky | Some applications of statistics to meteorology[END_REF]). On the other hand, the distribution of the turbidity with a positive skewness has large tails that extend to the right. The distributions of temperature and sea level are almost symmetric. This observation is consistent with the corresponding values of skewness, since the skewness of the symmetric normal distribution is zero.

In Boulogne-Sur-Mer, the temperature of the sea remains almost the same during August and September as well as during February and March. Whereas the passage from the wintry temperatures to the summer temperatures is very quick and the offseason is very short. The net sheet of the heat due to the ocean-atmosphere exchanges is then a seasonal curve with a sinusoidal shape. Therefore, the histogram of temperature shows a bimodal pattern with two modes around 9°C and 18°C. Also, the sea level histogram shows a bimodal pattern where the two peaks correspond to slack tides, i.e., high tide and low tide around 3 m and around 7 m respectively. The width of each peak is related to the Spring-Neap tidal cycle.

Before estimating the statistical models, it is important to analyze whether the series are stationary or not. Stationarity issues are addressed in the next section.

Section 3: Stationarity tests

The following definition of stationarity is taken from Challis and Kitney, 1991. "Stationarity is defined as a quality of a process in which the statistical parameters of the process do not change with time". The "weak" stationary condition requires that the mean, variance and covariance structure of each variable do not change over time. It is important to analyze whether the series are stationary or not since it is not valid to use standard statistical tests in the negative case [START_REF] Nelson | Trends and random walks in macroeconomics time series: some evidence and implications[END_REF]. Estimating the correlation between nonstationary variables can give spurious dependencies and conclusions about the development of temperature can be erroneous. If the dependent variable has a deterministic trend, a trend-adjusted series could be used instead. In case the variable has a unit-root, thus it has to be differenced. Further details about stationarity and its implications are available in Emery and Thompson (2014).

A preliminary analysis of stationarity is based on the graphs of the time series as well as the correlograms. Correlograms were used in early work such as the seminal paper by [START_REF]Why do we sometimes get nonsense correlations between time series? -A study in sampling and the nature of time series (with discussion)[END_REF], the monograph by [START_REF] Wold | A study in the analysis of stationary time series[END_REF] and the analysis of agricultural price series by [START_REF] Kendall | Oscillatory movements in English agriculture (with discussion)[END_REF]. As shown in Figs. 1-3, the graphs of temperature, turbidity and salinity are characterized by:  A non stationarity in mean showed by a slight linear upward or downward trend. For instance, the red linear fitting plotted in Fig. 1 For the sea level plotted in Fig. 6, the time series is stationary since there is no systematic change in mean (no trend) and variance. Regarding the correlogram, Fig. 7 shows that it has few significant spikes at small lags (lag <5); and decreases drastically at 8 or 9 lags. Therefore, the correlogram confirms that the sea level series is stationary. Nevertheless, the correlogram of the turbidity in Fig. 8 is not declining and shows persisting significant values, indicating that the time series is nonstationary. Differencing is a transformation of the series to a new time series where the values are the differences between consecutive values. This procedure may be applied iteratively more than once, giving rise to the "first differences", "second differences", etc. Differencing makes the time series stationary, as the correlogram has only a couple of significant spikes at lags 1 and 2. Similar results are obtained for salinity and temperature. After a preliminary inspection of the graphs, we move to a quantitative study of the stationarity. Indeed, using non-stationary time series data in environmental models may lead to poor understanding and forecasting. The solution to the problem is to transform the time series data so that it becomes stationary. If the non-stationary process is a random walk with or without a drift, it is transformed to stationary process by differencing. On the other hand, if the time series data analyzed exhibits a deterministic trend, the spurious results can be avoided by detrending. Furthermore, the nonstationary series may combine a stochastic and deterministic trend at the same time and to avoid obtaining misleading results both differencing and detrending should be applied, as differencing will remove the trend in the variance and detrending will remove the deterministic trend. Sometimes it can be hard to tell the difference between a series that is trend-stationary and one that is difference-stationary. Therefore, the ADF unit root tests were used to get a more definitive answer. These are statistical hypothesis tests of stationarity that are designed for determining more objectively whether differencing is required. That is why we have applied Augmented Dickey-Fuller (ADF) tests [START_REF] Dickey | Distribution of the estimates for autoregressive time series with a unit root[END_REF] for testing of unit root and stationarity. Unit root corresponds to the application of the difference operator  to induce stationarity. Given an observed time series

( 1 Y , 2 Y , …, N Y
), Dickey and Fuller consider three differentialform autoregressive equations to detect the presence of a unit root:

1.          p j t j t j t t Y Y Y 1 1    2.           p j t j t j t t Y Y Y 1 1     3.            p j t j t j t t Y Y t Y 1 1     
where t is the time index,  is an intercept constant called a drift,  is the coefficient on a time trend,  is the coefficient presenting process root, i.e. the focus of testing, p is the lag order of the first-differences autoregressive process and t  is an independent identically distributes residual term. The difference between the three equations concerns the presence of the deterministic elements  (a drift term) and



. More detailed, the null and alternative hypotheses corresponding to the models above are:

( 1 h ) 0 H : t Y   0  is random walk 1 H : t Y   0  is stationary process ( 2 h ) 0 H :   t Y    0 , 0   is random walk around a drift 1 H :   t Y    0 , 0   is level stationary process ( 3 h ) 0 H :   t Y    0 , 0   is random walk around a trend 1 H :   t Y    0 , 0  
is trend stationary process An important practical issue for the implementation of the ADF test is the specification of the lag length p . If p is too small then the remaining serial correlation

A C C E P T E D M A N U S C R I P T
in the errors will bias the test. If p is too large then the power of the test will suffer.

The results of the Monte Carlo experiments suggest that it is better to error on the side of including too many lags. [START_REF] Ng | Unit Root Tests in ARMA Models with Data-Dependent Methods for the Selection of the Truncation Lag[END_REF] (1989), is

                4 1 max 100 12 N p
where  

x denotes the integer part of x .

To estimate the significance of the coefficients in focus, the modified t (Student)statistic (known as Dickey-Fuller statistic) is computed and compared with the relevant critical value [START_REF] Mackinnon | Numerical distribution function for unit root and cointegration tests[END_REF]: if the test statistic is less than the critical value then the null hypothesis is rejected. Each version of the test has its own critical value which depends on the size of the sample.

The ADF tests results for salinity are presented in Fig. 9. The test of constant and a deterministic trend occurring results in the computed absolute ADF statistic (3.12), which is smaller than the tabulated critical value (11.62) at 90 % level of confidence. Therefore, the hypothesis that the time series is nonstationary and has a trend is accepted. The hypothesis of unit root with a drift is also accepted since the computed absolute ADF statistic (2.56) is smaller than the tabulated critical value (11.55) at 90 % level of confidence However, the first ADF test without a constant is rejected since the test values are greater (in absolute value) than the critical value (0.36) for different confidence levels.

Similar results are obtained for temperature and turbidity data: all the critical values [START_REF] Mackinnon | Numerical distribution function for unit root and cointegration tests[END_REF] are the same and only the test values are different. Thus, the unit-root tests exclude that the series are pure random walk processes (first model). Nevertheless, these series are random walk processes with a drift (second model). Moreover, they are nonstationary and have a deterministic trend (third model). In the next section, we perform power spectral analysis of these data.

Section 4: Power spectral analysis

Harmonic analysis

Harmonic analysis takes advantage of the fact that one knows a priori all the frequencies at which tidal energy will be found in a data time series. Most of these frequencies are astronomically caused but many are due to the nonlinear hydrodynamic effects of shallow water. Indeed, the shallower the water depth is the more the tidal wavelength will shorten and the faster the tidal characteristics of a waterway will change with horizontal (geographic) distance.

In classical harmonic analysis, the tidal signal is modeled as the sum of a finite set of sinusoids at specific frequencies related to astronomical parameters. The purpose of the tide analysis is to determine the amplitude and phase (the so-called tidal harmonic constants) of the individual cosine waves, each of which represents a tidal constituent identified by its period in mean solar hours. The tidal response is modelled as:

k A k  k B k 
where N constituents are used. Each constituent has a frequency σ k which is known from the development of the potential, and complex amplitudes (A k and B k ) which are not known. A possible offset and (optional) linear drift is handled by the first two terms α and β.

One widely-used program for performing harmonic analysis of oceanic tides is is the set of programs in MATLAB created by [START_REF] Pawlowicz R | Classical tidal harmonic analysis including error estimates in MATLAB using T_TIDE[END_REF] to perform classical harmonic analysis for periods of about 1 year or shorter. They have rewritten the FORTRAN package created by Mike [START_REF] Foreman | Manual for tidal heights analysis and prediction[END_REF], incorporating many useful extras. Confidence intervals for the amplitudes and phases of the tidal constituents are calculated from the estimated variance. Confidence intervals were computed from the average of power spectral density within frequency bands centered on 0; 1; 2; ... cycles/day, ±0.1 cycles/day (M 0 and M 1 ), ±0.2 cycles/day (M 2 to M 5 ), ±0.21 cycles/day (M 6 ), 0.26 -0.29 cycles/day (M 7 ), and 0.3 -0.5 cycles/day (M 8 ) [START_REF] Pawlowicz R | Classical tidal harmonic analysis including error estimates in MATLAB using T_TIDE[END_REF].

We have performed harmonic analysis of the data sets presented in section 2 "Presentation of the experimental database". For the sea level, the harmonic analysis is able to predict 99.7% of the variance of the original signal (see Fig. 9). The highest tidal amplitudes and the corresponding waves can be found in Table 2. Thus, the semidiurnal wave 2 M is the most energetic one. For the salinity, the harmonic analysis is able to predict 18.4% of the variance of the original signal. This is due to freshwater contents (received from the rivers) and a seasonal cycle with higher variability linked to the regime of rivers and residuel (subtidal) circulation. Regarding the turbidity, it contains tidal components, but also a biological signal such as the one induced by the plankton in the spring, high-frequency effects linked to storms and a cycle of tidal streams (suspension and advection) [START_REF] Burchard | The formation of estuarine turbidity maxima due to density effects in the salt wedge. A hydrodynamic process study[END_REF][START_REF] Chant | Particle trapping in a stratified flood-dominated estuary[END_REF][START_REF] Goni | Fluxes and sources of suspended organic matter in an estuarine turbidity maximum region during low discharge conditions[END_REF]. Such components make the time series non stationary. This is the reason why the harmonic analysis is even worse: only 8.3% of the original signal's variance is predicted. Hence, the harmonic analysis is not relevant except for the sea level. An adequate spectral analysis is introduced hereafter and allows comparing different periodograms.

Amplitude

Blackman-Tukey spectral analysis

The most popular method to compute powerspectra in earth sciences is the method introduced by Blackman and [START_REF] Blackman R | The measurement of power spectra from the point of view of communication engineering[END_REF]. This method performs autospectral analysis in three steps: calculation of the autocorrelation sequence, windowing and finally computation of the discrete Fourier transform. For the modified periodogram, several window shapes are available such as Bartlett (triangular), Hamming (cosinusoidal) and Hanning (slightly different cosinusoidal). An enhancement was made by the application of the Welch's method (1967). First, the time series is divided into overlapping segments. After the data is split up, the overlapping segments are windowed, like for the modified periodogram. Then, the Welch method includes computing the powerspectrum for each segment and finally averaging the powerspectra. To be more specific, we divide an point data sequence , , into segments of samples each. The Welch periodogram is an modified or windowed periodogram and can be defined as where , the energy in the window is The dynamics of turbidity is illustrated in Fig. 10 showing two powerspectra in log-log plot: the modified periodogram using a Hamming window and the Welch powerspectrum. Both reveal a very complex dynamics, with a turbulent-like behaviour at high frequencies, with a regime following roughly a power-law. This analysis is coherent with the power spectrum of turbidity obtained by [START_REF] Schmitt | Statistical properties of turbidity, oxygen and pH fluctuations in the Seine river estuary (France)[END_REF]. For low frequencies, deterministic forcings are visible, associated with several energetic spikes such as the diurnal tidal cycle 1 K , the semidiurnal cycle 2 M and the corresponding propagation of

A C C E P T E D M A N U S C R I P T
2 4 2M M 
. The diurnal period 1 K (23,93 h) is very close to the day duration and the corresponding spike may hide biological effects on turbidity like the nychtemeral vertical migration of zooplankton. Compared to the modified periodogram, the Welch method improves the signal-to-noise ratio of the spectrum versus a loss of resolution. For example, there are not any spikes in the Welch powerspectrum for the monthly ( Mm ) and semi monthly ( Mf ) lunar tides due to an important loss of resolution at these low frequencies. Even more, the annual cycle ( Sa ) is not detected at all. Indeed, for the Welch power spectrum, the raw data has been divided into several segment, the lowest frequency that could be resolved thus depends on the length of the segment.

In earth sciences, time series are often unevenly spaced. Missing values represent a common characteristic property of autonomous monitoring data bases due to routine maintenance, inaccessibility, vandalism, removal of biofouling, and failure of the measuring devices. The Blackman-Tukey method, however, requires evenly-spaced data. Therefore, we have interpolated the turbidity time series with 13% of missing data in order to generate the powespectra in Fig. 10, as done in some studies [START_REF] Ibanez | Prediction of missing values and detection of 'exceptional events' in a chronological planktonic series: a single algorithm[END_REF][START_REF] Paparella | Filling gaps in chaotic time series[END_REF]. Nevertheless, interpolation introduces numerous artifacts to the data, both in the time and the frequency domain. For this reason, an alternative method of time-series analysis has become increasingly popular in earth sciences, the Lomb-Scargle algorithm [START_REF] Lomb | Least-Squares Frequency Analysis of Unevenly Spaced Data[END_REF][START_REF] Scargle | Studies in Astronomical Time Series Analysis. I. Modeling Random Processes in the Time Domain[END_REF][START_REF] Scargle | Studies in Astronomical Time Series Analysis. II. Statistical Aspects of Spectral Analysis of Unevenly Spaced Data[END_REF][START_REF] Scargle | Studies in Astronomical Time Series Analysis. III. Fourier Transforms, Autocorrelation Functions, and Cross-Correlation Functions of Unevenly Spaced Data[END_REF][START_REF] Press | Numerical Recipes in Fortran 77[END_REF][START_REF] Schulz | SPECTRUM: Spectral Analysis of Unevenly Spaced Paleoclimatic Time Series[END_REF].

Lomb-Scargle powerspectrum

The Lomb-Scargle approach was done in the field of astronomy, corresponding to a framework which can be quite different from time series in coastal oceanography. In contrast to the Blackman-Tukey method, the Lomb-Scargle algorithm evaluates the data of the time series only at times that are actually measured. Detailed discussions of the Lomb-Scargle method are given in [START_REF] Scargle | Studies in Astronomical Time Series Analysis. III. Fourier Transforms, Autocorrelation Functions, and Cross-Correlation Functions of Unevenly Spaced Data[END_REF] and [START_REF] Press | Numerical Recipes in Fortran 77[END_REF]. An excellent summary of the method and a TURBO PASCAL program to compute the normalized Lomb-Scargle powerspectrum of paleoclimatic data has been published by [START_REF] Schulz | SPECTRUM: Spectral Analysis of Unevenly Spaced Paleoclimatic Time Series[END_REF]. A comfortable MATLAB algorithm lombscargle to compute the Lomb-Scargle periodogram has been published by Brett Shoelson (The MathWorks Inc.).

We suppose that a physical variable is measured at a set of time . Then, the Lomb-Scargle peridogram as a function of angular frequency is given by where is defined by

In the next, we consider the salinity time series over the year 2005, where 18.76% of the values are missing. Instead of filling the discontinuities with interpolated data, the spectral analysis is performed using the Lomb-Scargle algorithm (see Fig. 11). Similarly to the turbidity, a turbulent-like stochastic behaviour is visible at high frequencies, whereas the dynamics at low frequencies are dominated by deterministic forcings. Several energetic spikes are observed for the original time series of salinity. The first period corresponds roughly to the year ( Sa ), then the monthly ( Mm ) and semi monthly ( Msf ) lunar tides. Furthermore, the diurnal and semidiurnal cycles are identified as for the previous analysis of turbidity. The peak between Mm and Msf seems to be more pronounced (

2 10 5    f days).
It corresponds to a period of 20 days which is longer than an atmospheric disturbance. Thus, it could be related to the river discharges.

A C C E P T E D

M A N U S C R I P T The method is adapted for unevenly spaced data, however, one might ask the following questions. What are the limits? If the time series has a high rate of missing values, would the Lomb-Scargle powerspectrum be significant? Would the method always detect the most energetic frequencies? Or would it suffer above a certain threshold? To answer these questions, we have removed intentionally and randomly other data from the original time series as follows:

1. Case 1: Remove randomly one month of measures, representing 8.22% of additional missing values. Thus, the rate of missing values is increased from 18.76% to nearly 27%. 2. Case 2: Remove randomly 3 months of measures, representing 24.65% of additional missing values. Thus, the rate of missing values is increased from 18.76% to 43.41%. 3. Case 3: Remove randomly 6 months of measures, representing 49.31% of additional missing values. Thus, the rate of missing values is increased from 18.76% to nearly 68%. 4. Case 4: Remove randomly 90% of the data. 3 provides the number of significant frequencies above the 95% confidence level, dropping down from 131 for the original data of salinity to only 14 when 90% of the data are randomly eliminated. Several seasonalities are no more visible in the spectra and even a whole range of high frequency disappears progressively. Beyond 3 months of discarded data, the Lomb-Scargle algorithm really suffers from a high rate that approaches 50% of missing measures. Although the Fourier transform is a traditional spectral analysis method, it is predicated on a priori selection of basis functions. The Fourier transform is suitable for analyzing periodic signals like ocean tides or Earth tides. The signal is compared with harmonic oscillations of different frequencies. The Fourier Transform gives information about the frequency content of the signal without any time localization. Thus, time-varying amplitudes and/or periods can hardly be detected. Consequently, the Fourier analysis is also constrained by the assumptions of stationarity and linearity. To accommodate the variety of data generated by nonlinear and nonstationary processes, [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF][START_REF] Huang | A new view of nonlinear water waves: the Hilbert spectrum[END_REF] developed a new adaptive data analysis method designated by the NASA as the Hilbert-Huang Transform (HHT) and introduced hereafter. Unlike the Fourier transform, which transforms information from the time domain into the frequency domain, the HHT does not lose temporal information after transformation, i.e. energy-frequency information is maintained in the time domain.

Section 5: Hilbert Huang Transform

While the natural physical processes are mostly nonlinear and nonstationary, the available data analysis methods are either for linear but nonstationary processes such as wavelet analysis and the Wagner-Ville distribution, or nonlinear but stationary processes (see, for example, Tong 1990; [START_REF] Kantz | Nonlinear Time Series Analysis[END_REF][START_REF] Diks | Nonlinear Time Series Analysis: Methods and Applications[END_REF]. The development of the HHT was motivated by the need to describe nonlinear distorted waves in detail, along with the variations of these signals that naturally occur in nonstationary processes [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]. The HHT consists of the following steps:

 sifting, that is, empirical adoption of the Principal Component Analysis (PCA) for multi-components in the signal to rearrange the signal in terms of local bases that are nearly orthogonal each other;  physically based construction of instantaneous frequency whose concept is applicable to nonstationary and nonlinear signals;  complexification of the signal via the Hilbert Transform to characterize the signal in terms of the modulated amplitude and the associated instantaneous frequencies that appear to represent both interwaves and intrawaves;  reconstruction of the signal and the Hilbert spectrum (energy-frequency) and the multi-component frequency-time relations.

The HHT technique for analyzing data consists of two components: a decomposition algorithm called Empirical Mode Decomposition (EMD) and a spectral analysis tool called Hilbert spectral analysis.

Empirical Mode Decomposition

The key part of the approach HHT is that any complicated dataset can be decomposed with the EMD method into a finite and small number of Intrinsic Mode Functions (IMFs), which represent different scales of the original time series and physically meaningful modes. An IMF is defined as a function having the same number of extrema and zero-crossings. It has also symmetric upper and lower envelopes defined by the local maxima and minima respectively.

The IMFs represent different scales of the original time series and physically meaningful modes. Indeed, the EMD divides the time series into a series of modes, but unlike a Fourier transform where each component has a constant frequency and a cosine or sine shape, the EMD allows each mode to have a time-dependent frequency and amplitude. In contrast to almost all the previous methods, the EMD works directly in temporal space rather than in the corresponding frequency space. With the definition of the IMFs, one can then decompose any function   t x through a sifting process:  First identify all the local extrema (minima) and then connect all the local maxima (minima) by a cubic spline line to form the upper (lower) envelope. The envelopes should encompass all the data between them.  Calculate 1 m , the running mean of the envelope.

 Treat the difference   1 1 h m t x  
as the data in the next iteration and repeat the above sifting processing until

1 h becomes an IMF. o ) 1 ( 1 ) 1 ( 1 1 h m h   .
o After k times, the first IMF is derived:

1 ) ( 1 ) ( 1 ) 1 ( 1 c h m h k k k     . The
  1 1 c t x r  
and carry out the same process described above.  Repeat the procedure until n r becomes a constant, a monotonic function or a function with only one maximum and one minimum from which no more IMF can be extracted:

o 2 2 1 r c r   o … o n n n r c r   1
 By summing all the equations, we can get a decomposition with n IMFs and a residue n r as follows:

          n j n j t r t c t x 1 .
Due to a dyadic filter bank property of the EMD algorithm (Flandrin et al., 2004;[START_REF] Huang | An amplitude-frequency study of turbulent scaling intermittency using Hilbert spectral analysis[END_REF][START_REF] Wu | A study of the characteristics of white noise using the empirical mode decomposition method[END_REF], usually in practice, the number of IMFs modes is less than

  N 2 log
, where N is the length of the data set. Note that the sifting process serves two purposes: (1) to eliminate background waves on which the IMF is riding and (2) to make the wave profile more symmetric. Unlike Fourier based methodologies such as Fourier analysis and wavelet transform, the EMD method does not define a basis a priori (Flandrin and Gonçalvès, 2004;[START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF][START_REF] Huang | A new view of nonlinear water waves: the Hilbert spectrum[END_REF] and the IMFs obtained by sifting processes constitute an adaptive basis. Therefore, EMD is very suitable for nonstationary and nonlinear time series analysis [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF][START_REF] Huang | A new view of nonlinear water waves: the Hilbert spectrum[END_REF]. One way to express the nonstationarity is to find instantaneous frequency and instantaneous amplitude. This was the reason why Hilbert spectrum analysis was included as a part of HHT.

HHT has already been applied in marine sciences. For example, [START_REF] Dätig | Performance and limitations of the Hilbert-Huang transformation (hht) with an application to irregular water waves[END_REF] applied HHT to show excellent correspondence between simulated and recorded nonlinear waves. [START_REF] Schmitt | Analysis of velocity fluctuations and their intermittency properties in the surf zone using empirical mode decomposition[END_REF] applied the HHT method to characterize the scale invariance of velocity fluctuations in the surf zone. The EMD scheme was used in studying sea level rise [START_REF] Ezer | Gulf Stream's induced sea level rise and variability along the U.S. mid-Atlantic coast[END_REF]. [START_REF] Yin | Coastal-trapped waves in the East China Sea observed by a mooring array in winter 2006[END_REF] also applied the method and identified three kinds of low-frequency waves using some observations in the coastal water of the East China Sea.

In this paper, the EMD algorithm has been applied to the Marel Carnot data sets acquired from 2005 to 2009 and already presented in section 2 (Presentation of the experimental database). The analysis for the salinity records is presented below.

EMD of the salinity time series

After EMD decomposition of the salinity data sets, 19 IMFs modes are obtained plus the residual, as shown in Fig. 14 and Fig 15 . The time scale is increasing with the number of the IMF mode, the first IMF thus corresponding to the highest frequency. The residual from EMD algorithm has been recognized as the trend of the given data [START_REF] Moghtaderi | Trend filtering: empirical mode decompositions versus l1 and Hodrick-Prescott[END_REF][START_REF] Wu | On the trend, detrending, and variability of nonlinear and non-stationary time series[END_REF].

From EMD, several tidal waves are identified through the IMFs: semidiurnal ( 6IMF , 12 hours such as M2), fourth-degree diurnal (

    ) ( sin 2 24 ) ( sin 2 1 1 hours days f inertial        (Monin, 1990).
Here, the inertial wave is

    51 sin 2 24 44 . 5 1 
hours. The period of

7

IMF is close to the inertial period but the bottom friction induced by the strong tidal currents is likely to dissipate these inertial oscillations and this period may be due the non linear interaction between tidal and high frequency atmospheric forcing. Table 4 summarizes the mean period of the IMF modes and some tidal waves obtained after EMD decomposition of the salinity. Furthermore, the contribution of each IMF to the total energy is measured by the variance. Table 4 gives also the contribution of each IMF to the total energy (measured by the variance), percentage of total and cumulative percentage contribution to the sum of the variances. For the high frequencies, the firstseventh IMFs account for nearly 25% of the total energy. Thus, the main part of energy is not carried by the very high frequency components but by intermediate and low ones, especially the annual wave 17 IMF . It can be seen that the biggest contribution comes from four most energetic IMFs: the annual wave ( 17 IMF ), the semiannual wave ( 16IMF ), the semimonthly wave 13 IMF and the monthly wave 14 IMF . These four IMFs account for over 52 % of the total energy.

While the Fourier expansion would require tens of modes to represent the whole data, the EMD method decomposes the time series into only 19 IMFs plus the residual. When all the IMFs are added back successively, we notice that all the energy is recovered, as shown in all the cases in [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF][START_REF] Huang | A new view of nonlinear water waves: the Hilbert spectrum[END_REF].

IMFs

Mean period (days) To analyze the variability, both HHT and Conitunous Wavelet Transform (CWT) can be used. Comparaisons of HHT and wavelets are investigated in the next subsection.

Continuous Wavelet Transformation (CWT)

The wavelet transform is widely used as a time-frequency analysis technique to deal with nonstationary signals. The choice of the mother wavelet is usually dependent on the type of data to deal with [START_REF] Torrence | A practical guide to wavelet analysis[END_REF]. HHT on the other hand does not require any convolution of the signal with a predefined basis function or mother wavelet. The process of decomposition is data-driven.

The CWT is a consecutive pass-band filter through time series, recognized as a very successful approach for the analysis of time series [START_REF] Addison | The Illustrated Wavelet Transform Handbook: Introductory Theory and Applications in Science, Engineering, Medicine and Finance[END_REF][START_REF] Avdakovic | Correlation between Air Temperature and Electricity Demand by Linear Regression and Wavelet Coherence Approach: UK, Slovakia and Bosnia and Herzegovina Case Study[END_REF][START_REF] Torrence | A practical guide to wavelet analysis[END_REF]. However, similar to the Fourier transform and in contrast to the EMD algorithm, the wavelet transform requires evenly-spaced data. Therefore, the data have to be interpolated.

Mathematically, a wavelet transformation decomposes a signal into some elementary functions derived from a mother wavelet by dilation and translation, where denotes the position (translation) and (

) the scale (dilation) of the wavelet [START_REF] Lau | Climate Signal Detection Using Wavelet Transform: How to make a Time Series Sing[END_REF]. The wavelet transform of the signal about the mother wavelet is defined as the convolution integral . There are many mother wavelets available in the literature. However, in geosciences the most popular wavelet is the Morlet wavelet, which is given by where is the non-dimensional time and is the wavenumber [START_REF] Torrence | A practical guide to wavelet analysis[END_REF]. The wavenumber is the number of oscillations within the wavelet itself.

In Fig. 16, the wavelet transform has been applied to the four time series. Many high frequencies are detected for temperature, salinity and turbidity records. More variations are visible for the turbidity time series, which confirms that the data sets are highly non-stationary compared to the other data. This analysis confirms also the previous results presented in 4.1, where the performed harmonic analysis allows to predicting only 8.3% of the original signal's variance. As explained before, this is due to biological contents induced by the plankton in the spring, high-frequency effects linked to storms and a cycle of tidal streams (suspension and advection).

However, if we look at the CWT of the sea level in Fig. 16, the graph shows horizontal clusters of peaks at 365 days and 0.5 day, which corresponds to a strong annual cycle and the semidiurnal wave. The powerspectrum also revels other significant periods such as the diurnal wave (1 day), the semimonthly and monthly waves; although these are sometimes not very clear.

As the signal content is now known in the frequency domain, the focus will shift to the application of an EMD based Time Dependent Intrinsic Correlation (TDIC), in order to study the possible links between physical observations. Figure 15. The Morlet wavelet spectrum of the data. The period in the y-axis is expressed in days.

Section 6: Cross correlation methods

The classical global expression for the correlation (defined as the covariance of two variables divided by the product of the standard deviation of the two variables) assumes that the variables should be stationary and linear. Applied to nonstationary time series, the cross correlation information may be altered and distorted. The limitations of the correlation coefficient are also obvious: it is unable to provide local temporal information, and it cannot distinguish the main cycles from noise when measuring correlation. Many scientists tried to address the problem of nonsense correlations through different ways.

An alternative is to estimate the correlation coefficient by means of a timedependent structure. For example, [START_REF] Papadimitriou | Local correlation tracking in time series[END_REF] applied a sliding window to localize the correlation estimations. Rodo and Rodriguez-Aria (2006) developed the scale-dependent correlation technique. Although these methods detected the correlation between two nonstationary signals by computing the correlation coefficient in a local sliding window, the main problem is to determine the size of this window. Recently, [START_REF] Chen | The time-dependent intrinsic correlation based on the empirical mode decomposition[END_REF] introduced an approach based on EMD. They proposed to first decompose the nonlinear and nonstationary data into their IMFs, then use the instantaneous periods of the IMFs to determine an adaptive window and finally compute the time dependent intrinsic correlation coefficients. [START_REF] Huang | Time dependent intrinsic correlation analysis of temperature and dissolved oxygen time series using empirical mode decomposition[END_REF] used TDIC to analyze temperature and dissolved oxygen time series obtained from automatic measurements in a moored buoy station in coastal waters of Boulogne-Sur-Mer (France). In the following, we give an insight into the TDIC method.

Time Dependent Intrinsic Correlation (TDIC)

The correlation between two data sets is considered here. Suppose the two time series   [START_REF] Huang | Time dependent intrinsic correlation analysis of temperature and dissolved oxygen time series using empirical mode decomposition[END_REF] or by considering the Fourier energy weighted mean frequency, i.e.,  

       df f X f df f X t T i j i j i j 2 2 where   f X i j is the
Fourier power spectrum of each IMF mode. Then, at time inst t , the sliding window is given by

                      2 , max , 2 , max 2 1 2 1 t T t T a t t T t T a t t j j in st j j in st win
, where a is any positive number. This window is different from classical sliding windows: it is based on the maximum of two instantaneous periods

      t T t T j j 2 1
, max and thus it is adaptive.

Results for cross analyses

The focus will shift now to illustrate the cross correlations between temperature and sea level time series. We consider a sample of one year for these records between 24 th May 2006 and 24 th May 2007. First, the time series are not correlated since the global cross correlation coefficient is only 0.025.

Then, the EMD algorithm is applied to the data sets for the same time period between 24 th May 2006 and 24 th May 2007. For the temperature records, there are 12 IMFs modes with one residual, while the sea level is decomposed into 6 IMFs with one residual. Thus, the data are represented in a multiscale way (Flandrin and Gonçalvès, 2004;[START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF][START_REF] Huang | Arbitrary-order Hilbert Spectral Analysis: Definition and Application to Fully Developed Turbulence and Environmental Time Series[END_REF] and the IMFs can be used for multiscale correlation. For instance, let us consider the IMFs modes with a mean period of 12 hours (semidiurnal) : Fig. 17 displays the measured TDIC and shows rich patterns at small sliding window. We note a decorrelation of the TDIC and we observe less rich patterns with the increase of the window size (as the variable a increases). Although the global cross correlation is too small, the TDIC detects several periods of high correlation between the modes. To focus on these periods, Fig. 18 shows a zoom of the IMFs modes between 15 th july and 15 th August 2006. These are positively correlated with each other on some portions and negatively correlated on others; showing rich dynamics. The red circles point out some portions where the IMFs co-vary. This means that sometimes the temperature records are highly influenced by the tidal cycles. For example, the correlation coefficient between

Section 7: Conclusions and perspectives

In the era of "Big Data", the time series are longer and longer. Time series analysis in marine sciences is becoming a challenge at the crossroads of marine observatories' data. Indeed, the observatories allow more and more marine data to be safely archived. They contribute to understand and forecast the most crucial marine processes. Therefore, it is necessary to analyze this increasing data flow in an optimal way. In environmental sciences, adequate and specific methods are needed to analyze the complex time series which are nonstationary and reflecting nonlinear processes.

To explore these methods, we have analyzed here four time series recorded at high frequency in the coastal waters of Boulogne-sur-Mer, using the MAREL Carnot monitoring station: sea level, salinity, turbidity and temperature data sets. In order to perform adequate statistical and spectral analyses, it is necessary to know the nature of the considered time series. First, a preliminary analysis of stationarity, based on the graphs of the time series as well as the correlograms, is carried out. In contrast to the sea level time series which is stationary, the graphs of temperature, turbidity and salinity are characterized by a non stationarity in mean and variance. The series are also perturbed by effects such as the tide and the earth's rotation. Regarding the corresponding correlograms, these these show strong persistence since there are significant spikes at several lags, indicating that the time series are nonstationary. Then, the stationnarity of the series and the occurrence of unit-root are addressed with the ADF tests. For the salinity data, the time series is nonstationary. It has a unit-root with drift and a deterministic trend. Similar results are obtained for temperature and turbidity data.

After stationarity issues, an adequate power spectral analysis is performed. The harmonic analysis involves the search for the amplitude and phase of certain known frequencies, used for the calculation of the astronomical tide's generator potential.

Apart from the sea level where the harmonic analysis is able to predict 99.7% of the

A C C E P T E D

M A N U S C R I P T original signal's variance, such analysis is not relevant for the other time series. Different spectra for the turbidity and salinity are investigated. They reveal a very complex dynamics, with a turbulent-like stochastic behaviour at high frequencies, whereas the dynamics at low frequencies are dominated by deterministic forcings.

Compared to the modified periodogram, the Welch method improves the signal-to-noise ratio of the spectrum versus a loss of resolution. Since the Blackman-Tukey method requires evenly-spaced data, we have interpolated the time series. However, interpolation introduces numerous artifacts to the data, both in the time and the frequency domain. Thus, the Lomb-Scargle algorithm is adapted for unevenly-spaced data and has been used as an alternative method. To know the limits of this algorithm, we have investigated the time series of salinity over the year 2005 and we have intentionally increased the rate of missing values by discarding randomly data from the original time series. The Lomb-Scargle powerspectra show that, beyond 50% of missing measures, few significant frequencies are detected, several seasonalities are no more visible, and even a whole range of high frequency disappears progressively.

The authors also explored an other type of interesting direction of study since the Blackman-Tukey spectral analysis is not adapted to nonlinear and nonstationary data sets. An adaptive data analysis method, called EMD, was developed to accommodate the variety of data generated by nonlinear and nonstationary processes. In this paper, EMD was used to analyze the fluctuations and the correlations between different time series. For example, the salinity observations were decomposed with the EMD method into a finite and small number of IMFs (19 IMFs) which represent different scales of the original time series. These are also physically meaningful modes since a group of tidal waves were detected. Compartsions with CWT are addressed. Furthermore, the decomposition into modes helps also to estimate how correlations vary among scales. In fact, the authors apply a recent methodology, based on EMD and called TDIC, in order to display patterns of correlations at different scales for different IMFs modes. Finally, to illustrate this approach, cross correlations between temperature and sea level time series are considered. For the semidiurnal wave, the TDIC is displayed and allows to detecting several periods of high correlation.

We can mention here some perspectives opened by this work. Techniques of filtering will be investigated to study the high frequency or to concentrate on seasonalities. Classical Butterworth filters can be implemented or demerliac filters, already used by the "Service Hydrographique et Océanique de la Marine" (S.H.O.M) to calculate the daily mean sea level. Moreover, we have studied each parameter separately. A development of this approach will be to study the couplings between some parameters. Such coupling will be studied considering joint or conditional moments, joint probability density function, and joint spectral analysis.

More generally, time series analysis methods could be implemented in an automated way in the environmental databases in order to facilitate the data discovery as well as the investigation tools and techniques.

Figure 1 .

 1 Figure 1. Temperature time series over 5 years (a), showing strong annual and seasonal cycles. A 15-day portion indicates a periodic component associated to the tide together with stochastic fluctuations (b). The red curve represents a downward trend.

Figure 2 .Figure 3 .

 23 Figure 2. Turbidity time series over 5 years (a). A 15-day portion of the data sets indicates a strong periodic component associated to the tide together with stochastic fluctuations (b).

Figure 4 .

 4 Figure 4. PDFs of a) salinity, b) turbidity, c) temperature and d) sea level.

  (b) shows how the temperature values are declining and represents a downward trend for the considered 15 day portion.  An increase or decrease in variance by slices, characterizing a non stationarity in variance. For example, Fig. 5 displays the variance by slices for the temperature time series over the year 2007. The continuous and random increase and decrease in variance by slices illustrates the non stationarity of the time series. Here, the variability at time scale of the week is mainly induced by the variability of current (tides and wind induced circulation) rather than by the seasonal cycle of ocean-atmosphere fluxes.  A regularly repeated phenomenon from the season to the tidal period (i.e., from about one hour to six months), which lets suppose the presence of seasonality.

Figure 5 .

 5 Figure 5. Variance by slices (window size of one week) for the temperature time series over the year 2007.

Figure 6 .

 6 Figure 6. A 30-day portion of the sea level data sets indicates a strong periodic component associated to the tide.

Figure 7 .

 7 Figure 7. Correlogram of sea level data sets.

Figure 8 .

 8 Figure 8. Correlogram of turbidity data sets before and after differencing.

  Figure 9. Use of t_tide toolbox for the harmonic analysis of sea level data: Spectral estimates before and after removal of tidal energy.

Figure 10 .

 10 Figure 10. Two powerspectra of the turbidity data sets: modified periodogram with a Hamming window versus Welch powerspectrum.

Figure 11 .

 11 Figure 11. Lomb-Scargle powerspectrum of the salinity data over the year 2005.

Figure 12 .

 12 Figure 12. Comparing different Lomb-Scargle powerspectra for the salinity while increasing the rate of missing data.

Fig. 12

 12 Fig. 12 displays the corresponding Lomb-Scargle powerspectra. The more data are discarded, the less energetic spikes are observed in the powerspectrum. Table3provides the number of significant frequencies above the 95% confidence level, dropping down from 131 for the original data of salinity to only 14 when 90% of the data are randomly eliminated. Several seasonalities are no more visible in the spectra and even a whole range of high frequency disappears progressively. Beyond 3 months of discarded data, the Lomb-Scargle algorithm really suffers from a high rate that approaches 50% of missing measures.

term 1 c

 1 contains the highest frequency component of the original signal. o Separate 1 c from the rest of the data:

IMF

  may represent the inertial wave with a theoretical period of

Figure 13 .

 13 Figure 13. The first-tenth IMF modes for salinity. The time scale is increasing with the mode index.

5

  IMF for temperature and 1 IMF for sea level data sets. The global correlation coefficient is 0.09 for the corresponding IMFs modes.

  time series of temperature and sea level respectively, is 0.59 from 15 th to 20 th July 2006, 0.61 from 4 th to 7 th August 2006, which is much higher compared to the global correlation coefficient between these IMFs, equal to only 0.09.

Figure 16 .

 16 Figure 16. The measured TDIC for the 12-hour mean period obtained after EMD decomposition of temperature and sea level between 24 th May 2006 and 24 th May 2007.

Table 1 .

 1 Statistical properties of the time series.

		Mean	Variance	Median	Mode	Skewness
	Salinity (µmol/L)	33.43	0,79	33.56	33.75	-2.49
	Turbidity (NTU)	12.33	203,61	7.70	0	2.73
	Temperature (°C)	12.65	20,86	18.30	12.70	-0.03
	Sea level (m)	4.91	4.82	4.85	2.81	0.04

  suggest the following data dependent lag length selection procedure that results in stable size of the test and minimal power loss. If the absolute value of the t-statistic for testing the significance of the last lagged difference is greater than 1.6 (see[START_REF] Ng | Unit Root Tests in ARMA Models with Data-Dependent Methods for the Selection of the Truncation Lag[END_REF] then set max p p  and perform the unit root test. Otherwise, reduce the lag length by one and repeat the process.  A common rule of thumb for determining max p , suggested by Schwert

	 First, set an upper bound max p for	p .
	 Next, estimate the ADF test regression with	p 	p	max	.
	

Table 3 .

 3 Number of significant frequencies above the 95% confidence level detected by the Lomb-Scargle method.

	Time series  significant frequencies
	Original data	131
	Case 1	90
	Case 2	79
	Case 3	53
	Case 4	14

Table 4 .

 4 Contribution of each IMF to the total variance for

			Tidal waves	Variance	% of total	Cumulative %
	IMF1	0,06 = 1,43 hours		0.0292	3.04	3.04
	IMF2	0,11 = 2,66 hours		0.0247	2.57	5.61
	IMF3	0,18 = 4,44 hours		0.0359	3.74	9.35
	IMF4	0,26 = 6,20 hours	M4	0.0423	4.40	13.76

2

Temp time series

Acknowledgments

The authors would like to thank the Région Bretagne for financial support of the postdoctoral fellowship (SAD MASTOC). They also thank the organizing committee of the MAREL Carnot symposium held in June 2014, Boulogne-Sur-Mer, France. Finally, we wish to thank the Editor, the Associate Editor, and the Referees for helpful comments.


Guillen Jordi
email: jordiguillen@hotmail.com

Cheilari Anna

Damalas Dimitrios

Barbas Thomas

Oil for Fish: An Energy Return on Investment Analysis of Selected European Union Fishing Fleets

Keywords: energy efficiency, energy return on investment (EROI), fisheries, fishing fleet, fuel consumption, industrial ecology

World food production has increased substantially in the past century, thanks mostly to the increase in the use of oil as input in the production processes. This growing use of fossil fuels has negative effects, both on the environment and the production costs. Fishing is a fuel consuming food production activity, and its energy efficiency performance has worsened over time. World-wide fisheries are also suffering from overexploitation, which contributes to the poor efficiency performance, adding more pressure and criticism on this economic activity. In this paper we analyzed the energy efficiency performance of more than 20,000 European Union (EU) fishing vessels for the period 2002-2008, using the edible energy return on investment (EROI) indicator. The vessels analyzed, grouped in 49 different fleets, represented 25% of the vessels and 33% of the landings of the EU fishing sector. These EU fishing fleets' average EROI for 2008 was 0.11, which translates to an energy content of the fuel burned that is 9 times greater than the edible energy content of the catch. Hence, the significance of this study arises from the use of time-series data on a relevant part of the EU fleet that showed stable or even slight improvements on the EROI over time. Moreover, results showed that the energy efficiency of the different fleets varied significantly (from 0.02 to 1.12), mainly depending on the fishing gear and the vessel length. The performance of the most efficient fleets, such as large pelagic trawlers and seiners, was comparable to many agricultural production activities. The plausible drivers behind these trends are further considered.

Introduction

World food production has increased substantially in the past century, as has calorie intake per capita [START_REF] Nellemann | The environmental food crisis -The environment's role in averting future food crises[END_REF]. But this has happened at the cost of increasing the use of oil as input in the food production process. "Eating Oil" is the title of a book published in 1978 that investigates the extent to which food supply in industrialized countries relies on fossil fuels [START_REF] Green | Eating Oil -Energy Use in Food Production[END_REF]). Not surprisingly, most studies on energy efficiency for food production were realized during the 70s and 80s, after the first oil crisis in 1973 (i.e. [START_REF] Leach | Energy and Food Production[END_REF][START_REF] Pimentel | Food, energy, and society[END_REF][START_REF] Pimentel | Energy efficiency of farming systems: organic and conventional agriculture[END_REF][START_REF] Cleveland | Energy and the U.S. economy: a biophysical perspective[END_REF][START_REF] Hall | Energy and Resource Quality: The Ecology of the Economic Process[END_REF]). Nowadays, food supply worldwide is even more dependent on oil, and a lot more studies have populated the scientific literature (i.e. [START_REF] Dutilh | Energy consumption in the food chain[END_REF][START_REF] Tyedmers | Fisheries Impacts on North Atlantic Ecosystems: Catch, Effort and National/Regional Datasets[END_REF][START_REF] Tyedmers | Fisheries and energy use[END_REF]Pelletier and Tyedmers 2007;[START_REF] Pelletier | Environmental performance in the US broiler poultry sector: Life cycle energy use and greenhouse gas, ozone depleting, acidifying and eutrophying emissions[END_REF][START_REF] Schau | Energy consumption in the Norwegian Fisheries[END_REF][START_REF] Fet | A Framework for Environmental Analyses of Fish Food Production Systems Based on Systems Engineering Principles[END_REF], Pelletier et al. 2010a, 2010b).

The energy return on investment (EROI) is the ratio of the energy delivered by a process to the energy used directly and indirectly in that process [START_REF] Leach | Energy and Food Production[END_REF][START_REF] Pimentel | Food, energy, and society[END_REF][START_REF] Pimentel | Energy efficiency of farming systems: organic and conventional agriculture[END_REF][START_REF] Cleveland | Energy and the U.S. economy: a biophysical perspective[END_REF][START_REF] Hall | Energy and Resource Quality: The Ecology of the Economic Process[END_REF]). The EROI is generally used for measuring the efficiency of different technologies and sources of energy, expressed as the ratio of the gross energy return to the amount of energy invested [START_REF] Mulder | Energy return on investment: toward a consistent framework[END_REF].

Other names used for this concept are energy profit ratio, energy efficiency, energy gain [START_REF] Tainter | Resource transitions and energy gain: contexts of organization[END_REF], net energy [START_REF] Odum | Energy, ecology, and economics[END_REF]) and energy payback [START_REF] Keoleian | Application of life cycle energy analysis to photovoltaic module design[END_REF].

Although originally the EROI indicator was a measure of the efficiency of different energy sources, the EROI ratio has also been used in food production processes to compare the energy extracted from the edible content of the food and the energy in the form of input requested for their production. The EROI analysis provides a quantitative result that can be easily compared with other similar calculations [START_REF] Murphy | Order from Chaos: A Preliminary Protocol for Determining the EROI of Fuels[END_REF]. Therefore, the EROI indicator provides an objective measurement of the nourishment that can be produced by unit of energy spent (i.e. fuel combustion). This allows to compare the energy efficiency of different fleets without dependence on landing weights (that can provide partial vision due to large non edible contents, i.e. in mussels) or on their monetary values, and consequently, results can be objectively compared worldwide. It also allows comparison of the fishing sector with other sectors, such as livestock or aquaculture (i.e. [START_REF] Tyedmers | Fisheries and energy use[END_REF][START_REF] Tyedmers | Fueling global fishing fleets[END_REF]. Indeed, the ratio has been used in agriculture (i.e. [START_REF] Pimentel | Food, energy, and society[END_REF][START_REF] Pimentel | Energy efficiency of farming systems: organic and conventional agriculture[END_REF]), cattle farming (i.e., Pelletier 2008), fisheries (i.e. [START_REF] Tyedmers | Fisheries Impacts on North Atlantic Ecosystems: Catch, Effort and National/Regional Datasets[END_REF][START_REF] Tyedmers | Fisheries and energy use[END_REF] and aquaculture (i.e. Pelletier and Tyedmers 2007).

Fisheries have often been considered one of the less efficient food production activities [START_REF] Wilson | Fuel and financial savings for operators of small fishing vessels[END_REF]; however, this generalization is not supported by the work of several authors (i.e. [START_REF] Tyedmers | Fueling global fishing fleets[END_REF]. It has been estimated that fishing fleets around the world burn between 38 and 42.5 million tons of oil per year, accounting for about 1.2% of global oil consumption [START_REF] Tyedmers | Fueling global fishing fleets[END_REF]World Bank and FAO 2008). Moreover, a general tendency of energy performance becoming poorer over time is evident from a number of fisheries around the world (i.e. [START_REF] Tyedmers | Fisheries Impacts on North Atlantic Ecosystems: Catch, Effort and National/Regional Datasets[END_REF][START_REF] Tyedmers | Fisheries and energy use[END_REF]. Fisheries are suffering from poor profitability, mainly due to increasing fuel costs, as well as the increasing stocks overexploitation (i.e. [START_REF] Csirke | Global production and state of marine fishery resources[END_REF][START_REF] Abernethy | Fuelling the decline in UK fishing communities?[END_REF]. This implies that in order to extract the same amount of fish it is necessary to consume more oil. The energy efficiency of the different existing fleet segments varies significantly. Both the fishing gear/technique and the vessel length that characterize each fleet segment affect the energy performance of the fleets [START_REF] Thrane | Energy Consumption in the Danish Fishery: Identification of Key Factors[END_REF][START_REF] Ziegler | Environmental life cycle assessment of Norway lobster (Nephrops norvegicus) caught along the Swedish west coast by creels and conventional trawls -LCA methodology with case study[END_REF][START_REF] Driscoll | Fuel use and greenhouse gas emission implications of fisheries management: the case of the New England Atlantic herring fishery[END_REF].

Because fuel combustion has adverse environmental impacts and also raises production costs, it is advisable to reduce fuel consumption, or at least consume it in the most efficient way possible. This is also valid for fisheries, where fuel often is the most important cost (FAO 1995;[START_REF] Dahou | Socio-economic and Environmental Impact of Senegalese Fishery Support Mechanisms[END_REF][START_REF] Reddy | Economic analysis of artisanal fisheries in Fiji: issues of profitability and sustainability[END_REF][START_REF] Ferm | ITQs, ageing boats, and the price of fish. Profitability and autonomous adjustment in the South East Trawl Fishery[END_REF][START_REF] Sumaila | Rebuilding Hong Kong's marine fisheries: an evaluation of management options[END_REF]STECF 2010;[START_REF] Sumaila | Subsidies to high seas bottom trawl fleets and the sustainability of deep-sea demersal fish stocks[END_REF]. Under the Europe 2020 Strategy for smart, sustainable and inclusive growth (EC 2010), EU Member States have committed themselves to reducing greenhouse gas emissions (GHG) by 20% compared to 1990 levels , increasing the share of renewable sources in the EU's energy mix to 20%, and achieving a 20% increase in energy efficiency by 2020.

In this paper, by using public and comprehensive time series data we analyze the energy profiles of different EU fishing fleets with the edible Energy Return on Investment ratio. By doing this we do not just estimate the evolution of the EROI for different EU fleets, but we are in a position to identify the most efficient fishing techniques from the nourishment perspective.

Besides the obvious economic incentives, the use of more efficient fishing gears can contribute to the environmental sustainability.

Materials & Methods

Data included in the 2010 Annual Economic Report (AER) on the EU fishing fleet (STECF 2010) was used for this analysis. The AER provides the most complete compilation of statistics on the economic performance of the EU Member States fishing fleets. The report contains an economic and structural analysis of the European Union fishing fleets with data collected under the Data Collection Framework (DCF) (EC 2008). To the authors knowledge, there is no other publication worldwide that periodically reports fisheries production and cost data at such levels of detail.

For the current analysis, we considered only data for those fleet segments (a group of vessels belonging to the same length class and having the same predominant fishing gear during the year) for which full time series data for all the relevant variables (fuel consumption and weight of landings per species by fleet) were available. The analysis was performed at the fleet segment level and data from different countries were aggregated together. All main fleet segments were covered by the study: trawls and seines (demersal and pelagic), dredges, hooks, nets, pots and traps. The fleet segmentation follows the EU segmentation (EC 2001) by fishing technique and vessel length class (the codes used can be found below the table 2 and figures).

Fuel consumption and landing weights by species covering the years 2002 to 2008, were used for the calculation of the EROI. Data used originated from more than 20,000 fishing vessels grouped in 49 fleet segments of seven countries (Belgium, Germany, Denmark, France, Italy, Netherland and Sweden), representing 25% in terms of number of vessels and 33% in terms of landings' weight of the EU-27 for the time period studied.

For the calculation of the edible EROI we divided the energy extracted from the edible content of the fish i by the energy consumed in the production process (equation 1). The former is calculated by multiplying the landings from each fleet by the energy content and the edible content of each fish species. The energy content of fresh seafood per species was assigned according to [START_REF] Dorosz | Tabla de calorias[END_REF] and when no value for a species was reported, the value of 1,116 Kcal/kg was used as the average energy content [START_REF]La alimentación en España[END_REF]. Because not all the parts of a fish are edible, different edible fractions by species have been applied (Tyedmers, personal communication, 2010 ii ). Therefore, the boundary of this study is the EROI 1,d as classified in [START_REF] Murphy | Order from Chaos: A Preliminary Protocol for Determining the EROI of Fuels[END_REF] iii .

The edible EROI for a fleet on a given year is calculated according to the following formula:

∑ (1)
Where LW stands for the landings weight in Kg, FC for the fuel consumption in liters, EC for the energy content in Kcal per kg, EdC for the edible content of the fish expressed as a percentage, f indexes the fishing fleet, s the species and t the year. One liter of fuel is estimated to contain around 38.66 MJ [START_REF] Schau | Energy consumption in the Norwegian Fisheries[END_REF], and one Kcal is approximately 4.187 kJ.

Only the energy content of the edible content of fresh seafood has been accounted for in this study. This implies that co-product allocation (the use of fish parts not suitable or edible for human consumption, i.e. use in animal feeds) has not been taken into account. This is because data on the different consumption and processing patterns depend mostly on country/region and fleet (i.e. consumer preferences to buy the whole fish or fillets), as well as the processes coproducts are exposed. Consequently this data is rather difficult to obtain for such a diverse number of fleets. Therefore, the EROI estimated in this study corresponds to the EROI measure related to human-edible food energy return on industrially-mediated energy investment specified in [START_REF] Pelletier | Energy intensity of agriculture and food systems[END_REF].

Indirect energy use was excluded from the analysis due to the lack of data and only direct fuel consumption was used to calculate the relevant indicator. This is not a major shortcoming as direct fuel use in fisheries account between 75 and 90%, regardless of the fishing gear used or the species targeted [START_REF] Ziegler | Life Cycle Assessment of frozen cod fillets including fishery-specific environmental impacts[END_REF][START_REF] Tyedmers | Fisheries and energy use[END_REF][START_REF] Hospido | Life cycle environmental impacts of Spanish tuna fisheries[END_REF][START_REF] Pelletier | Energy intensity of agriculture and food systems[END_REF]. Depending on the fishery, the remaining 10 to 25% is generally composed of other energy inputs associated with vessel construction and maintenance, the provision of labor, fishing gear, bait and ice, fish processing and transportation to the market [START_REF] Rawitscher | Energy cost of nutrients in the American Diet[END_REF][START_REF] Watanabe | An estimation of direct and indirect energy input in catching fish for fish paste products[END_REF][START_REF] Tyedmers | Salmon and sustainability: The biophysical Cost of Producing Salmon through the Commercial Salmon Fishery and the Intensive Salmon Culture Industry[END_REF][START_REF] Andersen | Transport of fish from Norway: energy analysis using industrial ecology as the framework[END_REF][START_REF] Ziegler | Life Cycle Assessment of frozen cod fillets including fishery-specific environmental impacts[END_REF][START_REF] Fet | A Framework for Environmental Analyses of Fish Food Production Systems Based on Systems Engineering Principles[END_REF]. Moreover, we are interested in the comparison between fishing fleets and especially in trends in energy use and efficiency over time. Thus, the exclusion of indirect energy should not bias the results, but it should be noted that the estimated indicators are conservative with respect to the studies that consider also indirect energy use, by a magnitude of 10 to 25%.

Results

For the analyzed fleet segments, the landings weight and the fuel consumption have been decreasing over time, reaching in 2008 almost two thirds of the values recorded in 2002 (Figure 1). The decrease in the fuel consumption is almost monotonic, whereas in the weight of landings the decreasing trend appears more scaled. The distribution of the landings and fuel consumption by fishing technique for the studied fishing fleets for the year 2008 are presented in Figures 2 and3.

An overview of the EROI estimated for the different fleets is presented in Table 1. The weighted average edible EROI for 2008 was 0.11, ranging between 0.09 and 0.11 for the period of the analysis. That is to say that the average energy content of the fuel burned by the studied European fleets was nine times greater than the edible energy content of the resulting catch in 2008.

The different EU fishing fleets analyzed demonstrated important differences in the EROI levels (see table 2). Gears using hooks (0.04 in 2008), and beam trawlers between 24 and 40 meters and larger than 40 meters (0.02) had the lowest EROI values, followed by polyvalent passive gears (0.04), demersal trawls and seiners (0.05), and driftnets and fixed nets (0.07).

Pelagic trawlers and seiners larger than 40 meters (1.12) had the highest EROI values estimated. High values were also observed for dredges lower than 12 meters (0.20) and small vessels equipped with passive gears (0.38). These important differences in the ratio between fishing gears and vessel lengths was in accordance with [START_REF] Thrane | Energy Consumption in the Danish Fishery: Identification of Key Factors[END_REF], who presented different energy consumption levels for nine Danish fleets. The edible EROI values were in line with those reported in [START_REF] Tyedmers | Fisheries and energy use[END_REF] for 29 North Atlantic and Pacific fisheries from several studies. The values varied from 0.019 for the trawl fishing targeting shrimps and the trawl fishing targeting flatfish to an EROI of 0.56 for the purse seine fishing for herring and mackerel. Similarly, beam trawlers between 24 and 40 meters and beam trawlers larger than 40 meters had an EROI for 2008 of 0.020 and 0.016, respectively. These fleets target mainly sole, plaice (flatfishes) and shrimps which are species with high commercial value. On the other hand, pelagic trawlers and seiners larger than 40 meters, targeting low-value species such as herring, mackerel, sand-eels and sprat, had an EROI of 1.12.

From the results it can also be concluded that for each fishing technique, there was an optimal vessel length class, as can be seen in table 2. For pelagic and demersal trawlers and seiners, the larger vessels were more energy efficient; whereas the opposite was occurring for beam trawlers and dredges.

Discussion

The reduction of greenhouse gas emissions and the efficient use of resources have become important political objectives in the agenda of the European Union. Clear incentives to reduce fuel consumption, or at least consume it in the most efficient way possible, start to appear in many sectors, including fisheries and food production.

Continued technological advances in fishing fleets increase catch efficiency directly affecting fish stocks. It is commonly perceived that with the advent of modern fish-catching technologies, fleets would become more catch efficient. This phenomenon known also as technological creep or, simply, creep of the fleet is usually positively related to the increase in skipper skills, investments in auxiliary equipment, more efficient gear and materials, replacement of old vessels by new ones and, to a lesser extent, upgraded engines [START_REF] Rijnsdorp | Partial fishing mortality per fishing trip: a useful indicator of effective fishing effort in mixed demersal fisheries[END_REF].

Creep has been estimated to be as much as 10% per year (i.e. [START_REF] Ellis | Effects of fish density distribution and effort distribution on catchability[END_REF][START_REF] Marchal | Impact of technological creep on fishing effort and fishing mortality, for a selection of European fleets[END_REF]). Meanwhile, many fish stocks are at a very low level in relation to the sustainability criteria. [START_REF] Worm | Rebuilding global fisheries[END_REF] estimated that 63% of assessed fish stocks worldwide are overexploited, and so, they require rebuilding. In the EU, 88% of assessed stocks are overexploited, 30% of these stocks being outside safe biological limits (EC 2009). The European fishing industry is experiencing smaller catches and facing an uncertain future. As a result, changes in fleet catch efficiency may reflect shifts in effort and or trends in fish abundance. But in comparison to the beginning of the study, the decrease observed lately in the fuel consumption (-33%) is proportional to the one observed in the landings (-34%). The poor situation of the fish stocks is probably prevailing over advances in technology.

A general tendency for poorer outcomes on the energy performance of the fishing fleets over time has been documented [START_REF] Tyedmers | Fisheries Impacts on North Atlantic Ecosystems: Catch, Effort and National/Regional Datasets[END_REF][START_REF] Tyedmers | Fisheries and energy use[END_REF]. [START_REF] Mitchell | Resource Scarcity, Energy Use and Environmental Impact: A Case study of the New Bedford, Massachusetts, USA, Fisheries[END_REF], for the New Bedford fishery, studied the evolution of the EROI ratio in 20 years (between 1968 and 1988) and concluded that it went from 0.18 down to 0.03. Instead, our results show that the selected EU fleets efficiency remained practically stable or even had a slight increase iv . [START_REF] Tyedmers | Fueling global fishing fleets[END_REF] estimated the energy content of the fuel burned by the fisheries globally to be 12.5 times greater than the edible energy content of the catch. The European fishing fleet seems to be relatively more fuel efficient than the global fisheries, with an estimated energy content of the fuel burned nine times greater than the edible energy content of the resulting catch v . The difference between the two estimations could be attributed either to different energetic performances of the fishing fleets, or a different composition of the fishing fleets represented in the datasets herein and worldwide. This is an interesting topic that needs to be further investigated; one could deduce that given the current overexploitation of most EU fish stocks (EC 2009) the energy performance of the fleets exploiting them should be worse than the global one. This can at least in part be explained by the high fuel prices in the EU that works as an incentive to improve the energy efficiency [START_REF] Cheilari | Effects of the Fuel Price Crisis on the Energy Efficiency and the Economic Performance of the European Union Fishing Fleets[END_REF]).

Although fisheries have been considered a fuel intensive food production process as stated by [START_REF] Wilson | Fuel and financial savings for operators of small fishing vessels[END_REF], our results confirm [START_REF] Tyedmers | Fisheries and energy use[END_REF] and [START_REF] Tyedmers | Fueling global fishing fleets[END_REF] conclusions that fisheries are not among the less fuel efficient food production activities. This conclusion does not change even if we consider a potential 25% overestimation on our results due to not accounting for the indirect energy uses (Pelletier et al. 2007). Comparing to other food production activities (table 3), fisheries EROI ratios tend to be lower. However, certain fleet segments (i.e., pelagic trawlers and seiners larger than 40 meters, pelagic trawlers and seiners between 24 and 40 meters and vessels smaller than 12 meters using passive gears), are more energy efficient than some agricultural processes. While compared to livestock and aquaculture production EU fishing fleets show a better energy performance, as they have on average a higher EROI, even if the results from this study are underestimated compared to studies that consider embodied energy.

Different fleets have different energy profiles. Pelagic trawls and seiners targeting pelagic fish can be among the most energy efficient production systems. Based on these results one could be tempted to incentivize pelagic trawls and seiners in detriment of less efficient fishing gears. However, not all fishing gears target the same species. Thus, comparison could only be feasible for those gears with same target species, such as demersal trawlers versus beam trawlers or pelagic trawlers versus pelagic purse seiners. Unfortunately, these data is not yet available in detail, but improvements in the data availability have been achieved during the last years and the European Union publishes data on fuel consumption and landings' weights per fleet and country (i.e. STECF 2010). Data by vessel could prove very important to determine what factors can help to improve energy efficiency (i.e. particular investments or fishing strategies), however, such studies require a more limited sample than the one used in this study. Currently data on landings and fuel consumption by vessel can only be available at the national and regional administrations. Similarly, data on material and energetic demands of fishing fleets (Pelletier et al. 2007) for such a large number of fleets is also lacking (except on direct energy use). More efforts should be allocated to collect these data at a larger scale, and not only on a case study basis, in order to be able to extract sound conclusions and recommendations for fisheries practices and management. Nevertheless, we should not ignore that fishing is an economic activity, and fishers decide to use a fishing technique based on economic incentives. Because fuel costs are an important part of the total costs, fishers are also interested in energy efficiency. But fuel costs are not the only driver to this decision. Other important factors, such as the relative price of potential target species, catchability between techniques (selectivity), other costs and previous investments in vessel and gears, may define the choice, at least in the short term and should be taken into consideration during the development of management strategies that can reduce the fuel consumption of the fishing industry.

The proposal for the financial tool of the EU's Common Fisheries Policy (CFP), the European Maritime and Fisheries Fund (EMFF), has made provisions for supporting measures that will increase energy efficiency, reduce emissions and contribute to the Europe 2020 target on climate change. However, it is not clear yet which specific innovations will be eligible for funding within the new EMFF. A non-exhaustive list could include areas such as: engines with reduced emissions (e.g. hybrid engines); bio-fuels; vessel design and technology; vessel operation (maintenance of hulls and engines); use of alternative/renewable energy sources (wind, H2 fuel cells); efficient fishing gears (e.g. reduced gear drag); fishing tactics and techniques; fuel management systems; energy audits. Energy audit can serve as a tool for assessing the performance of the fleets as well as the success of the innovative techniques applied [START_REF] Notti | An Energy Audit tool for increasing fishing efficiency[END_REF]. So far energy audits have proven to be quite effective in reducing consumption and increasing efficiency (Danish Ministry of Food, Agriculture and Fisheries 2011; [START_REF] Sala | Audit Energetici nel settore della pesca (e-Audit)[END_REF]. Other more hypothetical measures to improve energy efficiency not related with this fund could be the inclusion of fisheries in a carbon emission quota system (i.e. the European Union Emissions Trading System) or changes in the taxation of fuel for fisheries. Then, the EROI indicator could be a good measure to monitor energy efficiency performance for vessels and fleets (when targeting different species), as well as one of the criteria to allocate funds. i The term fish refers to all aquatic organisms harvested by the fisheries: fish, mollusks, crustaceans.

Tables

Supplementary information file

ii Includes data from Wiviot and Mathews 1975[START_REF] Bykov | Marine fishes: Chemical composition and processing properties[END_REF][START_REF] Research | Yield and nutritional value of the commercially more important fish species[END_REF][START_REF] Crapo | Recoveries and yields from Pacific fish and shellfish[END_REF] iii Boundary for energy inputs is direct energy and material inputs, while the boundary for energy outputs is extraction. iv A 2.6% annual increase of the overall EROI is statistically significant at a 10% level, but not at a 5% level. v EROI values of the EU fleet taking into account the indirect energy use, which refers to the cumulative energy required in upstream processes associated with the production and delivery of such inputs, are

  Active gears: DTS: Demersal trawl and demersal seiner, PTS: Pelagic trawls and seiners, TBB: Beam trawl, DRB: Dredges. Passive gears: HOK Gears using hooks, DFN: Drift nets and fixed nets, FPO: Pots and traps, PGP: Polyvalent passive gears, PG: Passive gears (applicable only to vessels with length lower than 12 m). Polyvalent gears: PMP: Combining mobile & passive gears.VL0012: vessels less than 12 m in length, VL1224: vessels between 12-24 m in length, VL2440: vessels between 24-40 m in length, VL40XX: vessels greater than 40 m in length.

FiguresFigure 1 .Figure 2 .Figure 3 .

 123 Figures

Table 1 .

 1 EROI values by fleet (main gear) for the period 2002-2008 (Source: own calculations from STECF 2010 data).

	Fleet (gear)	2002	2003	2004	2005	2006	2007	2008
	Drift nets and Fixed Nets	0.058	0.052	0.055	0.047	0.042 0.060 0.067
	Dredges	0.126	0.092	0.120	0.106	0.083 0.087 0.071
	Demersal Trawl and
	demersal Seiner	0.036	0.037	0.039	0.046	0.047 0.045 0.051
	Pots and traps	0.105	0.124	0.163	0.119	0.095 0.084 0.122
	Gears using hooks	0.013	0.014	0.023	0.019	0.025 0.034 0.041
	Passive gears	0.285	0.391	0.317	0.276	0.288 0.392 0.376
	Polyvalent passive gears	0.038	0.040	0.040	0.045	0.051 0.050 0.044
	Combining mobile &
	passive gears	0.078	0.082	0.081	0.205	0.154 0.134 0.111
	Pelagic Trawls and Seiners	0.534	0.417	0.461	0.541	0.560 0.536 0.617
	Beam Trawl	0.021	0.021	0.022	0.022	0.023 0.024 0.027
	Total Average	0.102	0.091	0.097	0.110	0.111 0.103 0.114

Table 2 .

 2 EROI values by fleet segment for the period 2002-2008 (Source: own calculations from STECF 2010 data).

	Fleet segment	2002	2003	2004	2005	2006	2007	2008
	DFN VL0012	0.030	0.032	0.035	0.026	0.028	0.042	0.048
	DFN VL1224	0.083	0.070	0.072	0.073	0.061	0.082	0.094
	DRB VL0012	0.411	0.359	0.340	0.205	0.143	0.178	0.201
	DRB VL1224	0.092	0.065	0.089	0.083	0.069	0.070	0.047
	DTS VL0012	0.022	0.035	0.045	0.047	0.045	0.041	0.049
	DTS VL1224	0.031	0.032	0.035	0.042	0.040	0.038	0.047
	DTS VL2440	0.042	0.043	0.045	0.053	0.059	0.056	0.056
	DTS VL40XX	0.078	0.075	0.061	0.061	0.088	0.083	0.081
	FPO VL0012	0.105	0.124	0.163	0.119	0.095	0.084	0.122
	HOK VL0012	0.013	0.014	0.023	0.019	0.025	0.034	0.041
	PG VL0012	0.285	0.391	0.317	0.276	0.288	0.392	0.376
	PGP VL0012	0.033	0.036	0.036	0.040	0.047	0.047	0.041
	PGP VL1224	0.129	0.111	0.135	0.135	0.120	0.151	0.150
	PMP VL0012	0.094	0.101	0.155	0.221	0.164	0.159	0.112
	PMP VL1224	0.076	0.078	0.069	0.203	0.152	0.128	0.111
	PTS VL1224	0.169	0.140	0.175	0.183	0.221	0.209	0.208
	PTS VL2440	0.533	0.419	0.428	0.562	0.496	0.459	0.470
	PTS VL40XX	0.941	0.813	0.800	0.912	1.029	1.024	1.125
	TBB VL1224	0.045	0.047	0.052	0.060	0.053	0.052	0.059
	TBB VL2440	0.016	0.014	0.016	0.014	0.013	0.016	0.020
	TBB VL40XX	0.016	0.015	0.015	0.014	0.017	0.017	0.016
	Total Average	0.102	0.091	0.097	0.110	0.111	0.103	0.114

Table 1

 1 Supplementary. EROI values by fleet segment for the period 2002-2008, taking into account the indirect energy use.

	Fleet segment	2002	2003	2004	2005	2006	2007	2008
	Drift nets and Fixed
	Nets	0.044	0.039	0.042	0.036	0.032	0.045	0.051
	DFN VL0012	0.023	0.024	0.027	0.020	0.021	0.032	0.036
	DFN VL1224	0.063	0.053	0.055	0.055	0.046	0.062	0.071
	Dredges	0.096	0.070	0.091	0.080	0.063	0.066	0.054
	DRB VL0012	0.312	0.272	0.258	0.155	0.108	0.135	0.152
	DRB VL1224	0.070	0.049	0.067	0.063	0.052	0.053	0.036
	Demersal Trawl and
	demersal Seiner	0.027	0.028	0.030	0.035	0.036	0.034	0.039
	DTS VL0012	0.017	0.027	0.034	0.036	0.034	0.031	0.037
	DTS VL1224	0.023	0.024	0.027	0.032	0.030	0.029	0.036
	DTS VL2440	0.032	0.033	0.034	0.040	0.045	0.042	0.042
	DTS VL40XX	0.059	0.057	0.046	0.046	0.067	0.063	0.061
	Pots and traps	0.080	0.094	0.124	0.090	0.072	0.064	0.092
	FPO VL0012	0.080	0.094	0.124	0.090	0.072	0.064	0.092
	Gears using hooks	0.010	0.011	0.017	0.014	0.019	0.026	0.031
	HOK VL0012	0.010	0.011	0.017	0.014	0.019	0.026	0.031
	Passive gears	0.216	0.296	0.240	0.209	0.218	0.297	0.285
	PG VL0012	0.216	0.296	0.240	0.209	0.218	0.297	0.285
	Polyvalent passive
	gears	0.029	0.030	0.030	0.034	0.039	0.038	0.033
	PGP VL0012	0.025	0.027	0.027	0.030	0.036	0.036	0.031
	PGP VL1224	0.098	0.084	0.102	0.102	0.091	0.114	0.114
	Combining mobile &
	passive gears	0.059	0.062	0.061	0.155	0.117	0.102	0.084
	PMP VL0012	0.071	0.077	0.117	0.168	0.124	0.121	0.085
	PMP VL1224	0.058	0.059	0.052	0.154	0.115	0.097	0.084
	Pelagic Trawls and
	Seiners	0.405	0.316	0.349	0.410	0.425	0.406	0.468
	PTS VL1224	0.128	0.106	0.133	0.139	0.168	0.158	0.158
	PTS VL2440	0.404	0.318	0.324	0.426	0.376	0.348	0.356
	PTS VL40XX	0.713	0.616	0.606	0.691	0.780	0.776	0.853