Time series analysis of water surface temperature and heat flux components in the Itumbiara Reservoir ( GO ) , Brazil

Aim: Water temperature plays an important role in ecological functioning and in controlling the biogeochemical processes of the aquatic system. Conventional water quality monitoring is expensive and time consuming. It is particularly challenging for large water bodies. Conversely, remote sensing can be considered a powerful tool to assess important properties of aquatic systems because it provides synoptic and frequent data acquisition over large areas. The objective of this study was to analyze time series of surface water temperature and heat flux to advance the understanding of temporal variations in a hydroelectric reservoir. Method: MODIS water-surface temperature (WST) level 2, 1 km nominal resolution data (MOD11L2, version 5) were used. All available clear-sky MODIS/Terra images from 2003 to 2008 were used, resulting in a total of 786 daytime and 473 nighttime images. Time series of surface water temperature was obtained computing the monthly mean in a 3 × 3 window of three reservoir selected sites: 1) near the dam, 2) at the centre of the reservoir and 3) in the confluence of the rivers. In-situ meteorological data from 2003 to 2008 were used to calculate surface energy budget time series. Cross-wavelet, coherence and phase analysis were carried out to compute the correlation between daytime and nighttime surface water temperatures and the computed heat fluxes. Results: The monthly mean of the day-time WST shows lager variability than the night-time WST. All time series (daytime and nighttime) have a cyclical pattern, passing for a minimum (June July) and a maximum (December and January). Fourier and the Wavelet Analysis were applied to analyze this cyclical pattern. The daytime time series, presents peaks in 4.5, 6 12 and 36 months and the nighttime WST shows the highest spectral density at 12, 6, 3 and 2 months. The multiple regression analysis shows that for daytime WST, the heat flux terms explain 89% of the annual variation (RMS = 0.89 °C, p < 0.0013). For nighttime, the heat flux terms explain 94% (RMS = 0.53 °C, p < 0.0002). Conclusion: The daytime WST and shortwave radiation presents a good agreement for periods of 6 (with shortwave retarded) and 12 months (with shortwave advanced); For nighttime WST and longwave the good agreement is present for 1, 3, 6 and 12 months, all with longwave advanced in relation to WST.

Moreover, temperature differences between the water and air control the heat exchange in the air/ water boundary layer, and as a consequence, they are crucial to understanding the hydrological cycle.A study of the energy exchange between the lake and atmosphere is essential for understanding the aquatic system behavior and its reaction to possible changes of environmental and climatic conditions (Bonnet et al, 2000).The detection of trends and sudden changes in the aquatic system is dependent on both the availability of long-term time series of environmental data and their proper analysis (Stech et al., 2006).
Thermal infrared remote sensing applied to freshwater ecosystems has aimed to map surface temperatures (Crosman and Horel, 2009;Alcântara et al. 2010a), bulk temperatures (Thiemann and Schiller, 2003), circulation patterns (Schladow et al., 2004) and to characterize upwelling events (Steissberg et al., 2005).However, the application of thermal infrared images to the
Most of these satellites acquire data twice every 16 days, such as Landsat and ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer), at any given location.
The Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the Terra and Aqua satellites overcome these temporal resolution constraints, since MODIS data can typically be acquired daily due to their large scan angle (Justice et al., 1998).
Looking at reservoirs as thermodynamic systems, they can be approached as thermal engines which exchange heat with their surroundings; their functioning can be characterized by the development of their heat content, heat budgets, and water temperature.The solar energy and the thermal radiation reaching the water surface are quickly transformed into heat in the upper layers of the water column.At the same time, the wind-driven turbulent kinetic energy spreads this absorbed heat in all three dimensions of the water body (Moreno-Ostos et al., 2008).

