Spatial Interpolation of Two-Wavelengths Bio-Optical Models to Estimate the Concentration of Chlorophyll-a in A Tropical Aquatic System

Bio-optical models have been used to estimate and map the concentration of chlorophyll-a (chl-a) in aquatic systems. Bio-optical models’ algorithms try to infer the concentration of optically active components in water from their inherent optical properties (IOPs). We proposed the use of two single wavelengths to retrieve chl-a concentration in a tropical aquatic system. The results were compared to in situ measurements of chl-a. To spatialized the results of the bio optical modeling, we tested the spatial interpolation following two methods: (1) ordinary kriging technique was used to spatialize the calculated values of estimated chl-a for each model; (2) ordinary kriging was used to spatialize the estimated Rrs from 470nm and 700nm then the two wavelengths models were calculated by a map algebra using the spatialized Rrs. We generated four different spatial interpolations of chl-a concentration. They were compared to the spatialized reference based on the in situ chl-a collected in the reservoir of Itumbiara–Goiás in the same period. The comparison was performed through the "Spatial Language for Algebraic Geoprocessing" (LEGAL) implemented at SPRING software. Results showed a better accuracy for the procedure using the spatialization of Rrs and map algebra of them. Thus the spatializaton of proximal remote sensing measurements in order to retrieve the optically active components in water should be performed through the interpolation of the Rrs.


Introduction
Chlorophyll-a (chl-a) has been use as a parameter to calculate the trophic state index (TSIs) of an aquatic system.It also has been use as an indicator of photoautotrophic biomass which could be related to primary productivity (Huot et al., 2007).Blooms of chl-a are visible symptoms of eutrophication process in aquatic system (Mudroch, 1999).It is also related to the increased growth of algae in which cyanobacteria is a troublesome group.It is the responsible to cause severe oxygen depletion which causes fish mortalities.It is also the responsible for the death of cattle and other animals from ingestion of cyanotoxins (Melack, 2000).Another consequence of algal blooms is the high concentrations of dissolved organic carbon (DOC).Phytoplankton blooms also interferes in the decrease of light penetration through the water column (Boyer, Kelble, Ortner, & Rudnick, 2009).This process could depress seagrass growth and productivity.
Chl-a monitoring is important for water quality measures.Mainly because it may vary temporally and spatially in an aquatic system.However traditional field sampling methods used to estimate Chl-a concentration are expensive, time consuming (Duan, Ma, Xu, Y. Zhang, & B. Zhang, 2010) and difficult to make for large studies areas (Gons, 1999).These observations have highlighted the importance of the monitoring of chl-a concentration.Remote sensing has been use as a tool to estimate chl-a concentration in case 1 and case 2 waters.This water classification varies according to the optical properties of the water.In case 1 waters the optical properties are determined primarily by phytoplankton and related constituents as dissolved organic matter (CDOM) and detritus.However, in case 2 waters the optical properties can be influenced by other constituents such as mineral particles, CDOM, or microbubbles, whose concentrations do not covary with the phytoplankton concentration (Mobley, Stramski, Bisset, & Boss, 2003).

