Assessing Chlorophyll a Spatiotemporal Patterns Combining In Situ Continuous Fluorometry Measurements and Landsat 8/OLI Data across the Barataria Basin (Louisiana, USA)

Abstract

The acquisition of reliable and accurate data to assess environmental changes over large spatial scales is one of the main limitations to determine the impact of eutrophication, and the effectiveness of management strategies in coastal systems. Here, we used a continuous in situ Chl-a fluorometry sensor and L8/OLI satellite data to develop an algorithm and map Chl-a spatial distribution to assess the impact of freshwater diversions and associated high nutrient loading rates in the Barataria Basin (BB) complex, a coastal system in the northern Gulf of Mexico. We collected water quality samples at 24 sampling stations and high-frequency continuous fluorometry in situ [Chl-a] data along a ~87 km transect from 2019–2020. Field [Chl-a] values were highly correlated (r = 0.86; p < 0.0001) with continuous in situ [Chl-a] fluorometry values. These continuous in situ [Chl-a] values were significantly related to a surface reflectance ratio ([B1 + B4]/B3) estimated using L8/OLI data (exponential model; R² = 0.46; RMSE = 4.8, p < 0.0001). The statistical model replicated [Chl-a] spatial patterns across the BB complex. This work shows the utility of high-frequency continuous Chl-a fluorometry sampling coupled with L8/OLI image analysis to increase the frequency and number of field data sets to assess water quality conditions at large spatial scales in highly dynamic deltaic regions.

1. Introduction

Estuaries are coastal ecosystems characterized by a high biological diversity and productivity across latitudes [1]. However, despite their ecological and economic importance, estuaries are increasingly altered by anthropogenic impacts that modify biogeochemical cycling and ecosystem services (e.g., fisheries, water quality; [2,3]). Currently, coastal eutrophication is a pervasive environmental problem that alters the complex functional and structural properties of estuaries across fertility and salinity gradients [4,5,6,7,8]. Understanding the relationship between eutrophication processes and estuarine functional attributes has become a priority as managers seek to mitigate and regulate eutrophication at different spatiotemporal scales.

One of the major limitations to improving water quality and mitigating negative impacts caused by coastal eutrophication is the lack of data and information to evaluate the frequency, duration, and spatial extent of nutrient enrichment impacts on ecosystem primary productivity. Although eutrophication assessment initiatives have been developed since the 1950’s in the USA, the lack of extensive and intensive spatiotemporal data coverage still limits the development of more comprehensive regional evaluations, particularly the long term impacts of nitrogen (N) and phosphorous (P) loadings [9]. Since the 1990’s the National Estuarine Eutrophication Assessment (NEEA; 1999, 2007) has determined, for instance, the scale, scope and origin of nutrient enrichment conditions in 144 estuaries in the conterminous U.S. [9,10]. Yet, in coastal regions there are still data gaps to fully understand the environmental effects of increasing eutrophication mainly due to urban and industrial development that have critically altered the hydrology and geomorphology of watershed systems in a changing climate [11,12,13].

Nutrient enrichment in estuarine systems triggers a progression of eutrophic responses [14,15]. The most apparent is the significant increase of Chlorophyll a concentration (hereafter Chl-a) and macroalgal blooms [14] that are directly correlated to high N and P loading rates into estuarine and coastal waters [16,17]. Because phytoplankton rapidly responds to nutrient availability (i.e., hourly/daily), Chl-a is considered a good proxy to determine the eutrophic state of estuarine waters [18]. Indeed, since Chl-a is the main source of phytoplankton’s in vivo fluorescence [19,20], this biophysical property has promoted the development of remote sensing tools that use optical sensors to indirectly evaluate in situ fluorometry and determine phytoplankton productivity and biomass at large spatial scales (i.e., hundreds of km²; [18,19,21]). Yet, assessing rapid changes in Chl-a at large spatial scales in both coastal and inland waters using remote sensing instruments is still a challenge.

One of the major issues in the spatial extrapolation and interpolation of field data is separating the contribution of in situ water column optical properties from the actual optical signals (i.e., in vivo fluorescence/absorbance) of phytoplankton cells (i.e., “Case 2” waters; [22,23]). Although in vivo fluorescence is well correlated to Chl-a concentrations, there are interfering compounds or processes that can significantly affect this relationship. It is known, for instance, that taxonomic composition, physiological/morphological mechanisms (e.g., light, trophic status, nutrient stress), and turbidity are factors affecting the accuracy of Chl-a estimates via remote sensing [24,25,26,27]. Furthermore, in shallow systems, bottom reflectance often interferes with the phytoplankton signal in the water column [28].

To reduce this uncertainty when assessing Chl-a values, spatiotemporal trends of in situ fluorescence or absorbance measurements are needed to identify the relative contribution of other environmental factors to explain the total variance when converting fluorometric values to Chl-a values. To date, Chl-a datasets available for the calibration and evaluation of Chl-a/fluorometric models are limited in frequency and in spatial extend [21]. In this study, we use high-resolution, spatially explicit in situ data to increase the accuracy of Chl-a mapping in the Barataria Basin (USA; hereafter BB complex). The study site is a dynamic trophic estuarine system where optical water column properties are largely determined by the Mississippi River discharge fluctuations throughout the year [29,30].

