Do we need photonic-based instruments to characterize multiscale processes in marine environments ?

Do we need photonic-based instruments to characterize multiscale processes in marine environments ? Piera, J.1 Campbell, L.2 Sosik H.M.3 , Harred, L.2 ,Olson, R.J.3, Torrecilla, E.1, Ballabrera, J.1 1Marine Science Institute, Agencia Estatal Consejo Superior de Investigaciones Científicas, Barcelona, Spain. [email protected] 2Department of Oceanography, Texas A&M University, College Station, TX, USA 3Biology Department, Woods Hole Oceanographic Institution, Woods Hole, MA, USA Abstract Timeor space - series analysis is one of the main analytical methods in marine science to understand the dynamics of oceanic processes. It not only facilitates the identification of driving forces represented in sequences of observations, but also helps to forecast events by feeding models. When designing the sampling strategy for obtaining a time/space data series, one of the first questions to address is the sampling frequency and the number of samples to measure. In many cases the answers to these questions are based mainly on logistical or operational constraints (e.g., the maximum number of samples that can be processed, instrumental capabilities, etc.). However, it is important to take into account the principles of Information Theory in order to avoid potential artefacts derived from improper sampling design. Furthermore, in those cases where the processes are not stationary (a common situation in marine environments) the sampling frequency and the number of samples play an essential role in determining the time-frequency resolution required to fully characterize their dynamical properties. In this contribution, these problems of sampling requirements are illustrated using two example data sets: (a) a time series measured at a rate of 3 samples·h-1 by the Imaging FlowCytobot (IFCB), installed at the entrance to the Mission-Aransas estuary (Port Aransas, TX, USA) during 2008; and (b) a high resolution spatial series of fluorescence data obtained at a rate of 1 sample·m-1 with an autonomous underwater vehicle (AUV) in Alfacs Bay (Ebro Delta, Spain) during 2011. 1. Multi-scale analysis example: Time series of Dinophysis abundance The time series of Dinophysis abundance in the Port Aransas ship channel, measured by the Imaging FlowCytobot (Olson and Sosik, 2007; Sosik and Olson 2007) was used as a first example of a high resolution data set. Although the initial data was obtained at a rate of 3 samples·h-1, in order to reduce part of the variability derived from errors in the automatic identification, the data was averaged in 2 h integrated time intervals (i.e. 6 averaged samples for each interval). Further details of the sampling procedure can be found in Campbell et al (2010). The analysis was restricted to a relatively short time window, corresponding to the first 4 days of a bloom event (between February and March 2008, see Figure 1). This short period was chosen to ensure that the dynamical properties of the process remain constant and, may thus be considered stationary. Figure 1. Time series of Dinophysis abundance in the Port Aransas ship channel in early 2008. Automated analysis of IFCB data provided at a 2 h resolution (based on 6 averaged values sampled at a rate of 3 samples h-1) throughout most of the 4-month bloom. Under the assumption of stationarity, it is possible to model this selected interval on the basis of harmonic analysis, representing the time series as the linear combination of basic sinusoidal signals: x(nΔt)=A 0 + ∑ i=1 N A i cos(2πf i nΔt+ϕ i ) (1) Three different models were constructed (M1, M2 and M3) according to an increasing number of frequency components. M1 was constructed using 6 signal components in the low frequency range (N=6). M2 was based on the same frequency components of M1 with 6 additional components of higher frequency (N=12) and M3 was constructed based on M2 by adding 8 higher frequency components, thus using a total of 20 sinusoidal signals. Figure 2 shows the comparison of the model outputs with the original data. Figure 2. A1-C1: Time series reconstruction from harmonic analyses using the 3 proposed models (see text). Red crosses indicate the original observations. A2-C2 Comparison of the fitted results with the observations. Colour lines indicate the result of the linear regression estimation (black lines indicate the ideal result) in each plot the correlation coefficient is included. For practical purposes, we have considered that the errors in model 3 are negligible. We therefore use the associated mathematical expression of M3 as the full resolved data model. Although M1 and M2 do not completely resolve the dynamical properties of Dinophysis abundances, these models can also be considered valid approaches depending on the time scale we want to consider for further analysis. As we can see in Figure 3 (A-B), these models may provide valid information to predict averaged data values of Dinophysis abundances over different integrating time intervals. The analogy between data averaged at different scales in the time domain and low-pass filtering in the frequency domain (obtaining a partially resolved spectrum) is shown in Figure 3D. This simple analysis provides an example of a multiscale approach to characterize the dynamics of Dinophysis during the considered time frame of 4 days. It is a graphical demonstration that depending on the time scale we want to analyse, we do not need to resolve the complete spectral properties (M3) to evaluate the average properties in a particular integrating time interval. Figure 3. A-C. Prediction of averaged values at different times scales: (A) Model 3 is able to reproduce the dynamics at 2 h integrating time interval. Model 2 (B) and Model 1 (C) also reproduce the averaged values at 16 h and 24 h integrated time intervals respectively. (D) represents the normalized amplitude of the different frequency components to build models 1 to 3. 2. The requirements for selecting the optimal sampling frequency Obviously, depending on the integrating time interval, the number of averaged samples will change. The initial 48 samples at 2 h integrating intervals (Model 3, Fig 3A) become only four averaged samples at 24 h integration time (Model 1, Fig 3 C). In the frequency domain, there is also an associated reduction from the 20 frequency components (and maximum frequency of 0.204 h-1) in M3, to 6 components (and maximum frequency of 0.061 h-1) in M1. It is important to take into account that the characterization at different time scales (and the implicit reduction of the averaged time intervals and frequency components for larger time scales) does not imply the possibility to reduce the sampling frequency when dealing with analyses over larger scales. As an example, Figure 4 indicates the potential errors trying to characterize the daily dynamics of the Dinophysis abundance based only on a single measurement (of the original time series) in each 24 h interval. Different time series can be obtained, most of them with significant differences between the “unique” daily value (in cross symbols) and the averaged reference, (grey lines). Notice differences among the trends of global changes of the three time series considered, when single daily values are obtained at different intervals of the day. Figure 4. Examples of the potential artefacts derived from an improper sampling frequency (1 estimated sample in 24 h) to characterize the dynamics of the Dinophysis abundance on a daily basis. The grey lines correspond to the averaged value of all the experimental data within each integrating interval. Blue points correspond to the averaged values obtained during the four days at 9-11 h interval. Red and green correspond to the averaged intervals 11-13 h and 13-15 h respectively. The example shown in Figure 4 illustrates the importance to set the correct sampling frequency for a particular time scale. The results indicate that it is not possible to infer the daily dynamics of Dinophysis abundance based only on a single measurement per day. So how many samples per day do we need to address this problem? The theoretical framework to answer this question was already defined nearly a century ago by several authors (for historical review see Butzer and Stens (1992)) providing the known Sampling Theorem, which set the foundation of Information Theory. In essence, the theorem shows that a band-limited continuous signal can be reconstructed from an infinite sequence of discrete samples if the sampling frequency exceeds 2fmax samples per second, where fmax is the highest frequency of the original signal. Perfect reconstruction is mathematically possible for the idealized model (infinite number of samples) but only an approximation for realworld signals and sampling techniques, although in practice it is often a very good one. The constructive proof of the Sampling Theorem leads to an understanding of the aliasing that can occur when a sampling system does not satisfy the conditions of the theorem. A numerical experiment was designed to illustrate the aliasing effects on the example of Dinophysis abundance time series. The numerical experiment used the mathematical expression of M3 as the fully resolved data spectrum. Four different numerical experiments were conducted to create the corresponding time series. In these simulations, M3 was re-sampled at different sampling rates (and duration of the observation) according to Table 1. Ref Number of samples Sampling interval (frequency) Duration of observation E1 4096 2 h (0.5 h-1) 341 days E2 4096 12 h (0.083 h-1) 5.6 years E3 4096 24 h (0.042 h-1) 11.2 years E4 4096 48 h (0.021 h-1) 22.4 years The number of samples was chosen very high (212, 4096) in order to avoid effects of short time series. The simulations provided time series that, in the context of biological observations, can be considered as “very-high quality” in terms of temporal resolution and coverage. It would be difficult, for example, to find time series of biological field observations that cover more than 5 years with samples obtained regularly every 12 h, as it is the case of the simulation E2. However, only the first experiment E1 satisfies the conditions of the Sampling Theorem: the maximum frequency component in M3 is 0.204 h-1. According to the theorem, the minimum sampling rate should be twice this maximum frequency (0.408 h-1). In E1 the sampling rate is 0.5 h-1, but in the rest of the simulations, the sampling frequencies are far below the required threshold: 0.083 h-1, 0.042 h-1 and 0.021 h-1 respectively. The aliasing effects of not satisfying the conditions of the Sampling Theorem are illustrated in Figure 5 with the respective plots of the times series in the frequency domain. Figure 5. Aliasing effect of resampling model 3 with sampling periods. The spectra were computed using the FFT transform of the simulated time series. The first spectra (E1) is the only one resolved correctly (without aliasing) and has been used as a reference to indicate the elevated number of frequency artefacts in the range of low frequency processes (areas in grey) that are produced in the rest of the observations: E2 sampling interval 12 h, E3 sampling interval 24 h, E4 sampling interval 48 h. Discussion The aliasing effect sometimes is referred as “mirror effect” in the frequency domain: The false frequencies will appear as mirror images of the original frequencies around one half of the sampling frequency. According to the aliasing-mirror effect it is possible to classify the sampling in three general scenarios (Fig. 6): Figure 6. Classification of sampling scenarios according to the aliasing-mirror effect (see text). (A) No aliasing. (B) Partial aliasing. Grey bars indicate original frequencies and red bars aliased components. (C) Complete aliasing. Colour codes are the same in B and the green bars correspond to double mirrored aliasing components. (D) Detail of Partial aliasing indicating the spectral region that is not contaminated (0 to falias-free ) 1. No aliasing. When the sampling frequency is above 2 fmax and there is no mirror effect. In this case the spectrum is fully resolved. (Fig 6A) 2. Partial aliasing. When the sampling frequency is in the range of ( 2 fmax , fmax ). In this case part of the high frequencies are mirrored (Fig 6B). The spectrum is partially contaminated by the aliased frequencies but it is possible to recover the information corresponding to the low frequency range. 3. Complete aliasing. When the sampling frequency is below fmax. In this case there are multiple reflections of the mirror effect and the spectrum is completely contaminated (Fig 6C). In this case it is very complicated (and some times impossible) to recover correctly the dynamic properties in any temporal scale. Notice that the results of figure 5 indicate that the simulations E2 to E4 are examples of complete aliasing. The requirements to satisfy completely the conditions of the Sampling Theorem imposes serious challenges for environmental observations. However, in some cases it is possible to design a sampling strategy to obtain a time series with partial aliasing, retrieving only one fraction of the whole spectrum corresponding to the low frequency components. Nevertheless, we have seen in section 1 that resolving only part of the spectrum (models M1 and M2, fig 3B-C) can be also useful to predict dynamical properties at particular scales. For a partial aliasing strategy we need to define the spectral region of interest, identifying the maximum frequency without aliasing (falias-free in Fig 6D). We need to know which is the maximum frequency of the process (fmax ) as well. The sampling frequency (fs) can be estimated (according Fig 6D) using the following equations: falias-free > 0.5fs - ∆f (2) ∆f = fmax - 0.5fs (3) and combining (2) and (3) we obtain the general condition to estimate fs fs > falias-free + fmax (4) Notice that in the case that we want to resolve the complete spectrum (in this case falias-free is equal to fmax ), equation (4) becomes fs > 2 fmax (5) which is the general expression of the Sampling Theorem. Applying the equation (4) to the results obtained in the analysis of the Dinophysis abundance, it is possible to answer the issue raised in the previous section: which should be the required sampling frequency for a daily basis characterization? The falias-free can be estimated as the inverse of 24 h period (0.042 h-1) and the maximum frequency, fmax according to the model M3, is 0.204 h-1. The result (0.246 h-1) provides the low frequency sampling threshold, that corresponds approximately to a minimum of 6 samples per day. 4. Identifying (and minimizing) the maximum frequency component With the results obtained in the previous sections it is easy to conclude that in those systems with fast dynamics (i.e. high fmax ) the frequencies of interest (falias-free) play a minor role in determining the sampling requirements (equation 4). In these cases it is critical to determine the maximum frequency it is possible for the sampling system to measure. In many conventional sampling devices (where a sensor converts the parameter of interest in a continuous electrical signal that can be digitized further) the aliasing problem is solved by using an analog filter circuit (known as anti-aliasing filter) that blocks all the problematic high frequency components previously to the sampling process. The cut-off frequency of the filter is designed to satisfy the conditions of the Sampling Theorem. In environmental observations it is commonly very difficult to apply this solution for two main reasons: a) in many cases the process of sensing the parameter of interest does not involve any conversion to a continuous electrical signal, so it is not possible to apply an analog filter circuit to block the high frequency components. b) in those cases where electrical conversion could be possible (for example, in the case of a fluorometer), the cut-off frequency of the analog filter required in many of the conventional sampling scales (daily, weekly or monthly basis) would be so small that there are no standard electronic components (resistors and capacitors) available to implement the design of that filter. In these cases it is necessary to identify the maximum frequency that the observing system is able to measure and trying to adapt the sampling frequency accordingly. It is very difficult to estimate a precise value for fmax, as it depends not only on the dynamics of the system but also on the dynamic response of the instruments and methods we are using for sensing the parameters of interest. As a preliminary evaluation, we have considered in this discussion the sources of variability that may increase the maximum frequency of our observations. Basically there are two major components of variability in any (Eulerian) observation: 1) The local processes (in our example they will include processes such as cell division, grazing, …). Usually the local processes have associated a relative large time scales of variability (the order of hours or days) 2) The transport processes (advection and turbulence). The time scales of variability associated to these processes will depend on the velocity fields and the spatial distribution of the parameter under study. The higher the spatial gradients, the higher variability and the maximum frequency. As an example of spatial gradients, figure 7 shows the distribution at small scale of different parameters (Temperature, Salinity and Fluorescence as a proxy of total chlorophyll a) that were measured with an Autonomous Underwater Vehicle in Alfacs Bay (Spain). The results indicate strong gradients over short distances. These results, rather than being an exception, seems to be the general situation. Many field studies (among them Cowles et al., 1993; Seuront and Menu, 2006; Doubell et al., 2006; Mitchel et al. 2008) support the idea that high spatial biological gradients, that would yield in a high frequency variability in any Eulerian observation, are the normal case in coastal and oceanic waters.