Bio-Optical Models
Bio-optical models describe the interaction between water's optically active constituents and electromagnetic radiation.The application of these algorithms to case 2 waters still a challenge due to the presence of chl-a's non-covarying constituents (S.Mishra & D. Mishra, 2012).Empiricals and semi-analytical algorithms have been developed to estimate chl-a in case 2 waters (Moses, Gitelson, Berdnikov, Saprygin, & Povazhnyi, 2012;Duan et al., 2010;S. Mishra & D. Mishra, 2012).Gilerson et al. (2010) and S. Mishra and D. Mishra (2012) developed a two band algorithm for inland waters using red-near infrared (NIR) bands.Both studies used the bands 7 and 9 from the MEdium Resolution Imaging Spectrometer (MERIS) from the European Space Agency (ESA).
These two bands were centered on the following wavelengths: 665 and 708.75 nm (ESA, 2006).The 665nm wavelength represents a wide spectral absorption peak which is generally assigned to the absorption by chl-a pigment.The reflectance peak centered at 700 nm which is maximally sensitive to the variations in chl-a concentration in water (S.Mishra & D. Mishra, 2012).The parameters of both algorithms were empirically set by comparing radiometric data and in situ measured data from a particular geographical location under a particular seasonal regime.Therefore they only represent a specific region on Earth which has a unique environmental dynamic.These only demonstrate the need for studies in tropical regions since both models were developed in aquatic systems placed in latitudes higher than 30º.
Nevertheless there is no model developed for tropical case 2 waters, which is a distinct environment with a unique primary production dynamic.In order to analyze a tropical aquatic system, we adapted two models (Gilerson et al., 2010;S. Mishra & D. Mishra, 2012) using single wavelengths from proximal remote sensing.To spatialize the results of bio-optical modeling we tested two approaches of spatial interpolation using ordinary kriging algorithms.
The purpose of our study was to evaluate the accuracy of the two-wavelengths algorithms to estimate chl-a concentration in a tropical aquatic system.It was also a goal to analyze two procedures to interpolate the bio-optical modeling results: (1) by the spatialization of the estimated reflectance for 470 and 700 nm wavelengths followed by an map algebra of them; and (2) by the spatialization of the calculated model values for the same two wavelengths.

Bio-Optical Modeling
Empirical and semi-analytical bio-optical models have been developed to estimate water constituents concentrations.Empirical models are based on calculation of statistical relation between the water constituent concentration and water leaving radiance or reflectance (Dekker, 1993).Semi-analytical models describe the reflectance from Inherent Optical Properties (IOPs) and water constituent concentration, by means of the radiative transfer equation.
To estimate the remote sensing reflectance (R rs ) from each wavelength (Gordon, Brown, & Jacobs, 1975;Gordon et al., 1988) proposed a relation where the R rs is calculated just beneath the water surface and can be expressed as: Where a(λ) and b b (λ) are the total spectral absorption and backscattering coefficients.γ is a proportionality constant which depends on the geometry of the field of light emerging from the water body.It can be estimate by 2 different ways: γ= f/Q or γ=f.Where f is the anisotropic factor of the downwelling light field (Kirk, 1994) and Q is the geometrical factor (Gons, 1999).
To determinate the f value, it was previous found that it was a function of the solar elevation angle (Kirk, 1994).
It was reasonably well expressed as a linear function of μ 0 -the mean cosine of the zenith angle-and can be estimated as: 0.975 0.629 (2) The value of μ 0 was calculated according to the sampling time, locations (latitude and longitude) and solar zenith angle (Martin & McCutcheon, 1999;Rees, 2001).For the Q value it was proposed (Gons, 1999) an empirical equation for turbid inland waters under different solar elevation angles, as expressed:

2.38
(3) A factor of 0.544 (Austin, 1980) was proposed to relate the radiance just above the surface to radiance just beneath the surface.Therefore, R rs just above the water surface can be determined as follows: , 0 0.544 (4) To calculate the a(λ) it is possible to use the following relation: (5) Where: a w (λ) is the absorption coefficient of pure water, and a ph (λ) , a tr (λ) and a CDOM (λ) are specific absorption coefficients for phytoplankton, tripton and CDOM (Colored Dissolved Organic Matter), respectively.
The chl-a bio-optical models proposed by S. Mishra and D. Mishra (2012) and Gilerson et al. (2010) use two bands from MERIS.Their models were given through a division of the reflectances.Equation 6shows the bio-optical model developed by S. Mishra and D. Mishra (2012).Equation 7shows the model developed by Gilerson et al. (2010). And:  35,75 • 19.3 .
(7) Where: C chl is the concentration of chl-a and Rλ n is the reflectance of the MERIS bands centered in the "n" wavelength.
If we apply Eauqtion 4 to estimate the reflectance to be use in both models, the fact of using a reflectance division, nulls the factors (0.544 and γ).Therefore the R rs for these two models can be estimate through the use of the coefficients of total absorption and backscattering as it is expressed in Equation 8.
Equation 8 was used to estimate the R rs in 2 wavelengths.It was possible because we used it on Equations 6 and 7 which refer to the two-bands chl-a bio-optical models.This process adapted the models which used MERIS bands to single wavelengths.
Through these procedures it is possible to calculate the R rs for the sampling points from the IOPs of the water constituents.It was also possible to calculate the bio-optical values for the same points.However there is no standard method to spatially interpolate these calculated values.This fact enhances the importance for our investigation about the spatialization of bio-optical models values for an aquatic system.