We compiled data sets from a continuous water quality sampling system (i.e., SWAMMP-R) and contemporaneous Landsat 8/OLI (L8/OLI) images to develop a Chl-a detection algorithm for time series of L8/OLI images. Given the OLI sensor configurations of Landsat 8 (i.e., Operational Land Imager) that collect data from nine shortwave bands, it is possible to pair and calibrate image data with continuous field in situ Chl-a fluorometry values collected in estuarine and coastal waters at the pixel level. Specifically, the objectives of this study were to (1) evaluate continuous in situ Chl-a fluorometry values acquired over a large spatial scale that were coupled with L8/OLI images acquisition to increase the accuracy of statistical models linking continuous fluorometric-based Chl-a values to L8/OLI surface reflectance and (2) use these models to discern Chl-a spatial variability along salinity gradients in the BB complex. We particularly focus on locations impacted by controlled freshwater diversions aiming to enhance sediment inputs for wetland restoration and can introduce nutrients that may promote algae blooms including toxic cyanobacteria species [31,32,33,34]. Further, freshwater diversions could potentially trigger complex responses over extended coastal areas. These include the effects of increasing flooding duration, changes in relative elevation and accretion rates in wetlands, and reduction in salinity due to freshwater inputs [31,32,33,34,35]. Because controlled freshwater diversions can also increase N and P loading rates, we also determine the potential impact on phytoplankton biomass at variable spatial scales defined by differences in fertility gradients [35,36]. We hypothesize that the impact of freshwater water diversions on the BB complex plankton biomass reflected in Chl-a values can be accurately and efficiently assessed using a combination of calibrated statistical models and spatially intensive sampling using the continuous flow-thru SWAMMP-R (Spatially explicit WAter Monitoring system to inform Management Priorities and Risks) multiprobe. Understanding the spatiotemporal patterns and the interactions among hydrology and primary production at a landscape level in the BB complex will provide critical data and information to determine and quantify the effects of freshwater diversions on this complex estuarine system undergoing extensive wetland and hydrological restoration since the 1990′s [33]. This assessment is urgently needed to gauge both the present BB complex eutrophication magnitude and future changes as freshwater diversions are implemented in brackish and saline regions within the next 10 years [32,35,37,38].

2. Materials and Methods

2.1. Area Description

The Barataria Basin (i.e., BB complex) is considered one of the most productive estuarine systems in the northern Gulf of Mexico (nGOM) [39,40,41,42]. It is located between the west bank of the Mississippi River (MR) and Bayou Lafourche in southeastern Louisiana, USA (range: latitude, 29°5′ to 30°8′ N; longitude, 89°21′ to 91°3′ W) (Figure 1). This area is part of the Mississippi River Delta Plain (MRDP) and characterized by high levels of primary production due to high nutrient loading rates from the MR and agriculture/urban areas in the upper basin. The MR is the longest river system in North America (7854 km; including the Missouri river) and it ranks sixth and seventh globally in terms of water (610 km³ yr⁻¹) and sediment (210 Mt yr⁻¹) discharge [43,44,45]. The BB complex has a total area of approximately 6600 km² with the highest elevations ranging from ~4–5 m along natural levees and the lowest slightly above sea level in wetland areas [46]. Vegetation communities are well-defined in the basin from north to south with cypress-tupelo swamp defined by Taxodium distichum and Nyssa aquatica and freshwater marsh characterized by Panicum hemitomon, Eleocharis spp., and Sagittaria lancifolia at the upper basin; an intermediate-brackish marsh dominated by Spartina patens and Distichilis spicata at the middle basin; and saltwater marsh habitat dominated by Spartina alterniflora at the lower basin [46,47,48]. The total land area in the BB complex comprises agriculture fields (765 km²), urban areas (430 km²) and extensive wetlands (3587 km²) [49]. The major water bodies in the BB complex include Lac des Allemands and Lake Cataouatche at the upper basin; Lake Salvador, Bayou Perot, and Little Lake in the middle basin; and Barataria Bay proper at the lower basin (Figure 1). Mean water depth in this region ranges from 1.5–2.0 m [50,51]. The water residence time is variable among water bodies: 50–70 days in the middle of Lac des Allemands; 100–120 days in Lake Cataouatche; 100 days in Lake Salvador; 8–15 days in Bayou Perot; 10–30 days in Little Lake; and 1–4 days in Barataria Bay [52,53].

Water quality properties across the BB complex are determined by daily and seasonal drivers including tidal fluctuations, atmospheric factors controlling meteorological events (i.e., cold fronts, storms, and tropical cyclones), river discharge and anthropogenic activities. The BB complex is a micro tidal system with mixed diurnal and semi-diurnal tides. The daily tidal amplitude average is 0.32 m at the coast and 0.03 m at the upper basin; this tidal range is significantly influenced by wind during the winter season when cold fronts occur at 3–7 days duration while tidal oscillation is the main driver during the summer season [49,54,55]. Cold fronts in coastal Louisiana are characterized by a strong wind shear, an abrupt drop in atmospheric pressure, and a fast reduction in air temperature [56]. Although cold fronts are not as severe as tropical cyclones, wind direction/velocity control water levels and circulation, sediment transport (e.g., fluxes, wetland accretion/erosion, sediment resuspension), and the exchange of organic matter and nutrients between the BB complex and the GOM [38,56]. Because of the BB complex location in subtropical latitudes, maritime tropical air masses maintain persistent warm and humid conditions along coastal Louisiana [38,57]. Thus, mean annual temperature ranges from 19 °C (January) to 32 °C (July) (mean: 20 °C) while the minimum average ranges from 7 °C (January) to 22 °C (July) [38,57].

