A Multivariate Analysis Framework to Detect Key Environmental Factors Affecting Spatiotemporal Variability of Chlorophyll-a in a Tropical Productive Estuarine-Lagoon System

Here, we demonstrate how a combination of three multivariate statistic techniques can identify key environmental factors affecting the seasonal and spatial variability of chlorophyll-a (Chl-a) in a productive tropical estuarine-lagoon system. Remote estimation of Chl-a was carried out using a NIR-Red model based on MODIS bands, which is highly consistent with the in situ measurement of Chl-a with root mean square error (RMSE) of 15.24 mg m−3 and 13.43 mg m−3 for two independent datasets used for the model’s calibration and validation, respectively. Our findings suggest that the river discharges and hydraulic residence time of the lagoons promote a stronger effect on the spatial variability of Chl-a in the coastal lagoons, while wind, solar radiation and temperature have a secondary importance. The results also indicate a slight seasonal variability of Chl-a in Mundaú lagoon, which are different the from Manguaba lagoon. The multivariate approach was able to fully understand the relative importance of key environmental factors on the spatiotemporal variability of Chl-a of the aquatic ecosystem, providing a powerful tool for reducing dimensionality and analyzing large amounts of satellite-derived Chl-a data.


Introduction
Coastal lagoons are transitional ecosystems that exhibit large spatial variability and seasonality controlled by both continental and oceanic drives, such as river discharge, salinity, water temperature, and hydrodynamics [1,2].The limited water exchange promotes effective use and recycling of the nutrient inputs, contributing to high levels of biological productivity (e.g., phytoplankton, zooplankton, mussels and fishes) in coastal lagoons [3].On the other hand, the low water exchange rates also makes these water bodies vulnerable to eutrophication and contamination.This degradation process usually results in an increase of water turbidity due to blooms of cyanobacteria or green algae, affecting the entire trophic structure [4].Thus, the advancement of methodological frameworks to quantify eutrophication in different scales in time and space are important to improve the understanding of the ecosystem dynamics of coastal lagoons and to develop tools for accurate decision-making [5][6][7].16 36 * tidal prism is the volume of water in an estuary or inlet between mean high tide and mean low tide, or the volume of water leaving an estuary at ebb tide.

Mundaú lagoon. Additional physical characteristics of the Mundaú and Manguaba lagoons can be observed in
Figure 1.Mundaú-Manguaba Estuarine-Lagoon System (MMELS) study site, land-use cover in surrounding adjacent basins in MMELS (adapted from Guimarães Júnior, et al. [29]), as well as spatial distribution of sampling stations, which were used to collect water samples.
The eastern margin of Mundaú lagoon is predominantly occupied by the urban area of Maceió city (ca.900,000 inhabitants), the western margin is dominantly occupied mostly by coconut crops and two small urban areas (Santa Luzia do Norte and Coqueiro Seco with a population of approximately 7000 and 5000 inhabitants, respectively).Moreover, sugarcane crops cover the coastal tablelands surrounding the Mundaú lagoon.In Manguaba lagoon, both eastern and western land margins are intensively used for agriculture (i.e., sugarcane and coconut crops), grassland, and vegetation fragments (forest and mangrove areas).The urban areas of Pilar (ca.30,000 inhabitants in the northern margin) and Marechal Deodoro (ca.46,000 inhabitants in the southwestern margin) contribute with a significant nutrient loading and organic matter deposit into the Manguaba lagoon.Besides the surrounding cities, high organic and nutrient sources in both lagoons are the result of the Figure 1.Mundaú-Manguaba Estuarine-Lagoon System (MMELS) study site, land-use cover in surrounding adjacent basins in MMELS (adapted from Guimarães Júnior, et al. [29]), as well as spatial distribution of sampling stations, which were used to collect water samples.[28].