Study Site
The study was conducted in Itumbiara Reservoir (18°25' S, 49°06' W), located in Brazilian West-Central between the States of Minas Gerais and Goiás (Figure 1).The region is typically a tropical grassland savanna, and the reservoir was formed by damming of Paranaíba River resulting in the flooding of its main tributaries, Corumbá and Araguari Rivers.The geomorphology of the basin resulted in a reservoir with a dendritic pattern covering an area of approximately 814 km² and a volume of 17.03 × 109 m³ (Alcantara et al., 2010).The surface of reservoir is at 520 meters above the sea level.The length of major axis is 30 km and the maximum width is about 15 km.The depth of reservoir ranges from 0.5 m to 78 m with mean depth of 32 m.
According to Köppen (1931) Climate Classification System (CCS) the climate on the Itumbiara Reservoir region is classified as "tropical savanna" with two well defined seasons, dry (May-October) and wet (December-April).Monthly precipitation ranges from 5 mm (dry season) to 250 mm (wet season).The air temperature is high during the wet season (24 °C to 26 °C) and a little lower in the dry season (20 °C.The relative humidity has a similar pattern, but with a small shift in the minimum value towards September (47%).The wind intensity is lower during the wet season (1.6 m s -1 ) than in the dry season (3.3 m s -1 ) (Curtarelli, 2012).

Data Collection
The proximal remote sensing data and water samples for laboratory analysis were collected from 12 to 13 of May, 2009, obtained between 10:00 and 14:00 local to provide representative daily readings (Nascimento, 2010).
The samples were collected in 21 different locations in the Itumbiara's Reservoir (Figure 1).Absorptions coefficients were calculated from surface water samples were collected in 1 L Niskin bottles and immediately filtered onto 0.7 μm Whatman GF/F filters.The volume of water filtered were 250 ml.Particulate absorption coefficient, ap(λ), and absorption coefficient of detrital matter, ad(λ), were determined using standard quantitative filtration technique (QFT) as described in Mitchell, Kahru, Wieland, and Stramska (2003).A UV-2450 Shimadzu spectrophotometer with an integrating sphere was used to measure absorbance of the samples within a spectral range from 350 to 700 nm.For a CDOM (λ) the water samples analysis were filtered immediately after collection through 0.2 μm nucleopore membrane filters and read in a U-3010 Hitachi, spectrophotometer.Backscattering coefficients for 470 and 700nm were calculated using a Satlantic profiler (Satlantic, 2000) at 1m and were used to estimate R rs.
Concentration of limnological parameters-concentration of chl-a, suspended organic and inorganic matter and dissolved organic and inorganic carbon-were obtained from analytical process in the laboratory.The collected water samples from the subsurface-approximately 10 cm from the surface-were at cool temperature until the delivery at the laboratory for the analysis without duplicates.The method of chlorophyll-a analysis (Nush, 1980) consisted in filtered the collected samples in the same day of the collect and then analyzed the filters through a spectrophotometer.Total Suspended Solids (TSS) were determined based on Wetzel and Likens (1991) Water physical parameters as water temperature, pH and turbidity were also measured in situ using a YSI Inc and Horiba multi-parameter sensors.Euphotic zone was estimated multiplying the Secchi Depth by 2.7 (Cole, 1994).
The euphotic zone-or euthotic depth-is defined as the compartment where the solar radiation levels are at least 1% of the levels measured just below the free surface.It is in this depth where the primary production of lakes and reservoirs concentrates (Tundisi & Matsumura-Tundisi, 2008).

Bio-Optical Model Calculation and Validation
Using IOPs obtained from the profiler and QFT, we calculated the R rs for the two wavelengths (470 nm and 700nm) through Equation 8. Then the two-wavelengths bio-optical models values were calculated based on Equations 6 and 7 for the estimated R rs.
Although the two-band models (Eauqtions 6 and 7) used red-NIR bands, we tested wavelengths at 470 nm and 700 nm.We chose them because of the chl-a reflectance peak at 400-450 nm and 650-700 nm (Kirk, 2011).It was also considered departures from detritus form at 440 and 470nm associated with phytoplankton pigments made the absorption signals small (Ruddick, Gons, Rijkeboer, & Tilstone, 2001).At 700 nm the radiance spectrum correlated strongly with the chl-a concentration (Gitelson, 1992).It is the reflectance peak due to the minimum sum of absorption of phytoplankton, particulate and dissolved organic matter and pure water (Gitelson, 1992).
To validate the adapted algorithms we used a linear pearson correlation analysis.It was performed among the calculated two-wavelength bio-optical values and chl-a concentration.All of the values were normalized by minimum-maximum normalization.It was also calculated indicators that provide a quantitative estimate of the differences between model and the reference.Bias, Mean Square Error (MSE), Mean Absolute Error (MAE) and Root Mean Square Error (RMSE) were used evaluate the adapted models and were calculated according to Table 1.
Table 1.Error estimators and their formulas

Spatial Interpolation
The spatial interpolation of the models were conducted using a geostatistic technique called ordinary kriging.It is a generic name for a family of generalized least-squares regressions algorithms (Goovaerts, 1997).The calculation of the Kriging weights is based upon the estimation of a semivariogram model.The semivariogram was fitted using the the Geographical Information System (GIS) SPRING, version 5.2, developed at the National Institute for Space Research (INPE) (Camara, Souza, Freitas, & Garrido, 1995).It was chosen since it was possible to apply several theoretical models (spherical, exponential, Gaussian, linear and power) using the weighted least square method.
The adapted models were interpolated by two different procedures: (1) kriging technique was used to spatialized the calculated values of chl-a for each model; (2) kriging was used to spatialized the estimated R rs from 470nm and 700nm then the models were calculated by a map algebra using the spatialized R rs .
In order to compare the results, the spatial interpolated data was normalized by the minimum and maximum values in a range from 0 to 1.The bio-optical models were applied by the Spatial Language for Algebraic Geoprocessing (LEGAL) for SPRING (Camara et al., 1995).
To evaluate the precision of these two procedures to accurately retrieve the chl-a concentration, kappa statistic (Cohen, 1960) was used.The reference data was use to calculate the statistic through the SPRING GIS, and the prediction error plots were examined.
The use of Kappa statistic is a way to measure the magnitude of agreement between estimations results and the reference -the in situ data (Cohen, 1960).Then, the kappa index varies between -1 to 1, within 1 indicating perfect agreement and ≤0 indicating agreement equivalent to chance.For the calculation of kappa statistic, the estimated chlorophyll-a concentration probabilities are classified as not correlated or correlated status (≤0 or 1).It was also used the Mean Error (ME) in order to realize a spatial analysis.ME is equal to zero if the predictions were completely unbiased, i.e. centered on the measurement values.

Bio-Optical Models Validation
The results from the validation of our adapted models are expressed in Table 2. Adaptation 1 used S. Mishra and D. Mishra (2012) structure and adaptation 2 used Gilerson et al. (2010) structure.

Spatial Interpolation Results
The first method of spatial interpolation (calculated bio-optical models) to the estimated values of chl-a by the adapted bio-optical models and the reference are showed in Figure 2. Figure 2a shows the reference values spatialized by ordinary kriging to all the reservoir of Itumbiara.Figures 2b and 2c show the estimated concentration of chl-a from the adaptation of two-wavelengths bio-optical models presented.These results represented the first procedure using the calculated models values to estimate the chl-a.
Figure 3 shows the results from the second method of spatial interpolation (spatialization of R rs and then map algebra of them).The reference (Figure 3a) and estimations of chl-a concentration by two-wavelenghts models (Figures 3b and 3c) were interpolated for the spatial distribution of the reservoir.
Both figures (Figure 2 and Figure 3) illustrate the distribution of chl-a in the Itumbiara Reservoir.It revealed 2 regions which have different dynamics.The first region (1) indentified is located near the entrance of Corumbá River.It was a region with low concentration of chl-a, probably due to its clean water and low input of nutrients.
The second region (2) had a high concentration of chl-a.It was located near the entrance of Paranaíba River-in the border of the reservoir.The high concentration probably happened due to the fact that the river follows a meandric shape.This shape was the responsible for the accumulation of the organic matter in decomposition which were transported into the meander.Kappa statistical results suggest that the four proposed spatial interpolations were not correlated to the reference.
Furthermore they also suggest that the spatialization of R rs and map algebra method for the estimation got the best kappa statistic when compared to the calculated model values method.These results revealed that the adaptation 2 was more appropriated for a tropical and deep reservoir.Probably due to the structure of the algorithm which is more accurate than the first adaptaion.This fact can be explained by the algorithms themselves.The second algorithm (Gilerson et al., 2010) was calibrated with data from previous works and from the literature.In the other hand, the first algorithm (S.Mishra & D. Mishra, 2012) was a normalized difference index proposal, so it just depends on the reflectance of the wavelength.
We also tried to understand the reason for this poor statistical performance which is probably related to the reflectance estimation.It could be affected by innumerous factors like the instability on the aquatic system surface.It is also possible that Gordon's algorithm (1988) is not useful for case 2 waters since it was developed for case 1 waters.
A previous study (Nixdorf, 1994) reported that the reflectance on the surface of a water body depends on physical agents.These agents can be related to the water column depth and the bottom type which are the responsible for the magnitude and spectral composition of the backscattering flux.
To analyze the spatial changes among the estimations and the reference we used a mean error (ME) approach.It was calculated through a simple spatial difference between the estimations and the reference using LEGAL.
Figure 4 shows the mean error for the two adaptations models following the two spatial interpolations methods.
The results from the first method (calculated models values) are represented in Figure 4a and 4b.The ME results for the second method (spatialization of R rs and map algebra) are represented in Figures 4c and 4d.
As well as the kappa analysis, adaptation 2 got the best results for both methods.However the spatialization of R rs and map algebra method had the best results for the braces of the reservoir.The calculated model values method had the best results to the main body of the reservoir.Furthermore both two adaptation and methods had registered a good accuracy in the region with a high concentration of chl-a (region 2 in Figures 2, 3, and 4).But they had a low accuracy in the middle of the reservoir (region 3 in Figure 4) which is the region that had the lowest levels of chl-a concentration.
The differences in performance at each region can be related to the instability of the aquatic system' surface.The instable areas can be correlated to reservoirs' bathymetry (Figure 5) since it influences the water column stability.Nixdorf (1994) established that stability of a water column is determined by its depth and the gradients of temperature, i.e. the differences in density.Stability is one of the various factors contributing to the stratification and mixing of an aquatic system.It is also a measure for physical forcing on phytoplankton communities.Changes in the aquatic system surface can spoil the estimation of remote sensing reflectance.Besides the deepness of the reservoir, another important factor to the stability of the water column and to R rs estimation is the wind intensity.For a water surface illuminated by irradiance propagating in a direction, it is necessary to compute the irradiance reflectance as a function of wind speed (Cox & Munk, 1954).A Monte Carlo simulation was developed to estimate the upwelling radiance varying the angles and wind speeds (Morel & Gentili, 1996).Through it the authors circumvented their results by using the results for a wind speed of 0. It was noticed that it could lead to significant error for large viewing angles and wind speeds.
Figure 6 shows the spatialization of the euphotic zone-or euphotic depth-for the reservoir.We noticed a dependence of chl-a to the highest values of the euphotic zone-regions 2 and 4 (in Figure 6).In region 1, the euphotic zone was shallow mainly due to the entrance of sediments from the Corumbá River.Euphotic zone can also be related to low concentration of chl-a in aquatic systems, since it also have a daily variation due to meteorological events as cold fronts, circulation periods and water column stability (Tundisi & Matsumura-Tundisi, 2008).These results demonstrate the importance for a parameterization of a bio-optical model for tropical complex aquatic systems.The lack of studies for tropical inland waters is important because of the complexity of the environment.They have different atmospheric dynamics, dissolved matter and depth when compared to conventional studies in shallow, small and controlled reservoirs or lakes in latitudes higher than 30º.
The tropical reservoir of Itumbiara is not optically complex as others Brazilian reservoirs.However, due to the instability on the aquatic system surface-due to wind intensity and the reservoir's morphology-we could not estimate the chl-a concentration in some deep regions of the reservoir.The opposite occurred to shallow regions.This implies that the highest errors of the adaptations and even the original models we have tested had a correlation to the most instable regions of the reservoir.
The spatialization of the reference values of chl-a concentration showed a dependence to the instability in the water column.It is related to the mixing process in aquatic systems which generally depends on thermal stratification and wind-induced circulation (Alcântara et al., 2011;Tundisi et al., 2004Tundisi et al., , 2010)).This process can be accepted as a hydrological option for controlling phytoplankton blooms since it physically breaks the blooms.
In the other hand, to the radiometric data, the instability in the water column is totally miscorrelated to the accuracy of the retrieval of R rs .The equations proposed to estimate the R rs from aquatic systems were adequate to estimate the shallow-water reflectance and not deeper waters (Marioneta, Morel, & Gentili, 1994).The spectral shape of the surface reflectance was a function of the spatial variability of bottom albedo, the turbidity of the water and to depth (Burns, Taylor, & Sidhu, 2010).