Study area
The Itumbiara hydroelectric reservoir (18° 25' S, 49° 06' W), classified as case 2 water type (Jerlov, 1968;Nascimento et al., 2011), is located in a region stretched between Minas Gerais and Goiás States (Central Brazil) that was originally covered by tropical grassland savanna.The damming of the Parnaiba River flooded its main tributaries: the Araguari and Corumbá rivers.The basin's geomorphology resulted in a lake with a dendritic pattern covering an area of approximately 814 km 2 and a volume of 17.03 billion m 3 (Figure 1).
The climate in the region is characterized by an average precipitation ranging from 2.0 mm in the dry season (May -September) to 315 mm in the rainy season (October -April).In the rainy season the average wind intensity ranges from 1.6 to 2.0 m/s and reaches up to 3.0 m/s in the dry season; the preferential wind direction is southeast.The air temperature in the rainy season (October -April) study of surface water temperatures in hydroelectric reservoirs is scarce and, in Brazil, is being attempted for the first time.
Several satellites have been launched with spatial, temporal and radiometric resolutions for the study of surface water temperatures with relative accuracy (Steissberg et al., 2005).However, most of these satellites acquire data twice every 16 days, such as Landsat and ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer), at any given location.Also, Landsat does not record data at night (unless by special request).The Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the Terra and Aqua satellites overcome these temporal resolution limitations.MODIS data can typically be acquired daily due to their large scan angle (Justice et al., 1998).
Based on this, the aim of this research is to analyze the trends in water surface temperature and their periodical relationship with the heat flux using thermal remote sensing time series.ranges from 25 to 26.5 °C, decreasing to 21 °C by the beginning of the dry season (June).The relative humidity has a pattern similar to that of the air temperature, but with a small shift in the minimum value towards September (47%).Moreover, during the rainy season the humidity can reach 80% (Alcântara et al., 2010a).

Hydrometeorological data
Four hydrometeorological variables were measured from 2003 to 2008 period: the daily mean air temperature (°C), relative humidity (%), wind intensity (m/s) and precipitation (mm).These data were obtained from a meteorological station (see Figure 1 for location) near the dam.The daily mean of each variable was converted into monthly means to adequate it to the time scale of the satellite data.