Features Mundaú Manguaba
Area (km The eastern margin of Mundaú lagoon is predominantly occupied by the urban area of Maceió city (ca.900,000 inhabitants), the western margin is dominantly occupied mostly by coconut crops and two small urban areas (Santa Luzia do Norte and Coqueiro Seco with a population of approximately 7000 and 5000 inhabitants, respectively).Moreover, sugarcane crops cover the coastal tablelands surrounding the Mundaú lagoon.In Manguaba lagoon, both eastern and western land margins are intensively used for agriculture (i.e., sugarcane and coconut crops), grassland, and vegetation fragments (forest and mangrove areas).The urban areas of Pilar (ca.30,000 inhabitants in the northern margin) and Marechal Deodoro (ca.46,000 inhabitants in the southwestern margin) contribute with a significant nutrient loading and organic matter deposit into the Manguaba lagoon.Besides the surrounding cities, high organic and nutrient sources in both lagoons are the result of the release of untreated effluents from the Mundaú and Paraíba do Meio river basins.As result, both lagoons are characterized by eutrophic conditions and dominated by phytoplankton [30].

Water Sampling and Analysis
Shipboard data were collected during eleven field campaigns (Table 2).In each field campaign, a set of sampling stations spatially well distributed was defined in order to collect water samples for laboratory analysis (see Figure 1), representing a total of 171 surface water samples.Each 2-liter bottle was collected 0.2 m bellow the water surface immediately before satellite overpass.The samples were stored in a cooler with ice in the dark, and taken back to the laboratory for analyzing the concentration of Chl-a and total suspended solid (TSS).For each water sample, Chl-a was extracted from Whatman GF/F filters into 90% ethanol for 18 h in an amber flask and measured by the spectrophotometric trichromatic method [31].For TSS determination, Whatman GF/F filters were dried, before and after water filtering, to a constant weight at 103 to 105

Satellite Data and Image Processing
In this study, we used MODIS-Terra (MOD09A1) and MODIS-Aqua (MYD09A1) 500-m spatial resolution Level-3 Surface Reflectance 8-day composites, which are available at http://LPDAAC.usgs.gov.These MODIS products were selected in this study, since they are the best daily observation within 8-days interval for each pixel according some image quality criteria, such as viewing angle, presence of clouds, and aerosols quantity.In addition, these MODIS products use a global atmospheric correction algorithm, which corrects for the effects of gaseous and aerosol scattering and absorption as well as adjacency effects caused by variation of land cover, Bidirectional Reflectance Distribution Function (BRDF) and atmosphere coupling effects, and contamination by thin cirrus [32].Image download and processing were made using the GetMODIS and MODIS River Reflectance Retrieval (MOD3R), software developed by the Institut de Recherche pour le Développement (IRD).GetMODIS was used to download the MODIS products between 2000 and 2016 for all tiles that overlap MMELS domain, mosaic in case of multiple tiles for a scene, clip to the domain extent, and resample to the domain grid cells using a mask polygon.Pixels with horizontal distance between their centers and margin below to 500 m (image resolution) were not considered in this analysis in order to avoid or minimize the effect of margins and bottom sediment on spectral reflectance.In addition, MOD3R was used to obtain the surface reflectance from the first five spectral bands of the MOD09A1 and MYD09A1 products in each pixel of the MMELS domain, as well as to assess image quality.

NIR-Red Model Based on MODIS Bands to Retrieve Chl-a
An empirical algorithm based on MODIS spectral bands was established to retrieve Chl-a concentrations in MMELS.In order to evaluate the accuracy and stability of the best models based on MODIS remote sensing reflectance, the Chl-a dataset was divided into independent calibration and validation subsets according to Table 2. Despite the high number of collected water samples, only 28% of the sampling stations were used to calibrate and validate the NIR-Red band ratio algorithm to retrieve Chl-a in MMELS either due to the presence of clouds, shadows, the proximity of margins on the measurement sites, or due to the bad quality of the images from MODIS quality product control.This high number of discarded sampling stations reveals how challenging it is to develop bio-optical algorithms based on satellite data for coastal tropical inlands, which exhibit a high cloud cover index throughout the year.
An algorithm was developed in order to identify the optimal relationship between band ratio (λ 2 /λ 1 ) of MODIS remote sensing reflectance and Chl-a concentration.This algorithm consists of testing all possible sets of band ratio best-fit functions (e.g., linear, potential, exponential) between observed and estimated Chl-a within the five spectral bands of the MOD09A1 and MYD09A1.Finally, a NIR-Red band ratio model (Rrs (859)/Rrs (645)) for retrieving Chl-a in MMELS was confirmed (Figure 2a), showing a higher correlation between measured and estimated Chl-a concentration for the MMELS calibration dataset (r 2 = 0.83 and RMSE = 15.24mg m −3 ) compared to other tested models, such as the standard Blue-Green band ratio model (r 2 = 0.11 and RMSE = 40.25 mg m −3 ), which is usually applied to Case 1 waters.All calibration data in the Chl-a retrieval algorithm for MMELS were in the 95% confidence line of the linear regression.Due to the particular optical characteristics of waters in each lagoon, the calibration dataset allowed us to propose a unique model for the entire system with a wider range of Chl-a concentrations.The validation result indicates that the NIR-Red band ratio model may retrieve Chl-a successfully with relatively high accuracy (r 2 = 0.75 and RMSE = 13.43 mg m −3 , Figure 2b).However, considering the validation dataset, this model seems to underestimate Chl-a values in the low range, which was not observed in the calibration dataset.Thus, the usage of this model to estimate low values of Chl-a (<25 µg/L) may be associated to higher uncertainties.

NIR-Red Model Based on MODIS Bands to Retrieve Chl-a
An empirical algorithm based on MODIS spectral bands was established to retrieve Chl-a concentrations in MMELS.In order to evaluate the accuracy and stability of the best models based on MODIS remote sensing reflectance, the Chl-a dataset was divided into independent calibration and validation subsets according to Table 2. Despite the high number of collected water samples, only 28% of the sampling stations were used to calibrate and validate the NIR-Red band ratio algorithm to retrieve Chl-a in MMELS either due to the presence of clouds, shadows, the proximity of margins on the measurement sites, or due to the bad quality of the images from MODIS quality product control.This high number of discarded sampling stations reveals how challenging it is to develop bio-optical algorithms based on satellite data for coastal tropical inlands, which exhibit a high cloud cover index throughout the year.
An algorithm was developed in order to identify the optimal relationship between band ratio (λ2/λ1) of MODIS remote sensing reflectance and Chl-a concentration.This algorithm consists of testing all possible sets of band ratio best-fit functions (e.g., linear, potential, exponential) between observed and estimated Chl-a within the five spectral bands of the MOD09A1 and MYD09A1.Finally, a NIR-Red band ratio model (Rrs (859)/Rrs (645)) for retrieving Chl-a in MMELS was confirmed (Figure 2a), showing a higher correlation between measured and estimated Chl-a concentration for the MMELS calibration dataset (r 2 = 0.83 and RMSE = 15.24mg m −3 ) compared to other tested models, such as the standard Blue-Green band ratio model (r 2 = 0.11 and RMSE = 40.25 mg m −3 ), which is usually applied to Case 1 waters.All calibration data in the Chl-a retrieval algorithm for MMELS were in the 95% confidence line of the linear regression.Due to the particular optical characteristics of waters in each lagoon, the calibration dataset allowed us to propose a unique model for the entire system with a wider range of Chl-a concentrations.The validation result indicates that the NIR-Red band ratio model may retrieve Chl-a successfully with relatively high accuracy (r 2 = 0.75 and RMSE = 13.43 mg m −3 , Figure 2b).However, considering the validation dataset, this model seems to underestimate Chl-a values in the low range, which was not observed in the calibration dataset.Thus, the usage of this model to estimate low values of Chl-a (<25 µ g/L) may be associated to higher uncertainties.

Time-Series of Chl-a and Environmental Variables
The monthly mean Chl-a value for seventeen years (2000-2016) was calculated in each pixel of the MMELS domain from MODIS-derived 8-day Chl-a data using the NIR-Red band ratio model, which was presented in the previous section.
In order to assess cause-and-effect relationships between environmental variables and spatiotemporal variability of Chl-a in MMELS, we obtained the historical time-series of air temperature, wind intensity, solar radiation, precipitation, tidal levels, and main flow tributaries.Meteorological data were available from a meteorological station located 7.5 km northeast of the Mundaú lagoon center.This station is maintained by the Instituto Nacional de Meteorologia (INMET), providing hourly meteorological data.Information on the hourly tidal levels was provided by the Navy Oceanographic Center (CHN) and was based on a harmonic analysis of water level data from Maceió Harbor (9 • 40.968 S, 35 • 43.424 W), located 9 km from the MMELS inlet.Two river gauge stations located near the Mundaú and Paraíba do Meio river mouths were selected, providing daily flow data.For each environmental variable, the monthly mean values were calculated based on their original time-series.Finally, an eigenvector-based filtering technique [33] was used to fill missing data in these time-series.

Spatiotemporal Patterns of Chl-a and Their Key Environmental Factors
A combination of three multivariate statistical techniques was used in order to identify how key environmental factors affect seasonal and spatial variability of Chl-a in MMELS.These techniques generally require a complete time series of input maps without data voids.Therefore, a method to eliminate missing data based on the eigenvector-based filtering technique decomposition was applied to a large dataset to obtain complete a time series of Chl-a maps.However, in order to minimize the effect of the replaced values on the spatiotemporal results of Chl-a, initially we applied a quality control to discard satellite-derived monthly Chl-a maps exhibiting more than 25% missing pixels in the MODIS images over the MMELS.Moreover, for each pixel, time series exhibiting more than 5% of missing data were discarded, resulting in 53 and 79 valid pixels (i.e., pixels with a complete time series of Chl-a) for Mundaú and Manguaba lagoons, respectively.The remaining missing data (less than 1% of total data) were replaced using an eigenvector-based filtering technique, which is a statistical method based on principal component analysis according to Ibanez and Conversi [33].
Principal Component Analysis (PCA) was used to identify the main Chl-a patterns of spatial variability.This analysis is useful for reducing the information of a large number of variables into a smaller set of orthogonal variables (modes or components) that represent the major part of the overall data variance.The complete theory of PCA may be found in Jackson [34].Each PCA mode is composed by an eigenvector and an eigenvalue.The eigenvectors represent the co-variability of each variable, whereas the eigenvalues express how much of the original variance is explained by each eigenvector.PCA was applied to the satellite-derived Chl-a data, containing the estimated time-series of Chl-a in each valid pixel for both lagoons (i.e., 53 pixels × 139 months and 79 pixels × 83 months for Mundaú and Manguaba lagoons, respectively).Arrays of eigenvectors of the first four components, which explained more than 70% of the total variance, were evaluated to identify areas with similar patterns of Chl-a variability.In addition, the principal components were rotated using the varimax procedure in order to relax the orthogonality constraints of the PCA and to improve the understanding of the spatial patterns [35].Environmental factors that may explain the observed patterns in these first components are discussed.
A cluster analysis was performed to identify Homogeneous Spatial Groups (HSG) of Chl-a with similar variability, complementing PCA analysis.Before this analysis, the data was filtered using the rotated principal components so that at least 70% of the original data variability was maintained.This filtering was made in order to minimize the interference of noise on the cluster analysis [36].The method of Ward (i.e., a criterion applied in hierarchical cluster analysis) was applied to establish the numbers of HSG [37], and the clustering validation was based on the physical understanding of the geographical expression of the cluster in the MMELS.
Canonical Correspondence Analysis (CCA) was used to assess the relationship between environmental variables and satellite-derived Chl-a data (i.e., time-series with monthly average values in each valid pixel).A complete description of CCA can be seen in Braak [38].Specifically, CCA was used to identify which environmental variables are important in the spatiotemporal distribution of Chl-a.The derived CCA scores (eigenvectors) were classified considering temporal seasonality (rainy and dry seasons) and HSG obtained from cluster analysis.
All multivariate analyses were performed using the software Palaeontological Statistics (PAST), which is freely available to download at https://folk.uio.no/ohammer/past/, and output maps were produced in ArcGIS ® 10.2.

Constituent Concentrations
Notably, the Chl-a and TSS concentrations in the water samples exhibited high variability (Table 3), which is a particular characteristic of this system (i.e., high spatial and temporal gradients of Chl-a).The levels of Chl-a differed between the lagoons, ranging between 0.97 and 48.9 mg m −3 (mean and SD values of 12.86 and 9.72 mg m −3 , respectively) in Mundaú lagoon, and between 5.99 and 117.54 mg m −3 in Manguaba lagoon (mean and SD values of 42.77 and 24.22 mg m −3 , respectively).In contrast to the Chl-a concentrations, the TSS concentrations were higher in Mundaú lagoon (range of 15.2-61.0mg L −1 and mean value of 32.8 mg L −1 ) in comparison with Manguaba lagoon (range of 9.0-44.0mg L −1 and mean value of 22.7 mg L −1 ).In addition, Chl-a and TSS concentrations showed low correlation in the MMELS (r 2 = 0.08; data not shown), suggesting MMELS can be characterized as an optically complex environment (Case 2 water).

Spatiotemporal Variability of Chl-a and Driven Factors
The first four modes of the rotated PCA explained together at least 70% of the total variance of the satellite-derived Chl-a for both lagoons (Table 4).The first three principal components (i.e., PC1, PC2 and PC3) were found to be the most significant ones, both statistically and physically.The PCA identified individual spatial patterns of the eigenvectors for each principal component in both lagoons (Figure 3).For Mundaú lagoon, the first rotated component of the PCA (Figure 3a) reveals a transverse gradient (east-west direction), consisting of negative eigenvector values throughout the lagoon.The second rotated component shows a longitudinal gradient (north-south direction), with positive values in the southern and central part of the lagoon and negative values concentrated in the northern lagoon (Figure 3b).The third rotated component is mostly characterized by negative values of the eigenvectors, which are concentrated in the southern and northern regions of the lagoon, close to the main water inflows, while the positive values mainly highlight the central part of the Mundaú lagoon (Figure 3c).
For Manguaba lagoon, the first rotated component of the PCA is constituted by negative values of eigenvectors crossing a well-defined transversal gradient (east-west direction), which covers the center and south part of the lagoon (Figure 3d).In addition, northern lagoon shows a slightly gradient in comparison to other regions.The second rotated component shows a longitudinal gradient with positive values in the center and south part of the lagoon, as well as negative values mainly highlighted in the northern Manguaba lagoon (Figure 3e).The third rotated component is majority characterized by positive values of eigenvectors, which are concentrated in the central and southern regions of the Manguaba lagoon (Figure 3f).In addition, there is no significant variability of eigenvector values in the north center part of Manguaba lagoon (Figure 3f).
The cluster analysis using the first four components of PCA permitted identifying four well-defined HSG of Chl-a variability in each lagoon (Figure 4).All HSG identified by the cluster analysis are coherent with the spatial patterns presented by the rotated PCA.In general, the groups were spatially disposed in the longitudinal direction in both lagoons (Figure 4), following a similar characteristic pointed out by the second mode of the rotated PCA.For Mundaú lagoon, the stations that compose HSG1 are located in the northern region and two large groups (HSG2 and HSG3) comprise stations located on the central region of the lagoon.In addition, a fourth group was identified in the southern lagoon (HSG4), with a few cells distributed over the west margin.There was no significant gradient of average values of Chl-a among the groups in Mundaú lagoon, only a slight increase of average values of Chl-a in HG1 and HG4 in comparison to central HSG (Table 5).In Manguaba lagoon (Figure 4b), two large groups (HSG2 and HSG1) have been found in the central and northern part of the lagoon.The cluster analysis also properly identified two small groups in the southern lagoon (HSG4) and in the west margin (HSG3).A longitudinal Chl-a gradient was observed from south (HSG4 with higher values) to north (HSG1 with lower values) part of the lagoon (Table 5).
The monthly variability of the average values of estimated Chl-a in each HSG may be observed in Figure 5.In general, Manguaba lagoon has a more definite seasonality in comparison to Mundaú lagoon, with higher Chl-a values being observed in the wet season (April to June).In Mundaú lagoon, higher values of Chl-a were observed in the dry season (October to March).The northern Manguaba lagoon (HSG1) has a similar temporal variability of Chl-a to Mundaú lagoon.
The derived CCA scores indicated a slight seasonal variability of Chl-a in Mundaú lagoon (Figure 6a), making it difficult to identify the relative importance of environmental factors on the temporal variability of Chl-a.Environmental factors such as precipitation, river discharge, and tide variability are mostly associated with spatial groups located at northern portions of the Mundaú lagoon (groupings 1 and 2), whereas wind and solar radiation affect a large portion of points in HSG3.CCA analysis did not reveal a notable relationship between environmental factors and variability of Chl-a in the HSG4, being air temperature seems to have a larger relative importance in this HSG.
Differently from Mundaú lagoon, CCA indicated more defined seasonal variability of Chl-a in Manguaba lagoon (Figure 6b).In the dry season, environmental factors such as wind, solar radiation, and air temperature have a major influence on Chl-a while precipitation, river discharge, and tide variability play an important role on Chl-a variability in the rainy season.Chl-a variability in the northern portion of the Manguaba lagoon is mainly affected by precipitation and river discharge.CCA analysis did not reveal a clear relationship between environmental factors and variability of Chl-a in the central portion of the Manguaba lagoon (groupings 2 and 3).CCA scores also indicated that air temperature, wind, and solar radiation influence the southern portion of Manguaba lagoon.Temporal scores were split into dry (o) and rainy (+) seasons, respectively.Spatial scores (•) were split by the HSG obtained from cluster analysis (see Figure 4).

Discussion
The first component of the rotated PCA indicated a transverse gradient of Chl-a variability between the west and east margins for both lagoons.This main pattern of variability may be associated to the Hydraulic Residence Time (HRT), since it is a determining factor in spatial distribution processes and nutrient availability, affecting the phytoplankton development, and, therefore, it may explain the spatiotemporal changes in Chl-a concentration on the water surface [39].
In general, the HRT in Mundaú lagoon is lower than in Manguaba lagoon [28,40], since Manguaba lagoon has a lower influence of tidal forces due to the narrower channels compared to channel systems that connect Mundaú lagoon to the inlet.Thus, the flow rate of seawater to Manguaba lagoon is reduced.In Mundaú lagoon, the lower HRT values (i.e., ranging from 6 h in the south to 18 days in the north) are found near the west margin, which indicates a larger circulation in this region from the mouth of the Mundaú River to the entrance of the channel that links the lagoon to the ocean [41].On the other hand, higher HRT values (between 23 and 27 days) were found in the most central region and close to east margin, decreasing water circulation and promoting less variability in Chl-a concentrations (see standard deviation in Table 5).The spatial variability of the Temporal scores were split into dry (o) and rainy (+) seasons, respectively.Spatial scores (•) were split by the HSG obtained from cluster analysis (see Figure 4).

Discussion
The first component of the rotated PCA indicated a transverse gradient of Chl-a variability between the west and east margins for both lagoons.This main pattern of variability may be associated to the Hydraulic Residence Time (HRT), since it is a determining factor in spatial distribution processes and nutrient availability, affecting the phytoplankton development, and, therefore, it may explain the spatiotemporal changes in Chl-a concentration on the water surface [39].
In general, the HRT in Mundaú lagoon is lower than in Manguaba lagoon [28,40], since Manguaba lagoon has a lower influence of tidal forces due to the narrower channels compared to channel systems that connect Mundaú lagoon to the inlet.Thus, the flow rate of seawater to Manguaba lagoon is reduced.In Mundaú lagoon, the lower HRT values (i.e., ranging from 6 h in the south to 18 days in the north) are found near the west margin, which indicates a larger circulation in this region from the mouth of the Mundaú River to the entrance of the channel that links the lagoon to the ocean [41].On the other hand, higher HRT values (between 23 and 27 days) were found in the most central region and close to east margin, decreasing water circulation and promoting less variability in Chl-a concentrations (see standard deviation in Table 5).The spatial variability of the ) were split by the HSG obtained from cluster analysis (see Figure 4).