Conclusion
A proximal remote sensing was used to calculate the bio-optical and R rs for the water surface of Itumbiara's reservoir.Two approaches to interpolate the proximal remote sensing data were proposed and analyzed.The results showed that the spatialization of the reflectance had better results when compared to the spatialization of bio-optical model's values.The difference noted between the both approaches is the range used in the spatialization.Since the range of the R rs values were smaller than the range of the bio-optical models values.
More extensive work needs to be done on aquatic systems from various geographic locations, depths and dissolved matters.However the use of proximal remote sensing will allow us to fully understand and quantify the limits of the potential universal applicability of the Remote Sensing products to tropical optically complex aquatic systems.

Figure 1 .
Figure 1.Localization of Itumbiara Reservoir and the distribution of the chl-a samples points

Figure 2 .
Figure 2. (a) Chl-a reference data.(b) Spatial interpolated chl-a concentration based on the spatialization of calculated adaptation 1 model values.(c) Spatial interpolatedf chl-a concentration based on the spatialization of calculated adaptation 2 model values

Figure 4 .
Figure 4. (a) Error of chl-a estimation from adaptation 1 by calculated model values method.(b) Error of chl-a estimation from adaptation 2 by calculated model values method.(c) Error of chl-a estimation from adaptation 1 by spatialization of R rs and map algebra method.(d) Error of chl-a estimation from adaptation 2 by spatialization of R rs and map algebra method

Table 2 .
Adapted models validation