Satellite data
MODIS water surface temperature (WST) level 2, 1 km nominal resolution data (MOD11L2, version 5) were obtained from the National Aeronautics and Space Administration Land Processes Distributed Active Archive Center (Wan, 2008).All available clear-sky MODIS Terra imagery between 2003 and 2008 were selected by visual inspection, resulting in a total of 786 daytime images and 473 nighttime images (Figure 2).
The WST-MODIS data were extensively validated for inland waters and were considered accurate (Oesch et al., 2005(Oesch et al., ,2008;;Reinart and Reinhold, 2008;Crosman and Horel, 2009).A shoreline mask to isolate land from water was built using the TM/Landsat-5 image in order to isolate some anomalously cold or warm pixels remaining at some locations near the shoreline of the reservoir (Sentlinger et al., 2008).
The mask was built based on the water level fluctuations in the reservoir; that is, if the MODIS data is from January, then a TM/Landsat-5 was processed take into account the lake level.Each masked MODIS image was checked to assure that each pixel in the border is removed.The mask was built using the NDWI (Normalized Difference Water Index) algorithm proposed by McFeeters (1996) as (Equation 1): where NIR is the near infrared spectral band.where φ sf is the sensible heat flux (Wm -2 ), r a = 1.2 kgm -3 is the air density, c p = 1.005×10 3 Jkg -1 K -1 is the specific heat capacity of air, c H = 1.1×10 -3 is the coefficient of turbulent exchange and V  = surface wind speed (m/s).
The latent heat flux was computed as follows (Large et al., 1997) (Equation 7): where φ lf is the latent heat flux (Wm -2 ), c E = 1.1×10 -3 is a coefficient of turbulent exchange, L = 2.501×10 6 Jkg -1 is the vaporization of latent heat and p a is the atmospheric surface pressure (mb).
The energy exchange also occurs through precipitation, chemical and biological reactions in the water body, and conversion of kinetic to thermal energy.These energy terms are small enough that should be omitted.Many researchers agree that omitting energy budget components with small values does not significantly affect the results (Bolsenga, 1975;Sturrock et al., 1992;Winter et al., 2003).The sensible and latent heat flux was calculated for daytime and nighttime using the monthly mean surface water temperature derived from MODIS data.
Due to the complexity of these fluxes and the limitations of the atmospheric data available for the area of study, some constraints were imposed.The air temperature (T a ) and wind intensity ( V  ) were considered the same for the whole reservoir because only one meteorological station was available.Also, the incident solar radiation (φ s ) was calculated through Equation 2 while considering that the calculated incoming shortwave radiation was the same for the whole reservoir.All heat flux components were computed previously and can be accessed in detail in Alcântara et al. (2010a).

Correlation between the WST and the heat fluxes
Correlation and multiple regression analyses were carried out to better understand the relationship between the surface water temperature (daytime and nighttime) variability and heat flux terms.

Time series
We used a 3x3 moving window to smooth original maps and estimate mean values, minimizing noises on output maps and to generate time series in three different reservoir regions: (1) dam, p1, (2) central area, p2, and (3) confluence of the main rivers, p3 (See Figure 1 for location).

Surface energy budget
The surface energy budget for the Itumbiara reservoir was estimated using lake surface temperature data from the MODIS-WST and surface air temperature, wind speed, and relative humidity data from the hydro-meteorological surface meteorological station.
The exchange of heat across the water surface was computed using the methodology described by Henderson-Sellers (1986) as (Equation 2): where φ N is the surface heat flux balance, φ s is the incident short-wave radiation, A is the albedo of water (=0.07,MacIntyre et al., 2002;Laird and Kristovich, 2002), φ ri is the Longwave flux, φ sf is the sensible heat flux and φ lf is the latent heat flux.The units used for the terms in Equation 2are Wm -2 .
The incident short-wave radiation φ s was calculated using the following Equation 3: where a 1 = 0.79 and b 1 = 1.15 are two calibration parameters determined by a comparison with the radiometer data, φ 0 = 1390 Wm -2 is the solar constant, d is the hour angle and C is the cloud cover index, which was estimated using the Reed (1977) empirical relation.
The net longwave radiation corresponding to the outgoing flux minus the incoming flux (LW ↑↓ = LW ↑ -LW ↓), φ ri is computed as given by Fung et al. (1984) (Equation 4): where ε = 0.97 is the thermal infrared emissivity of the water, σ is the Stefan-Boltzmann constant, T s = water surface temperature (K), T a = surface air temperature (K), λ = 0.8 is the Reed (1977) correction factor, and e a = partial pressure of vapor (mb), which was calculated as (Equation 5): where r is the relative humidity and e sat is the saturation vapor pressure.This was calculated using the polynomial approximation of Lowe (1977).
The non-radiative energy term φ rad accounts for the sensible heat flux and latent heat flux.The sensible heat flux was estimated as (Large et al., 1997) (Equation 6): The time frequency space of WST time series was analyzed using the wavelet transform.Wavelet analysis is becoming a common method for analyzing localized variations of power within an environmental time series (Meyers et al., 1993;Kumar and Foufoula-Georgiou, 1997;Massei et al., 2006).Decomposing an environmental time series into time-frequency space allows for the determination of the dominant modes of variability and their modes of variation in time.
Continuous wavelet analysis was applied to WST time series from satellite data using the Morlet wavelet as reference function (the so-called "mother wavelet").The Morlet wavelet is the most common wavelet transform and consists of a Gaussianwindowed complex sinusoid defined in the time and frequency domains by (Equation 13): Where ω 0 is the non-dimensional frequency, here assumed to be 6 to satisfy the admissibility condition (Farge, 1992); η is a non-dimensional time parameter and ψ 0 (η) is the wavelet function.
Wavelet spectral power at different scale (ω) and time location (τ) can be calculated by (Equation 14): Where W is the wavelet transform described below.
The continuous wavelet transform of a discrete sequence x n is defined as the convolution of x n with a scaled and translated version of ψ 0 (η) (Torrence and Compo, 1998).The wavelet transform W for a given time series x is calculated by (Equation 15): Where *indicates the complex conjugate; δt is the uniform time step, k' is an integer from 1 to N (number of data points) and x n is the time series data.The scale-averaged spectral-power-based wavelet analysis reflects the average variance for different time scales (frequency or period).The calculation procedures of continuous wavelet analysis were coded in Matlab 6.5 (The MathWorks, Inc., Natick, MA).

Cross wavelet, coherence and phase
To show the relationship between WST and the time series identified as statistical significant by the 2.6.Time series analysis

Fourier power spectrum
The WST time series was analyzed using the Fourier power spectrum.The function that implements the transform is given by (Equation 8): , N is the number of points and X is the discrete Fast Fourier Transform of the WST time series.
Data windowing was carried out to avoid leakage in the estimation of the periodogram.Data windowing modifies the relationship between the spectral estimate at discrete frequencies and its continuous (periodogram) spectrum at nearby frequencies.There are many window functions used to prevent spectral leakage, most of them rising from zero to a peak and then falling again.This paper uses Hanning and Hamming windows (Press et al., 1992).The first was used to prevent leakage while the last was used to smooth the resulting spectra.The Hanning and Hamming window coefficients are show in (Equation 9 and 10), respectively: Power spectrum smoothing reduces the variance and increases the statistical confidence, or reduces the confidence limits.In the algorithm implementation there must be a compromise between strong (more confidence but stronger bias) and weak smoothing (less confidence but less bias).In this work smoothing was carried out with a variable length Hamming window.The inferior and superior confidence limits for the spectra area, respectively, are given by (Equation 11 and 12): where χ is the Chi square distribution, df is the number of degrees of freedom, a = 1 -p / 100 and p is the confidence level (95%).
Also is easy to observe in all time series (daytime and nighttime) a cyclical pattern, passing for a minimum (June -July) and a maximum (December and January).This cyclical pattern can be analyzed through the Fourier and the Wavelet Analysis.

Fourier power spectrum
The Figure 4 shows the Fourier Power Spectrum for both daytime or nighttime WST.The daytime time series, p1 (Figure 4a) and p2 (Figure 4b) presents two peaks with same periods, that is the annual cycle and the semiannual cycle.The period of 4.5 months appear only in the p1 time series.The maximum density for the time series p1 and p2 occurs for periods of 6 and 12 months, respectively.The time series p3 (Figure 4c) shows a high density for periods of 6 months, followed by 12 and 36 months.
The analysis of nighttime WST for p1 (Figure 4d), p2 (Figure 4e) and p3 (Figure 4f ) shows that the highest spectral density were 12, 6, 3 and 2 months.As was emphasized by Kaiser (1994) the Fourier transform can hidden periods with significant spectral density, by imposing a scale or response interval dependent of the time series size (in our case 72 months).Based on this Kaiser suggested a time-frequency analysis using the Wavelet Power Spectrum.

Wavelet analysis
The time series obtained from the MODIS/Terra images shows that the principal range of variability is from 8 to 24 months.The time series p1 (Figure 5a) shows that for the year 2003 (from June to July) the maximum peak of energy occurs for periods ranging from 8 to 12 months.From year 2004 (months of May, September and December) to 2005 (March and April) the higher energy peak occurs for periods bellow 8 months (probably 6 months, as showed by the Fourier power spectrum).
In 2006 the energy peaks occurs in May, June, November and December with periods until 24 months; whereas in 2007 the peaks occurs in May, June and July with periods of significantly variability up to 12 months.In 2008 only a small peak of energy, if compared with others years, was observed during May for periods smaller than 12 months.
The time series p2 (Figure 5b) shows that highest energy peaks not exceed periods of 12 months in all years of analysis.The months which showed these highest or maximum spectral density were the same Pearson's correlation coefficient, the cross wavelet analysis was carried out and the coherence and phase analyzed (Grinsted et al, 2004).
The cross wavelet transform of two time series x n and y n is defined as W xy = W x W y* , where *denotes complex conjugation, so the cross wavelet power could be defined as |W xy | (Grinsted et al. 2004).The interpretation of complex argument arg (W) can be interpreted as the local relative phase between x n and y n in time frequency space.The cross wavelet power theoretical distribution of two time series with background power spectra P x k and P y k is given in Torrence and Compo (1998) as (Equation 16): where Z v (p) is the confidence level associated with the probability p for a probability density function (pdf ) defined by the square root of the product of two χ 2 distributions.Cross wavelet power reveals areas with high common power.
Another useful measure is how coherent the cross wavelet is in time frequency space.The wavelet coherence could be defined following Torrence and Webster (1999) as (Equation 17): where S is a smoothing operator; Notice that this definition closely resembles that of a traditional correlation coefficient, and is useful to think of the wavelet coherence as a localized correlation coefficient in time frequency space.The calculation procedures of cross wavelet and wavelet coherence were coded in Matlab 6.5 (The MathWorks, Inc., Natick, MA).

Water surface temperature (WST) time series
The time series of WST measured by satellite are show in Figure 3.The WST measured during daytime and nighttime in Figure 3a,b, respectively.
The monthly daytime mean WST shows larger variability than the nighttime WST.According to Chapra (1997) the parameters that can explain the WST are shortwave radiation, wind speed and duration.During daytime it is possible to observe a small distinction between the three time series (p1, p2 and p3) which is not clear in the nighttime WST. that was previously presented for the time series p1 (Figure 5a).
For the time series p3 (Figure 5c) there is energy peaks up to 60 months (five years) occurring in 2004 for months of September and December.Also as was verified before in time series p1 and p2, energy peaks of 24 months occur in December 2006.We also found peaks for periods of 60 months in May 2006 and February 2007.
For nighttime WST, the time series are show in Figure 5d-f where can be observed energy peaks of up to 24 months.Nevertheless, these peaks are

Statistical model for daytime and nighttime water surface temperatures
The Pearson correlation computed between the daytime and nighttime WST derived from the MODIS image against the heat flux terms is shown in Table 1.For daytime temperatures, the only significant variable correlated heat flux was the shortwave radiation; for nighttime temperatures, it was the longwave radiation, latent heat flux and sensible heat flux.
The multiple regression analysis shows that for daytime WST, the correlated heat flux terms explain 89% of its annual variation (RMS = 0.89 °C, p < 0.0013).For nighttime, the heat flux terms explain 94% (RMS = 0.53 °C, p < 0.0002).The representative Equations of these relationships are presented in Equations 18 and 19: best defined and with smaller energy than those of daytime time series.Also it is a fact that, in general way, the months with maximum density energy includes the periods from May to August and from November to April of the subsequent year.In these cases a little delay were recorded in some samples.
The time series of nighttime WST of highest energy was the p1 (Figure 5d) followed by the time series p2 (Figure 5e) and with a smallest spectral density the p3 (Figure 5f ).In the case of p3 the highest energy was observed for small periods, less than 16 months and small energy for periods higher than 20 months.
An analysis of correlation between the WST and the heat flux components was computed to check out the relative importance of each heat flux in the WST time series.when the shortwave radiation increases, the air temperature increases and transfers heat in the water column.However, the pattern of wind intensity and direction during the day is also important; with high wind intensity, the surface water loses heat through convection.
At night, the processes of convection and advection acting at the surface water and the balance of longwave radiation are important to drive the water temperature.Because of this, Equation 18is a subtraction of the joint effects of these three important fluxes.where T daytime is the daytime water surface temperature (°C), T daytime is the nighttime water surface temperature (°C), φ S is the short wave radiation (Wm -2 ), φ ri is the long wave radiation (Wm -2 ), φ lf is the latent flux (Wm -2 ) and φ Sf is the sensible flux (Wm -2 ).
Only the incoming shortwave radiation was needed to model the daytime WST.This is because  The cross wavelet analysis between the nighttime WST and φ sf is show in the Figure 8. Four high common power were identified, (1) band period of 2-3 months localized from September to December 2004 (Figure 8a), with φ sf retarded 45° [~9 days] in The cross wavelet between the daytime WST and the φ S shows two band period of high common power, the first in the band period of 4.5-7 (regions 1 and 3) months, in special from years of 2004 and 2007; the second with the highest common power between 9-15 months (region 2) also with special attention in the years of 2004 and 2007 (Figure 6a).
The coherence in those common power periods (Figure 6b) shows that for the region 1 (period of 4-7 months) the φ S is retarded 45° [~21 days] in relation to WST (from January 2004 to July 2005).That is, for a period from 4-7 months when the φ S presents a peak of variability the WST needs 21 days to present like a response to this solar enhancement.The El-Niño was present during 2004-2005 years and can affect the regional climate, principally the air moisture and the rainfall regime.
For the region 2 (period of 8-16.5 months) the φ S is advanced 45° [~1.5 months] in relation to WST.In this case, for periods higher than 8 months the WST tends to reach their maximum variability but the φ S will present its higher variability 1.5 months later.The region 3 (band period of 4.5-6.5) is representative of the 2007 year (La-Niña year) with the φ S is retarded 90° [~1.25 months] in relation to WST.

Nighttime water surface temperature versus longwave radiation (φ ri )
The cross wavelet analysis between the nighttime WST and φ ri is show in the Figure 7. Five high common power was identified: (1) band period of 1-2 months localized between November 2004 and April 2005 (Figure 7a), with φ ri advanced 45° [~6 days] in relation to nighttime WST (Figure 7b); (2) band period of 3-4 months localized between The WST variability is not only due to heat flux (shortwave, longwave, sensible and latent), but also the physical processes that occur into the water column (Alcântara et al., 2011b).According to Alcântara et al. (2010b)  The cross wavelet analysis between the nighttime WST and φ lf is show in the Figure 9 and were identified four high common power, (1) band period of 1-3 months localized from August 2004  modifying the vertical temperature profile and consequently the temperature distribution in the water surface.

Conclusions
The results obtained allow reaching the following conclusions: • If we assume that the water surface temperature (WST) could be explained only the heat fluxes then for daytime WST only the shortwave radiation explain 89% of the temperature variability; to nighttime WST the longwave, sensible and latent flux explain 94%; • The daytime WST and shortwave radiation  , vol. 23, no. 3, p. 245-259

Figure 1 .
Figure 1.The Itumbiara hydroelectric reservoir and their bathymetry shows the location of SIMA buoy, meteorological station and the samples of water surface temperature time series (p1, p2 and p3).

Figure 2 .
Figure 2. Information about the MODIS thermal infrared images: (a) hour of passage and number of satellite imagery available for daytime (b) and nighttime (c).

Figure 3 .
Figure 3.Time series (3x3) of monthly mean water surface temperature (WST) in three selected stations (p1, p2 and p3) from 2003 to 2008: (a) daytime WST, (b) nighttime WST.Where p1 was taken near the dam, p2 in central region and p3 in the river confluence (see Figure 3 for location).

Figure 5 .
Figure 5. Wavelet analysis of the WST time series: (a) daytime WST time series p1, (b) daytime WST time series p2, (c) daytime WST time series p3, (d) nighttime WST time series p1, (e) nighttime WST time series p2 and (f ) nighttime WST time series p3.Cross-hatched regions on either end indicate the "cone of influence," where edge effects become important.

Table 1 .Figure 6 .
Figure 6.Cross wavelet between daytime WST and φ S (a) and the coherence and phase (b).
the passage of cold front over the reservoir can modify the heat flux and consequently the WST.The cold front enhances the convective processes in the water column, relation to nighttime WST (Figure 8b); (2) band period of 5-7 months localized from April 2004 to April 2005, with φ sf advanced 135° [~2 months]; (3) band period of 5.5-7 months between May 2007 to March 2008, with φ sf advanced 135° [~2 months]; and (4) band period of 9-15 months localized between February 2004 until November 2007, with φ ri advanced 45° [~2 months].3.5.4.Nighttime water surface temperature versus latent flux (φ lf )

Figure 8 .Figure 7 .
Figure 8. Cross wavelet between nighttime WST and φ sf (a) and the coherence and phase (b).
presents a good agreement for periods of 6 (with shortwave retarded) and 12 months (with shortwave advanced); • For nighttime WST and longwave the good agreement is present for 1, 3, 6 and 12 months, all with longwave advanced in relation to WST; • The nighttime WST and sensible flux high common power for band periods of 2 (retarded), 6 (advanced) and 12 (advanced); • Finally, the nighttime WST and latent flux with band periods of 2 (retarded), 6 (retarded) and 12 months (WST and latent flux in antiphase).

Figure 9 .
Figure 9. Cross wavelet between nighttime WST and φ lf (a) and the coherence and phase (b).