Discussion
The first component of the rotated PCA indicated a transverse gradient of Chl-a variability between the west and east margins for both lagoons.This main pattern of variability may be associated to the Hydraulic Residence Time (HRT), since it is a determining factor in spatial distribution processes and nutrient availability, affecting the phytoplankton development, and, therefore, it may explain the spatiotemporal changes in Chl-a concentration on the water surface [39].
In general, the HRT in Mundaú lagoon is lower than in Manguaba lagoon [28,40], since Manguaba lagoon has a lower influence of tidal forces due to the narrower channels compared to channel systems that connect Mundaú lagoon to the inlet.Thus, the flow rate of seawater to Manguaba lagoon is reduced.In Mundaú lagoon, the lower HRT values (i.e., ranging from 6 h in the south to 18 days in the north) are found near the west margin, which indicates a larger circulation in this region from the mouth of the Mundaú River to the entrance of the channel that links the lagoon to the ocean [41].On the other hand, higher HRT values (between 23 and 27 days) were found in the most central region and close to east margin, decreasing water circulation and promoting less variability in Chl-a concentrations (see standard deviation in Table 5).The spatial variability of the HRT in Manguaba Lagoon was recently estimated by Souza [40].However, this study neglected the effect of the wind on water circulation.Previous studies reported that wind might play an important role on the mixing and spatial patterns of hydrodynamics in Mundaú and Manguaba lagoons [28,42].Therefore, further efforts are needed as an attempt to associate the spatial distribution of Chl-a and HRT in Manguaba Lagoon.
Besides HRT, another factor that may influence the variability pattern of Chl-a is the lateral contribution of organic load due to point and non-point sources from the land use and urban occupation surrounding the lagoons.Costa, et al. [43] reported that the main external sources of particulate organic matter to the lagoons (allochthonous source) during the rainy season are untreated sewage discharged from Santa Luzia, Coqueiro Seco, and Maceió cities for Mundaú lagoon, and fertilizer and nutrient loading, probably coming from sugarcane crops, for Manguaba lagoon.That work also identified that sewage discharge was practically nonexistent into Manguaba lagoon, although minor wastewater discharges have been observed near the mouths of the Sumaúma and Paraíba do Meio Rivers.Despite the external sources of particulate organic material in suspension, the major source is generated within the lagoons (autochthonous source), resulting from the phytoplankton mortality [43].This suggests that the conditions of water quality are being affected more strongly by the HRT than by the nutrient loads, as identified by Wainger, et al. [44] in small sub-estuaries of Chesapeake Bay.
The second principal component indicated a longitudinal variability of Chl-a concentrations in both lagoons, since the values of eigenvectors pointing in opposite directions from the northern (close to the main tributaries) to the southern portion (close to ocean inlets) of the lagoons.In Mundaú lagoon, this component may be capturing the effects of fluvial discharge and tides at the north and south extremities, respectively, where highest Chl-a concentrations have been found.Indeed, CCA has also shown that Chl-a variability in the northern region of Mundaú Lagoon is strongly influenced by the Mundaú River discharge and precipitation.Thus, land-use and land-cover of the Mundaú river basin play an important role in the input of freshwater, organic matter, and nutrients to the lagoon, favoring the development of phytoplankton species more adapted to less saline environments.As this region does not show an excessive algal growth, it is suggested that the main source of total suspended solids (TSS) (i.e., solids discharged by Mundaú River, see Oliveira and Kjerfve [28]) may be limiting the algal growth in this region.In the central eastern region of Mundaú Lagoon (groups HSG2 and HSG3), where HRT are higher, Chl-a concentrations showed a slight variation in comparison with regions occupied by groups HSG1 and HSG4.In this region, other environmental factors such as wind, temperature, and solar radiation seem to control this variability (see Figure 6a).Despite the high HRT and constant nutrient supply (from domestic sewage and/or from phytoplankton decomposition), high concentrations of Chl-a (e.g., above 60 µg/L) were not observed in this central region.This limitation may be attributed to higher mortality and predation rates by bivalve mussels (particularly by the charru mussel, Mytella charruana, Maioli, et al. [45]) and also to a larger growth limitation due to light inhibition in the water, restricting the development of phytoplankton in this region.
Melo-Magalhães et al. [29] found that most of the phytoplankton species identified in the Mundaú and Manguaba lagoons are freshwater (51%) species, followed by species of the neritic or neritic/marine (47%), and only 2% of phytoplankton species have an estuarine origin.In Mundaú Lagoon, the authors also identified a domain of diatom species, mainly Skeletonema cf.costatum, in the regions with higher concentrations of salinity.It is well-known that the growth of diatom species is more difficult in environments with high levels of water turbidity [46] than other phytoplankton species such as cyanobacteria.Moreover, diatom species are the main source of food for the charru mussel [47], which is found only in Mundaú Lagoon [28].
In the southern area of Mundaú lagoon, which receives a larger contribution of sewage from the Maceió city and has a lower HRT, the causal relationship for Chl-a variability was not evidenced by CCA.In this region, a tide prism probably controls the patterns of variability of the Chl-a concentration.However, the time scale adopted in the study (monthly average) may not have been suitable to detect faster responses, such as diurnal and semi-diurnal variations in the tidal range (6 h between the highest and lowest levels), impairing the cause and effect analysis.Moreover, the CCA did not indicate a well-defined temporal variability of Chl-a between dry and rainy seasons over the analyzed period in Mundaú lagoon.This fact corroborates with the results obtained during the field campaigns and past studies carried out by other authors during the period considered in this study [48].
In Manguaba lagoon, where the tidal prism have a less marked effect than in Mundaú Lagoon [28], the longitudinal gradient of Chl-a seems to be associated with the hydrological seasonality (i.e., dry and rainy seasons).In the northern end of Manguaba Lagoon (HSG1), environmental factors such as wind, solar radiation, and temperature play an important role in Chl-a variability during the dry season (Figure 6b).During the dry season, the development of phytoplankton species, which are better adapted to freshwater, seems to be more closely associated to the lower input of inorganic particulate material from the Paraiba do Meio River and to the high autochthonous production of nutrients in comparison with Mundaú lagoon.In the rainy season, the enhancement of fluvial and solids discharge rates, possibly, limits phytoplankton growth in the northern region of Manguaba Lagoon when compared to the other regions.
The central region of Manguaba lagoon (HSG2 and HSG3) corresponds to a transition area that is influenced mainly by environmental factors closely associated with the dry season (i.e., temperature, solar radiation, and wind) and, to a lesser effect, by the tide.The central region shows much higher concentrations of Chl-a than in Mundaú Lagoon, which is probably due not only to the lower concentrations of suspended inorganic particulate matter here than in the northern region but also, and mainly, to the type of phytoplankton species dominating in this region.Melo-Magalhães et al. [29] reported that the dominant species in Manguaba lagoon are freshwater cyanophytes, mainly Anabaena spiroides in the dry season, Microcystis aeruginosa in the rainy season, and diatom species, predominantly Cyclotella meneghiniana.
In the southern Manguaba Lagoon (HSG4), Chl-a variability shows extremely high values in May associated to environmental factors of the rainy season (i.e., precipitation and flow).In this season, the higher Chl-a concentrations may be attributed to the higher inputs of nutrients and organic matter from the Sumaúma River and from the Marechal Deodoro city, as well as to lower turbidity and to the dominance of freshwater cyanobacteria species, particularly Microcystis aeruginosa, since salinity concentrations drop in this season.Microcystis species are well-known as strategist species, being capable to move vertically, enabling them to find better conditions of light in the water column, thus favoring their growth [49].Moreover, these species are also able to group together, forming colonies and thereby hindering herbivory by predators.