The main freshwater inputs into the basin are rainwater, surface runoff from surrounding areas, the Gulf Intracoastal Waterway, and the Davis Pond Freshwater Diversion (DPFD) located in the upper basin with a maximum operational discharge of 250 m³ s⁻¹ [53,54]. Although levees along the Mississippi River have reduced the hydrological connectivity between the river and the adjacent wetlands, freshwater discharge from the Southwest Pass is advected from the mouth into the estuary through tidal inlets. This exchange promotes inputs of freshwater (c.a., ~200 m³ s⁻¹; [58]), salinity, and nutrient concentrations at the coastal boundary [59,60,61]. The Mississippi River discharge starts increasing in November, and peaks between March and June to begin declining throughout July. Given the differences in regional coastal sediment transport and deposition due to shoreline erosion and inlet expansion, the total water exchange varies across the lower Barataria Bay’s four main tidal inlets: Barataria Pass (66%), Quatre Bayou Pass (18%), Caminada Pass (13%), and Pass Abel (3%) [58,62] (Figure 1).

Historically, the BB complex has been impacted by a combination of anthropogenic (e.g., construction of levees, canal cutting, pond creating, nutrient enrichment) and natural (e.g., sea level rise, cold fronts, hurricanes) disturbances [63]. One of the major impacts occurred when the BB complex became isolated from the Mississippi River after the closure of Bayou Lafourche in 1902 and the construction of extensive levees for flood control in the 1930s and 1940s; these major hydrological alterations resulted in a decrease of freshwater and sediment inputs into the basin [51,64]. The lack of freshwater and sediment inputs, along with increasing relative sea level rise, has driven sustained wetland loss since 1932. Although the land loss rate has decreased from 1995 to 2016 due to natural (e.g., lack of major storms since 2008, reduction of tectonic activity, geomorphology) and anthropogenic (e.g., shift in rates of hydrocarbons production since 1974, and navigation activities) factors, the BB complex is still undergoing wetland loss [65]. Overall, this land loss is presently the highest in the Unites States (c.a., 25.9 km² yr⁻¹) [65,66,67,68,69]. In an attempt to mitigate wetland loss, the DPFD, in service since 2002 [35] was constructed to reconnect the BB complex to the Mississippi River and allow the reintroduction of river flow into the interdistributary basin. The main DPFD management objective is to reduce saltwater intrusion, improve commercial and recreational fisheries and prevent and mitigate wetland loss by sediment reintroduction to allow marsh accretion [33,67] (Figure 1). However, despite the potential benefits, some studies suggest that the DPFD diversion could produce negative effects impacting the water quality in receiving lakes near to the diversions, and even promoting eutrophication in the estuarine coastal areas [30,49,70].

2.2. Field and Laboratory Methods

2.2.1. In-Situ Continuous Sampling

Continuous in situ Chl-a data (hereafter continuous in situ Chl-a) was collected throughout the BB complex using the flow-thru SWAMMP-R multiprobe based on the Dataflow flow-through sample acquisition platform [71,72,73,74]. The system (integrated by a Eureka manta +35 water quality multiprobe, a Pro Oceanus CO₂-Pro CV sonde, and a YSI ProSolo ODO optical dissolved oxygen meter) consists of sensors connected in series to measure temperature (°C), specific conductance (mS cm⁻¹), pH, dissolved oxygen (DO, mg L⁻¹), salinity (ppt), pCO₂ (ppm), and in situ fluorometers for Chl-a (µg L⁻¹) and Chromophoric Dissolved Organic Matter (CDOMf; µg L⁻¹)_. The Chl-a sensor is a Cyclops-7 fluorometer (Turner Designs) that detects Chl-a by using a light excitation wavelength at 460 nm (blue range) and a 685 nm fluorescence wavelength detection; the sensor is tested and calibrated using a liquid dye standard solution (Rhodamine WT Dye 400 ppb; Turner Designs) before field deployments. The CDOMf sensor is a Cyclops-7 fluorometer (Turner Designs) that utilizes UV excitation for CDOM detection, encompassing a wide range of CDOM types in coastal waters [26,75]. The system is mounted on a vessel, and while navigating, it continuously records data at 5-second intervals by automatically withdrawing water from the water column at a range from ~0.5 to 1.0 m depth. Each reading is georeferenced using an integrated GPS unit [72].

The sampling routes encompass a variable number of upstream-downstream transects that are selected a priori based on the spatial location of salinity gradients (Figure 1). The freshwater endmember begins at the Bayou Verret channel (hereafter Gate) located at ~5 km southeast from the DPFD, adjacent to the I-90 highway and before the flood control structure gate (Figure 1, “Gate station”). The channel is hydrologically connected to Lake Cataouatche. In a southeastern direction, the sampling transects cross both lakes Cataouatche and Salvador and continue into Bayou Perot, Little Lake, and via Grand Bayou; the transect ends in Barataria Bay close to the Quatre Bayou Pass (Figure 1, “Quatre Bayou Pass station”). Depending on time, weather conditions and tidal stage, the transect in the lower Barataria Bay is extended into two east-west transverse transects partially covering the central and lower estuarine regions in the wider bay area (Figure 1). The transect sampling takes an average of 5.5 h to survey and is assumed to represent a “snapshot” of environmental conditions within a tidal cycle for a particular date; the length of the transect from upstream (“Gate station”) to downstream (“Quatre Bayou Pass station”) is ~87 km. The length of the transverse east-west transects is variable (above: ~33.5 km; below: ~32.3 km) (Figure 1). In this study, we analyzed Chl-a data obtained from the continuous in situ Chl-a fluorometer sensor during the historic Mississippi River flood pulse (2019–2020) on 9 July 2019, 13 September 2019, 14 December 2019, 22 February 2020, and 30 May 2020.

2.2.2. In-Situ Discrete Sampling