Conclusions
Although the studied lagoons show different optical characteristics, it was found that an unique bio-optical model based on MODIS bands might be applied to estimate Chl-a concentrations for both lagoons, which may be used for time-series reconstruction of Chl-a concentrations.The higher Chl-a values found in Manguaba lagoon better fit the upper part of the model, while the lower Chl-a concentrations with a high interference of turbidity in Mundaú lagoon better fit the lower part of the model.Thus, despite the high rate of siltation and shallow depths, the high turbidity and/or productivity of these lagoons allowed the use of remote sensing techniques to estimate Chl-a concentrations without any interference of the bottom sediment on the spectral reflectance, since the depth of the photic zone (light penetration zone) is lower than the water depth in both lagoons.
The NIR-Red band ratio model for retrieving Chl-a in MMELS used only two MODIS bands (i.e., Rrs(859) and Rrs(645)), which are also available in other MODIS products with a higher spatial resolution (e.g., MOD09Q1 and MYD09Q1 with 250-m spatial resolution) .Therefore, we recommend performing a new model calibration considering more refined spatial resolution products in future investigations.Furthermore, new field campaigns are strongly recommended to develop and validate a local atmospheric correction algorithm, to improve the validation of the radiometric model obtained in this study, as well as the calibration of new models considering the radiometric and spatial resolutions of the newly released sensors, such as Sentinel-2/MSI, Envisat/MERIS and Sentinel-3/OLCI satellites.In addition, we recommend the use of airborne sensors (e.g., using drones) that allow sampling spectral reflectance with higher spatial coverage and radiometric resolution in future investigations.
As the presence of clouds in tropical estuaries is very frequent, satellite sensors with poorer temporal resolution may be not suitable for a most extensive temporal assessment of Chl-a in these systems.Thus, the satisfactory achievements obtained in this study considering a monthly temporal scale may be attributed to (a) the monitoring frequency of the MODIS, since a high number of images were discarded due to low quality information; and to (b) the high HRT for both lagoons, which provide a quasi-steady state for short-term scales (i.e., from days to weeks).
The statistical techniques used in this study were useful tools for identifying properly the spatiotemporal distribution patterns of Chl-a in the MMELS and the key environmental factors associated to this variability.The multivariate statistical methods were also able to identify the main components of variability of Chl-a.In the transverse direction, Chl-a variability was associated to the HRT and nutrient loadings from lateral sources.However, the role of HRT on the Manguaba Lagoon Chl-a variability still needs additional investigation.Our findings showed that precipitation and river discharge promote a stronger effect on Chl-a variability in the north of Mundaú lagoon and in the south of Manguaba lagoon, while wind, solar radiation, and temperature influence this variability more strongly in the northern area of Manguaba lagoon and the center-eastern area of Mundaú lagoon.This analysis did not show clearly the effects of the tide on the spatiotemporal patterns of Chl-a, which may be related to the time scale adopted in the study.

Figure 2 .
Figure 2. Near infrared-red (NIR-Red) band ratio algorithm to estimate Chl-a in MMELS: (a) semiempirical potential model using calibration data.The dashed lines represent 95% confidence lines of the linear regression; and (b) observed versus retrieved values of Chl-a using data entirely independent of the calibration data.The dashed line is a linear regression of data points and the solid line represents the discrepancy ratio between observed and predicted values.

Figure 2 .
Figure 2. Near infrared-red (NIR-Red) band ratio algorithm to estimate Chl-a in MMELS: (a) semiempirical potential model using calibration data.The dashed lines represent 95% confidence lines of the linear regression; and (b) observed versus retrieved values of Chl-a using data entirely independent of the calibration data.The dashed line is a linear regression of data points and the solid line represents the discrepancy ratio between observed and predicted values.

of 17 9 Figure 3 .
Figure 3. Isolines of the eigenvector values of the co-variance matrix of the (a,d) first; (b,e) second and (c,f) third rotated principal component for Mundaú and Manguaba lagoons.

Figure 3 .
Figure 3. Isolines of the eigenvector values of the co-variance matrix of the (a,d) first; (b,e) second and (c,f) third rotated principal component for Mundaú and Manguaba lagoons.

Figure 4 .
Figure 4. Homogeneous Spatial Groups (HSG) identified by Cluster Analysis using the first four principal rotated components from PCA in the (a) Mundaú and (b) Manguaba lagoons.