To evaluate the continuous in situ Chl-a fluorometer values and response accuracy, we used regression analysis between those data sets and Chl-a concentrations determined from discrete water samples (hereafter discrete in situ Chl-a) following the EPA-Method 445.0-2. Discrete water samples were collected at 0.5–1.0 m depth along the field transects at sampling stations (N = 23–24) in summer (9 July 2019), fall (13 September 2019), winter (14 December 2019) and spring (30 May 2020) to characterize spatiotemporal variations in water bodies throughout the BB complex (Figure 1). Water samples were collected in 1 L HDPE amber Nalgene bottles, which were pre-rinsed with distilled and deionized water at the laboratory and then with in situ water when the bottle was filled with water exiting the flow-thru system. Samples were immediately stored at 4 °C and transported to the laboratory for analysis in the same day of collection. Chl-a concentrations were measured using a 90% acetone extraction during 18 h [76]. Total suspended solids (TSS) was determined by filtering between 200–700 ml of each discrete water sample within 24 h after sample collection; water was filtered using a pre-weighed, ash-free filter (47 mm, 0.7 µm, Whatman GF/F), then this filter plus sediment sample was dried for 24–48 h until no change in weight for the determination of total weight by volume. Additional water quality analyses (data not shown but see pCO₂ and CDOM in supplemental Figure S1) included dissolved inorganic nitrogen (DIN: NH₄, NO₂ + NO₃), soluble reactive phosphorus (PO₄), total alkalinity (TAlk), dissolved inorganic carbon (DIC) and particulate organic carbon (POC).

2.2.3. Satellite Imagery Acquisition and Processing

Nearly cloud-free Level-1 L8/OLI satellite imagery were downloaded from the USGS website (https://earthexplorer.usgs.gov/). L8/OLI images are geometrically corrected and are projected in the UTM (Universal Transverse Mercator) map with the WGS 84 datum. We used ENVI 5.3 (Exelis Visual Information Solutions, Boulder, CO, USA) and ArcGIS 10.5 (ESRI, Redlands, CA, USA) for image processing and to derive surface reflectance data. Image processing included a radiometric calibration to convert digital numbers (DN) from the satellite images into radiance values per pixel (units: µW/(cm² × sr × nm) and to perform the required atmospheric correction. The following formula [77] was used to convert DN to radiance values:

L_λ = M_L × Q_cal + A_L

(1)

where L_λ is the spectral radiance (W/(m² × sr × μm; top of the atmosphere radiance), M_L is the radiance multiplicative scaling factor for the selected band (RADIANCE_MULT_BAND_n from the metadata), A_L is the radiance additive scaling factor for the band (RADIANCE_ADD_BAND_n from the metadata), and Q_cal is the pixel DN value.

Atmospheric effects were removed by using the FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubus) atmospheric correction module included in ENVI 5.3. This atmospheric correction accounts for factors that can affect the radiation detected by the satellite sensor (i.e., sky light, atmospheric sun glint, and radiance reflected by white caps; [77]) when retrieving surface reflectance (R_t) data. Previous studies have demonstrated that FLAASH is effective to determine the surface reflectance produced by water constituents [21,77,78,79]. The equation [77] to calculate surface reflectance (R_t) is:

R_t = 1 − ρe / A (L − B × (ρe/1−ρeS) − La)

(2)

where L is the spectral radiance at-sensor pixel, ρe is an average surface reflectance for the pixel and its surroundings, S is the atmospheric spherical albedo, A and B are coefficients that depend on atmospheric and geometric conditions but not on the surface, and La is the atmospheric back scattered radiance. The selected model parameters included both the atmospheric US-Standard and the maritime aerosol models to correct the L8/OLI images. Additionally, multispectral settings were readjusted for L8/OLI images, including the adjustment for open water areas based on the proximity to the water-dominated areas. The FLAASH output is the apparent reflectance data where pixel values ranging from 0 to 10,000 represent a corresponding range from 0–100% reflectance. In some cases, pixel values are outside this range as result of extremely dark/shadow areas or noise; thus, the Band Math Tool from the ENVI toolbox [80] was used to rescale the pixel’s values to a floating-point format (i.e., range 0–1).

2.2.4. Remote Sensing-based Algorithms to Estimate Chl-a Concentrations

Most of the remote sensing algorithms for Chl-a determination are based on the optical properties (i.e., sunlight absorption patterns) of algal pigment structure and quantity [21]. These empirical models were developed using regression analyses to correlate reflectance data to in situ Chl-a concentrations. Despite the simplicity of this approach, results from this statistical analysis are highly informative about water quality conditions when using Chl-a as indicator of eutrophication or levels of nutrient enrichment and potential nutrient sources. However, one of the criticisms (and limitations) is that these empirical algorithms are calibrated and validated based on a small number of in situ samples (c.a., ~20 samples) [21,81,82]. In this study, we used Level-1 Landsat 8/OLI images (path 22, rows 39 and 40) captured on 21 February 2020 and continuous in situ Chl-a fluorometry data collected approximately 20 h apart (22 February 2020; N = 1441 data observations) to increase the sample size (>1000 sampling points). Due to cloud coverage at the Barataria Bay coastal boundary, the transverse transects were not considered for the model’s calculation and validation (Figure S2). Similarly, due to geomorphological features (e.g., narrow man-made channel distribution, dredging), constructed structures (e.g., levees and gates), and hydrological variability (e.g., channel width and water depth) linked to major ecological differences (i.e., extensive floating mats of the water hyacinth (Eichhornia crassipes) and tree canopies), we excluded the Gate transect section for the model forecasting when compared to the main water bodies. Thus, for the algorithm development we used continuous in situ Chl-a data collected along salinity gradients from the most northern area of Lake Cataouatche (0.16 ± 0.01 ppt) to lower Barataria Bay (6.8 ± 5.3 ppt) (Figure S2).

L8/OLI bands were combined into different spectral bands followed by applying linear, polynomial, and exponential regression models to calculate a three bands ratio algorithm to characterize the relationship between the surface reflectance (R_t) data and the continuous in situ Chl-a (i.e., SWAMMP-R acquisition) concentrations. We used the continuous in situ Chl-a data and the best correlated R_t bands ratio combination to fit and test the model (Figure 2). After sequentially selecting the best correlated R_t bands ratio combination, the dataset was evaluated for noise and outliers due to sensor response and boat speed using studentized residuals criterion in JMP Pro 15.1.0 software [83]. This revised dataset (N = 1332) was then separated into two subsets containing half (N = 666) of the total observations for the corresponding model calibration and validation steps. To test the model predictive capabilities, we used regression analysis, analysis of variance (α = 0.05), Nash-Sutcliffe efficiency coefficient (NSE), and cumulative sum test (CUSUM). To evaluate the model predictive power, the NSE coefficient determines the residual variance compared to the observed data variance [84,85], while the CUSUM analysis helps to identify and interpret patterns in complex spatial and temporal series data [86]. Once the model validation was completed, we used the selected model equation to map Chl-a concentrations across the BB complex spatial domain (Figure 2).

3. Results

3.1. Discrete and Continuous Chl-a Data Statistical Summary

The average in situ Chl-a concentration (±SD) obtained from discrete water samples collected during the Mississippi River flood pulse (2019–2020) was 18.5 ± 15.2 SD µg L⁻¹. The lowest Chl-a value (3.97 µg L⁻¹) was measured at Bayou Perot during Winter and the highest (102.92 µg L⁻¹) along the lower Barataria Bay transect in Spring 2020 during the MR peak discharge (Figure 3A). No significant seasonal patterns were observed throughout the BB complex. However, the presence of a spatially discrete phytoplankton bloom was observed in the lower Barataria Bay during Spring 2020 near the Quatre Bayou Pass-station (Figure 3A). The average TSS concentration during the Mississippi River flood pulse (2019–2020) was 27.2 ± 20.5 SD mg L⁻¹; minimum values were observed in Lake Salvador (<1 mg L⁻¹) during summer 2019 and a maximum concentration in the east central Barataria Bay transect (79 mg L⁻¹) during fall 2019 (Figure 1 and Figure 3B).

There was a significant correlation coefficient (r = 0.86) and coefficient of determination (R² = 0.74; p < 0.0001) between the discrete in situ Chl-a values and the continuous in situ Chl-a (SWAMMP-R) values indicating that the continuous in situ Chl-a values directly measured through fluorometry techniques are robust estimations of Chl-a concentrations; thus, these values are suitable to interpolate values across the BB complex (Figure 4).

A regression analysis between TSS and Chl-a revealed seasonal patterns where Chl-a was a major constituent of particulate matter load during the summer (9 July 2019) and spring (30 May 2020) seasons (Figure 5).

During the 2019–2020 MR flood pulse the continuous in situ Chl-a data also showed consistent spatial patterns, where consistently higher Chl-a concentrations occurred at the upper basin (low salinity; 0.20 ± 0.03 SD ppt) and lower Chl-a values at the lower basin (higher salinity; 10.4 ± 5.9 SD ppt). Also, seasonal phytoplankton blooms were characterized during winter (14 December 2019) and spring (30 May 2020) when maximum continuous in situ Chl-a values were 82.4 µg L⁻¹ in Little Lake (mean: 17.8 ± 14.1 SD µg L⁻¹) and 52.8 µg L⁻¹ in Barataria Bay (mean: 14.2 ± 12.5 SD µg L⁻¹) (Figure 6).

3.2. Determining the Best Band Ratio and Algorithm Development Application

Linear, polynomial, and exponential regression analyses were performed using a total of 1441 observations of continuous in situ Chl-a collected on 22 February 2020 and different ratio combinations of R_t L8/OLI bands derived from the L8/OLI satellite image acquired on 21 February 2020 (Figure 2). The best relationship between these two variables was observed using the continuous in situ Chl-a values and the ratio between the R_t bands 1, 4 and 3 (i.e., R_t Ratio (B1 + B4)/B3;

Table 1. Linear, polynomial, and exponential band regression results between continuous in situ Chl-a and L8/OLI surface reflectance data (R_t). Modeling was performed using 1441 sampling points distributed throughout the BB complex along salinity gradients collected on 22 February 2020 (see methods). R², coefficient of determination; RMSE, root mean square error. All values are significant at p < 0.0001, except when p ≤ 0.0002 (*). Bold text and numbers indicate the selected model used for spatial analysis and mapping (see Figure 8).

Bands Ratios	Average Spectrum Wavelength (nm)	Linear		Quadratic		Cubic		Quartic		Exponential
		R2	RMSE	R2	RMSE	R2	RMSE	R2	RMSE	R2	RMSE
B3/B1	560/443	0.23	5.75	0.27	5.59	0.32	5.42	0.37	5.20	0.26	5.66
B5/B2	865/480	0.17	5.99	0.19	5.91	0.21	5.84	0.23	5.75	0.19	5.93
B5/B3	865/560	0.15	6.03	0.19	5.89	0.21	5.84	0.24	5.73	0.18	5.96
B5/B4	865/655	0.12	6.17	0.17	5.99	0.19	5.92	0.20	5.86	0.14	6.09
B4/B2	655/480	0.01	6.53 *	0.01	6.53 *	0.03	6.47	0.04	6.43	0.01	6.53
B3/B2	560/480	0.05	6.39	0.09	6.25	0.09	6.25	0.11	6.21	0.06	6.37
B2/B1	480/443	0.16	6.03	0.18	5.93	0.19	5.93	0.19	5.93	0.18	5.96
B4/B1	655/443	0.02	6.48	0.04	6.43	0.05	6.41	0.08	6.30	0.03	6.47
B4/B3	655/560	0.29	5.53	0.32	5.42	0.33	5.40	0.33	5.37	0.31	5.44
B5/B1	865/443	0.12	6.17	0.16	6.02	0.19	5.92	0.19	5.92	0.13	6.11
B2/B3	480/560	0.04	6.42	0.08	6.29	0.10	6.24	0.10	6.24	0.05	6.40
B1/B2	443/480	0.14	6.08	0.17	5.97	0.19	5.92	0.19	5.93	0.17	5.99
B1/B3	443/560	0.19	5.89	0.29	5.53	0.29	5.52	0.35	5.32	0.25	5.68
(B1 + B4)/B3	(443 + 655)/560	0.39	5.12	0.50	4.66	0.50	4.63	0.54	4.47	0.46	4.80

). Band 1 is centered at 443 nm corresponding to the first Chl-a strong absorption peak; band 4 is centered at 655 nm (low point of the second Chl-a absorption peak) while band 3 is centered at 562 nm within the green range corresponding to the minimum Chl-a absorption ([87]; Figure 2). Indeed, R_t bands 1, 4 and 3 are commonly used for spectral estimates of Chl-a as these bio-physical properties strongly shape its contribution to surface reflectance [81,87].

Prior to the model calibration and validation, the full dataset (N = 1441) was evaluated for noise and outliers due to sensor response and boat speed using a studentized residuals criterion in JMP Pro 15.0.0 software (SAS Institute, Cary, NC, USA) (Figure 2). After removing outliers, the dataset (N = 1332) was divided into two subsets for model calibration (N = 666) and validation (N = 666). Overall, the proportional relationship between these two variables was best characterized using an exponential model between the continuous in situ Chl-a values and the ratio (R_t [B1 + B4]/B3), expressed as:

[Chl-a] = 14027.14 × e^{(−3.36×[(B1 + B4)/B3])}

(3)

Although more complex polynomial regression models were significant (Table 1), they were biased (data not shown) towards the lower and middle Chl-a values (~50% of the total observations; Figure 7). This bias was compounded by the spectral signature of the Gate transect data (N = 68 observations), which showed major differences in water optical properties as result of hydrology and ecological processes when compared with the rest of the water bodies in the BB complex (Figure S1). Thus, these transect values, representing ~5% of the total number of continuous in situ Chl-a values (N = 1441) collected on 22 February 2020 were excluded in the regression analysis (see Section 4).

3.3. Algorithm Validation

The exponential model was statistically significant and able to predict Chl-a concentrations along the BB complex transect sections (r = 0.77). A total of 666 observations were used to validate the model predictive power (R² = 0.61; RMSE = 3.44; NSE = 0.61; p < 0.0001; Figure 8).

The Chl-a concentration was not significantly different (α = 0.05) between the continuous in situ Chl-a and the models predicted Chl-a for the entire BB complex (F(1, 1330) = 0.02, p = 0.89; Figure 9).

The selected model was able to replicate the frequency distribution of the validation dataset, where ~50% of the Chl-a observations are between 5–15 µg L⁻¹, and its predicted values were similar to the continuous in situ Chl-a values across the entire salinity gradient (Figure 10).

Chl-a concentrations were highly variable along the BB complex transects on 22nd February 2020; the highest average continuous in situ Chl-a value was recorded in Lake Cataouatche (26.60 ± 5.17 µg L⁻¹; Figure 6). Although the CUSUM test did not completely replicate the Chl-a patterns along salinity gradients, it confirmed the distinctive increase of Chl-a concentrations and presence of outliers in Lake Cataouatche (Figure 10). These spatial patterns are controlled by hydrological differences, freshwater inputs, water residence time, tidal exchange [53] among well-defined subsystems within the BB complex [52,53]. These hydrological differences control nutrient availability and transport among BB complex sub-basins, promoting high Chl-a concentrations (≥60 µg L⁻¹) in Lac des Allemands and Lake Cataouatche where the Chl-a maximum zones can be detected as plumes (Figure 11A–C).

4. Discussion

Our study demonstrated the value of using high frequency (SWAMMP-R) in situ Chl-a fluorometry-based measurements paired with Landsat 8/OLI to characterize spatially-explicit Chl-a distribution at large spatial scales throughout the BB complex. The regression model using the R_t ratio ([B1 + B4]/B3) was able to characterize Chl-a spatial patterns in water bodies with variable environmental conditions (i.e., depth, hydrology, and freshwater inputs). CUSUM analysis of continuous in situ Chl-a revealed Chl-a trends along salinity gradients with the highest average value observed in Lake Cataouatche (26.0 ± 5.3 µg L⁻¹) in the upper basin where mean salinity was 0.16 ± 0.01 ppt. The lower continuous in situ Chl-a values were registered throughout the middle basin in Little Lake (9.3 ± 1.6 µg L⁻¹; mean salinity: 0.20 ± 0.06 ppt) and the lower basin in Barataria Bay (12.8 ± 2.7 µg L⁻¹; mean salinity: 6.8 ± 5.3 ppt).

The exponential model generated Chl-a concentration patterns comparable to the continuous in situ Chl-a across the BB complex, predicting higher Chl-a concentrations in Lake Cataouatche (22.8 ± 6.6 µg L⁻¹), followed by lower values throughout the middle basin at Lake Salvador (11.6 ± 1.8 µg L⁻¹), Bayou Perot (11.8 ± 2.3 µg L⁻¹), and Little Lake (9.7 ± 1.9 µg L⁻¹) and the lower basin at Barataria Bay (10.9 ± 2.1 µg L⁻¹). The selected model also confirms that the highest Chl-a concentrations occur in the upper basin at Lac des Allemands (Figure 11B), which has been reported as a highly productive area during the last ~40 years due to major nutrient loading rates resulting from agriculture practices and urban development [29,30,47,49,88].

Due to the low correlation (R² = 0.02; p = 0.15) between TSS and Chl-a concentration (N = 84 samples), we propose that the BB complex became an optically complex system [Case 2 water; [81]] during the MR flood in 2019–2020. When this relationship was analyzed by season, we detected a positive correlation between TSS and Chl-a in spring and summer, but a negative trend in winter (Figure 5). These results suggest that Chl-a becomes the major constituent driving water color during summer and spring when nutrient concentrations, temperature and solar irradiance are high coinciding with the peak water column productivity [30,49,89,90]. Moreover, observations from both the discrete sample assays and field high-frequency fluorometry allow us to differentiate the optical properties among different water masses. We noted the most conspicuous difference in the water column along the Gate transect section, where mats of floating vegetation covered the channel surface due to the interaction of high nutrient concentrations and wind direction/velocity causing rapid plant growth and vegetation biomass accumulation. One significant difference between water masses among the Gate channel and the rest of the transects across the BB complex was the high pCO₂ concentrations in the water column (Figure S1). During the 22 February 2020 field sampling, high pCO₂ values in the Gate channel area ranged from 1459.6 to 4095.6 ppm, in contrast to the lower range (423.4–520.4) measured in the Barataria Bay transect (Figure S1). This high pCO₂ concentration reflects the high respiratory activity occurring in the narrow channel due to the decomposition of live and dead plant material (organic matter) accumulated in the channel bottom and the water column. Water oxygen concentration >3 mg L⁻¹ facilitates organic matter decomposition, releasing organic derived substances to the water column that influence the optical characteristics (i.e., surface reflectance) captured by the Landsat 8/OLI satellite in the Gate transect (Figure 7). This type of channel (i.e., narrow and running across wetland areas) is common throughout the BB complex, not only as result of natural sediment deposition patterns, but also as result of dredging for oil and gas exploration and navigation activities, particularly in the last 50 years [66].

Our analysis showed that within each season and during river flooding conditions, the average continuous in situ Chl-a concentrations are higher at the upper basin in the Gate transect (Figure 6; 26.8 ± 14.1 µg L⁻¹) and Lake Cataouatche (Figure 6; 21.8 ± 11.5 µg L⁻¹) while salinity increases in the middle of Lake Salvador (Figure 6; 7.7 ± 3.8 µg L⁻¹) and the lower Barataria Bay basins (Figure 6; 13.7 ± 6.2 µg L⁻¹). We attribute this spatial pattern to both the direct effect of the MR flood pulse, when high concentrations of inorganic nutrients enter through the DPFD into the northern region of Lake Cataouatche and the higher water residence time (>100 days; [53]) in these systems that promote high phytoplankton productivity [91]. This production pattern is similar in Lac des Allemands where the residence time is also high (50–70 days; [52,53]), and where there is also a constant input of inorganic nutrients due to runoff from surrounding sugar cane farms [49,88]. The local plume distributions of both TSS and Chl-a are clearly demarcated in both lakes as freshwater input enters in pulses, highlighting the relative importance of nutrient inputs for controlling water column primary productivity at local scales (Figure 11A–C). On the marine end-member of the salinity gradient downstream, the MR plume also enters the lower basin through tidal inlets as showed in previous studies [30,60]. This direct riverine input enters into the lower Barataria Bay during the peak of the MR discharge in combination with water currents running parallel to the adjacent continental shelf, causing an increase in Chl-a in Barataria Bay as underscored by water with high inorganic nutrient concentrations during the period from 1992 to 2016 (e.g., NO₃ + NO₂, 12.0 ± 15.2 µM) [49,92].

The interaction among the MS river discharge seasonality, differences in water residence time, and biological processes across the BB complex can influence Chl-a concentration and phytoplankton productivity. This interaction is apparent in the upper basin when nutrient loading is low during the summer [93]. Water residence time increases during the summer when the Mississippi River is at its lowest level reducing nutrient inputs into the estuary [52]. Sampling of inorganic nutrients along the transect (data not shown) indicates that water bodies in the BB complex might be nutrient limited during the summer when temperature increases, and peak phytoplankton productivity occurs. Similarly, low nitrogen concentrations have been associated with denitrification rates in benthic sediments thus promoting nitrogen limitation during the summer when denitrification rates are the highest [94,95]. Low Chl-a concentration in the upper and middle basins can also be attributed to consumption by benthic organisms where the presence of Rangia clam populations (Rangia cuneata) is extensive. Shallow water column depth and high water residence time facilitates the availability of organic matter, including phytoplankton, to benthic communities [52,96].

The observed salinity and predicted Chl-a gradients revealed well defined subsystems within the BB complex as underscored by other studies (i.e., [53]). These subsystem boundaries are supported by the high Chl-a spatial variability we documented within Lake Cataouatche when compared with Chl-a values in the middle and lower basins. Although all BB complex sections captured by our SWAMMP-R transects represent a hydrological continuum, these Chl-a patterns reveal the spatially and temporally local effects of freshwater inputs, residence time, geomorphology, and weather (i.e., precipitation, cold fronts, and major storms) as reflected by the year-round patchiness in Chl-a spatial distribution in Lake Cataouatche (Figure 6 and Figure 11B1,C1). This irregular and patchy Chl-a distribution is difficult to capture using only discrete sampling techniques along transects as shown by our results and previous studies [49,88].

The ecological responses at local/regional scales to extreme weather and climate patterns (e.g., precipitation, wind velocity and strength, tropical storms) in estuaries and lagoons along coastal Louisiana are not well understood, particularly under current climate variability. For instance, in August 2016, a stationary low pressure system caused more than 50 cm of cumulative precipitation in a three-day period affecting coastal areas in South Louisiana including record flooding levels inland [97]. This event was associated to the lowest average annual salinity recorded in 2016 during a span of 22 years (1994–2016) in the BB complex (Figure S3; [49]. Recently, the 2020 hurricane season was a record for coastal Louisiana with four major storm landfalls in one season (i.e., tropical storms and hurricanes). For example, hurricane Delta, (category-2) at the time of landfall in Creole, LA and ~335 km from the BB complex, caused a storm surge (1–1.5 m) triggering ocean water intrusion into the lower BB complex. A salinity survey using the SWAMMP-R system and performed on 18 October 2020, nine days after Delta’s landfall showed saltwater intrusion up to Lake Cataouatche where salinity increased three times the annual average when compared to measurements in previous years (Figure S4A,B); this increase in salinity underscores the regulatory pulsing control of these climate events on water quality in the BB complex.

The development of robust quantitative relationships between in situ measurements of environmental and ecological parameters and remotely sensed earth surface observations is especially useful for very dynamic ecosystems and landscapes, places characterized by high variation across broad ranges of temporal and spatial scales. Our exponential model (Equation (3)) reflects this complex environmental signature across the distinctive water bodies that form part of the BB complex. This strategy for the developing of empirical models has been widely used for specific conditions, regions and even for large scale ocean color descriptions including hypertrophic water systems (West Lake-Vietnam; [21]), highly turbid lakes and coastal waters (Dnieper river, Changjiang estuary; [98,99]), freshwater lakes (Great Lakes; [100]), coastal water areas (Florida shelf; [101]), and open ocean waters [102,103].

In open ocean waters, satellite images have been used to retrieve water quality data since the late 1970′s when the Coastal Zone Color Scanner satellite (CZCS) was launched in 1978 [81]. The main limitation for the development of robust algorithms, particularly for coastal areas (e.g., [92]), is the lack or reduced number of in situ observations to train and validate these quantitative models. While previous studies have demonstrated the potential use of Landsat 8/OLI hyperspectral data to determine Chl-a, few direct field measurements (N < 25) are available to validate satellite algorithms in large study areas (e.g., 97 km²) [104]. We addressed this limitation by increasing the number of in situ Chl-a values (N = 1332) used in model calibration and validation by acquiring continuous in situ Chl-a field fluorometry measurements across the BB complex (open water area ~1903 km²; [52]).

The development of new technologies and improvement of existing instrumentation open wide possibilities for monitoring water quality at large spatial and temporal scales, while also allowing the characterization of hydrological and biochemical processes at local scales. We anticipate enhanced understanding of critical ecosystem processes like primary productivity and biomass (i.e., Chl-a), that are driven by pulsing events such as cold fronts, river discharge, and storm surges after hurricane passage on primary productivity and biomass (i.e., Chl-a). In this study we showed that linking Landsat 8/OLI sensor observations, with a high frequency environmental sampling approach is a powerful strategy to develop effective and robust algorithms for the assessment of both water quality and optical properties. This is especially useful in spatially large-scale and highly dynamic estuarine systems (i.e., “Case 2” waters; [22,23,92]) like those found in the Mississippi River delta plain, where intricate physical and ecological dynamics cannot be captured by discrete sampling alone.

5. Conclusions

Our results demonstrated that the use of Landsat 8/OLI images provided reliable results to characterize the spatiotemporal Chl-a patterns in the Barataria Bay complex, a hydrologically delta-dominated coastal ecosystem. An exponential model was validated to evaluate Chl-a seasonal and spatial distribution using an intense acquisition of fluorometry measurements and a reflectance R_t ratio ([B1 + B4]/B3) from Landsat 8/OLI imagery/data. These results are supported by contemporary and historical data, especially the effects of local environmental variables controlling Chl-a spatial distribution and magnitude. Model interpolation and predictions showed that higher phytoplankton concentrations (i.e., Chl-a) are extensive and persistent in the upper basin, partially due to significant nutrient loadings from extensive agricultural lands and from the MR discharge through the seasonal and pulsing DPFD delivery. Our findings complement previous studies in which salinity is used as a proxy to characterize well-defined subsystems within the BB complex, and where despite the hydrological connectivity, distinctive biochemical processes, and productivity patterns define each subsystem. To further expand the validation of Chl-a spatiotemporal patterns using our proposed algorithm, we recommend the analysis of historical and new Landsat 8/OLI data to assess past, current, and future ecological and hydrodynamic changes. This is particularly pressing from a management perspective, since this deltaic estuarine system is undergoing large-scale hydrological restoration as it experiences increasing rates of relative sea level rise and wetland loss as a result of human disturbances, hurricane impacts, and climate change [105,106].