Figure 5 .
Figure 5.Long-term (2002-2016) monthly average values of Chl-a concentration for (a) Mundaú and (b) Manguaba lagoons, considering different HSG obtained in the Cluster Analysis.

Figure 4 .
Figure 4. Homogeneous Spatial Groups (HSG) identified by Cluster Analysis using the first four principal rotated components from PCA in the (a) Mundaú and (b) Manguaba lagoons.

Figure 4 .
Figure 4. Homogeneous Spatial Groups (HSG) identified by Cluster Analysis using the first four principal rotated components from PCA in the (a) Mundaú and (b) Manguaba lagoons.

Figure 5 .
Figure 5.Long-term (2002-2016) monthly average values of Chl-a concentration for (a) Mundaú and (b) Manguaba lagoons, considering different HSG obtained in the Cluster Analysis.

Figure 5 .
Figure 5.Long-term (2002-2016) monthly average values of Chl-a concentration for (a) Mundaú and (b) Manguaba lagoons, considering different HSG obtained in the Cluster Analysis.

Figure 6 .
Figure 6.Canonical correspondence analysis (CCA) ordination diagram with spatial and temporal scores (points) and environmental variables (arrows) for (a) Mundaú and (b) Manguaba lagoons.Temporal scores were split into dry (o) and rainy (+) seasons, respectively.Spatial scores (•) were split by the HSG obtained from cluster analysis (see Figure4).

Table 2 .
Data of each field campaign and the number of collected water samples in MMELS.From this total, 26 and 22 sampling stations were used to calibrate and validate the NIR-Red band ratio model, respectively.

Table 3 .
Statistics of Chl-a and TSS concentration measured on the water samples collected in the MMELS.

Table 4 .
Explained variance for the first four rotated principal components for both lagoon in MMELS.

Table 5 .
Statistics of estimated Chl-a considering different HSG for each lagoon.

Table 5 .
Statistics of estimated Chl-a considering different HSG for each lagoon.

Table 5 .
Statistics of estimated Chl-a considering different HSG for each lagoon.