Atmospheric Correction of Satellite Optical Imagery over the R í o de la Plata Highly Turbid Waters Using a SWIR-Based Principal Component Decomposition Technique

: Estimating water reﬂectance accurately from satellite optical data requires implementing an accurate atmospheric correction (AC) scheme, a particularly challenging task over optically complex water bodies, where the signal that comes from the water prevents using the near-infrared (NIR) bands to separate the perturbing atmospheric signal. In the present work, we propose a new AC scheme specially designed for the R í o de la Plata—a funnel-shaped estuary in the Argentine– Uruguayan border—highly scattering turbid waters. This new AC scheme uses far shortwave infrared (SWIR) bands but unlike previous algorithms relates the atmospheric signal in the SWIR to the signal in the near-infrared (NIR) and visible (VIS) bands based on the decomposition into principal components of the atmospheric signal. We describe the theoretical basis of the algorithm, analyze the spectral features of the simulated principal components, theoretically address the impact of noise on the results, and perform match-ups exercises using in situ measurements and Moderate Resolution Imaging Spectrometer (MODIS) and Visible Infrared Imaging Radiometer Suite (VIIRS) imagery over the region. Plausible water reﬂectance retrievals were obtained in the NIR and VIS bands from both simulations and match-ups using ﬁeld data—with better performance (i.e., lowest errors and offsets, and slopes closest to 1) compared to existing AC schemes implemented in the NASA Data Analysis Software (SeaDAS). Moreover, retrievals over images in the VIS and NIR bands showed low noise, and the correlation was low between aerosol and water reﬂectance spatial ﬁelds.


Introduction
In atmospheric correction (AC) of optically complex waters such as the Río de la Plata (RdP), the black water assumption in the near-infrared (NIR) bands on which traditional AC schemes rely [1,2] is usually invalid due to the high backscattering produced by the high concentration of suspended particulate matter (SPM) present in the water [3][4][5]. Furthermore, over RdP, the iterative schemes based on the NIR bands-such as [6,7]-usually do not provide good estimates in Moderate Resolution Imaging Spectrometer (MODIS) images due to either the saturation of NIR bands or limitations in the bio-optical model used to estimate the water reflectance in the NIR and visible (VIS) bands. Although other strategies based on the black water assumption in the shortwave infrared (SWIR) bands have already been proposed [8,9], in this work we explore the possibility of implementing this strategy but with a different extrapolative scheme based on the principal component analysis (PCA) of a set of radiative transfer simulations in which only the atmosphere and the air-water interface effects are considered (i.e., where water reflectance, ρ w , is equal to zero). The main difference with existing AC schemes is that the quantities used for

Algorithm Description
The PCA-SWIR AC starts from assuming a negligible contribution of the water-leaving radiance over the total radiance budget received by the satellite sensor at top-of-atmosphere (TOA) in the SWIR bands and assuming a high correlation of the aerosol signal between the SWIR bands and the VIS/NIR bands. Consider the following expression for the Rayleigh-corrected (RC) reflectance at top-of-atmosphere (TOA), ρ RC , [1]: where ρ a , Tρ g, and tρ w are the aerosol, glint, and water reflectance contributions at TOA, and T and t are the direct and total (direct + diffuse) transmittance factors, respectively. In this equation, ρ RC is defined as πL RC /[F 0 cos(θ s )] where L RC is the respective Rayleighcorrected radiance, F 0 the extra-terrestrial solar irradiance, and θ s the solar zenith angle, and ρ w is defined as πL w /E d where L w is the water-leaving radiance and E d the downwelling solar irradiance just above the surface. Assuming that the water reflectance in the SWIR is negligible (the so-called black pixel assumption) and selecting regions that are not affected by direct sun glint, the RC reflectance in the SWIR and VIS-NIR bands (subscript "VN") can be expressed as: ρ RC (λ VN ) ≈ ρ a (λ VN )+t(λ VN )ρ w (λ VN ).
We propose to model the aerosol term as a combination of principal components of a set of radiative transfer simulations in which the water reflectance is set to zero (ρ w = 0) (see Section 2.2). The principal components correspond as usual to an orthogonal basis Remote Sens. 2021, 13, 1050 3 of 26 of eigenvectors of the variance-covariance matrix of the set of simulated RC reflectance values. Assuming that most of the explained variance is expected to be contained in the first two or three principal components, the aerosol reflectance in the VIS/NIR bands, ρ a (λ VN ), can be obtained as: ρ a (λ VN ) ≈ a 1 e PCA 1,VN +a 2 e PCA 2,VN +ρ a (λ VN ), (4) where e PCA j,VN represents the VN component of the j-th principal component of the set of simulations, a 1 and a 2 represent the projections of the first N components, and <·> represents the mean value of the aforementioned ensemble of simulations.
Notice that the total number of eigenvectors is N + 1, i.e., the number of SWIR bands (N = 2 or N = 3 depending on the scheme) plus the VIS/NIR band to be corrected, meaning that the last principal component is assumed to have negligible contribution to the aerosol signal in Equation (4). Given that the N first components explain >95% of the variance in the considered set of simulations (see Section 4.1.1), this assumption is considered as reasonable. The projections of the N first principal components (hereafter N = 2 for brevity), a 1 and a 2 , can then be obtained by inverting the linear system composed of Equation (2) specified at the two SWIR bands and assuming the same decomposition as in Equation (4) The N × N matrix (2 × 2 in this example) that appears on the left hand side of Equation (5) is a subset of the variance-covariance matrix of size (N + 1) × (N + 1) and is the fundamental matrix of the scheme, herein M, being M j,λ = e PCA j,λ (j = 1, . . . ,N, and λ = SWIR1, . . . ,SWIR-N). Any error in the right-hand side of Equation (5) (sensor noise, calibration error, or error introduced in the obtention of ρ RC ), is potentially amplified by the linear inversion and the bound to such error amplification is measured by the condition number of M: where M = max x =0 M.x 2 x 2 , being x an N-tuple and · 2 the conventional Euclidean norm.
where m = sec(θ s ) + sec(θ v ) is the air mass factor and where an average anisotropy factor of 2/3 is assumed for aerosols. To estimate the transmittance factor, the following aerosol optical thickness τ a (λ) was used: This expression is a power law where the Angstrom exponent is fixed at 1, as in Steinmetz et al., 2011 [11], and the value at τ a (500 nm) = 0.06 corresponds to the mode of the set of values measured at the CEILAP-BA AERONET station located in Buenos Aires (approximately 3 km from RdP northwest coast). Naturally, the anisotropy factor and the aerosol optical thickness are not fixed in nature, but this simplified expression for the transmittance is preferred in this work given that the focus is on the PCA extrapolation procedure and the fact that the approximations used to obtain the transmittance factor typically have a low impact over the total error budget of the AC in comparison to other Remote Sens. 2021, 13, 1050 4 of 26 sources of error. It must be reinforced that these constraints are posed only over the transmittance factor (a second order correction factor) and not over the set of simulations from which e PCA j,λ is derived. Figure 1 schematizes the steps performed to calibrate the algorithm and to apply it to both images and simulated data. (approximately 3 km from RdP northwest coast). Naturally, the anisotropy factor and the aerosol optical thickness are not fixed in nature, but this simplified expression for the transmittance is preferred in this work given that the focus is on the PCA extrapolation procedure and the fact that the approximations used to obtain the transmittance factor typically have a low impact over the total error budget of the AC in comparison to other sources of error. It must be reinforced that these constraints are posed only over the transmittance factor (a second order correction factor) and not over the set of simulations from which e , is derived. Figure 1 schematizes the steps performed to calibrate the algorithm and to apply it to both images and simulated data. In red boxes are the sources/inputs: satellite data (Moderate Resolution Imaging Spectrometer (MODIS)/ Visible Infrared Imaging Radiometer Suite (VIIRS)), simulated data using the successive orders of scattering radiative transfer code (SOS) and field measurements (Analytic Spectral Device (ASD)/"Tri" Optical Sensors (TriOS)). In blue boxes the outputs: the principal component analysis (PCA) eigenvectors and the PCA-shortwave infrared (SWIR)-derived satellite and simulated water reflectance. Gray solid arrows symbolize each of the intermediate steps in each process, and the dashed black arrows show the comparisons that were performed to test the algorithm.

Radiative Transfer Simulations
To simulate the top-of-atmosphere (TOA) reflectance in the spectral bands of the sensors listed in Table 1, the Centre National d'Études Spatiales/Successive Orders of Scattering (CNES-SOS) radiation transfer code v5.0 [19,20] was used to produce two sets of simulations: Set 0: This set was used to calibrate the principal components of the atmospheric signal, i.e., the eigenvectors of the variance-covariance matrix of the set in the VIS/NIR/SWIR bands. It was composed of CNES-SOS simulations with null water reflectance (ρw = 0). Set W: This set was used to theoretically evaluate the performance of the AC in the VIS/NIR bands. It was composed of CNES-SOS simulations where a subset of 22 water reflectance spectra measured with an Analytic Spectral Device (ASD) was used as input instead of the black water condition of Set 0 ( Figure 2). The inclusion of water reflectance to the CNES-SOS code is done through the addition of a Lambertian (isotropic non-polarized) term in the surface Mueller matrix, which also contains the Fresnel reflectance term (for further detail see [19]). This subset of 22 measurements was selected in order to span a large range of measured water reflectance in the RdP with a pseudo-uniform step, i.e., trying to cover in the most uniform possible way all the range of water reflectance values Figure 1. Global scheme of the methodology developed in this study. In red boxes are the sources/inputs: satellite data (Moderate Resolution Imaging Spectrometer (MODIS)/ Visible Infrared Imaging Radiometer Suite (VIIRS)), simulated data using the successive orders of scattering radiative transfer code (SOS) and field measurements (Analytic Spectral Device (ASD)/"Tri" Optical Sensors (TriOS)). In blue boxes the outputs: the principal component analysis (PCA) eigenvectors and the PCA-shortwave infrared (SWIR)-derived satellite and simulated water reflectance. Gray solid arrows symbolize each of the intermediate steps in each process, and the dashed black arrows show the comparisons that were performed to test the algorithm.

Radiative Transfer Simulations
To simulate the top-of-atmosphere (TOA) reflectance in the spectral bands of the sensors listed in Table 1, the Centre National d'Études Spatiales/Successive Orders of Scattering (CNES-SOS) radiation transfer code v5.0 [19,20] was used to produce two sets of simulations: Table 1. Spectral bands of the sensors used to theoretically test the SWIR-PCA algorithm. The bands in which the black water assumption is valid are marked in green. In yellow, the bands in which this assumption is affected only with very high sediment concentrations. Set 0: This set was used to calibrate the principal components of the atmospheric signal, i.e., the eigenvectors of the variance-covariance matrix of the set in the VIS/NIR/SWIR bands. It was composed of CNES-SOS simulations with null water reflectance (ρ w = 0). Set W: This set was used to theoretically evaluate the performance of the AC in the VIS/NIR bands. It was composed of CNES-SOS simulations where a subset of 22 water reflectance spectra measured with an Analytic Spectral Device (ASD) was used as input instead of the black water condition of Set 0 ( Figure 2). The inclusion of water reflectance to the CNES-SOS code is done through the addition of a Lambertian (isotropic non-polarized) term in the surface Mueller matrix, which also contains the Fresnel reflectance term (for further detail see [19]). This subset of 22 measurements was selected in order to span a large range of measured water reflectance in the RdP with a pseudo-uniform step, i.e., trying to cover in the most uniform possible way all the range of water reflectance values that were measured in the region of interest in order to avoid sampling biases in the retrievals. The reflectance in each band of interest (Table 1) was calculated using the spectral response functions (SRFs) that are available from the NASA Ocean Biology Processing Group (OBPG) webpage [21] in the cases of MODIS and VIIRS, and using square-shaped SRFs defined by the centers and bandwidths reported by CONAE in the case of the SABIA-Mar sensor (Dr. Carolina Tauro Table 1).   Table 1).

Band
The remaining parameters used as input to CNES-SOS are detailed in Table 2. The Rayleigh scattering due to air molecules was parameterized following the optical thicknesses and molecular depolarization ratio reported by Bodhaine et al., 1999 [22], using an e-folding molecular height scale of 8 km, i.e., assuming an isothermal-barotropic atmosphere. The aerosol physical properties, characterized by their granulometry and complex refractive index, were parameterized assigning linear combinations of the continental, maritime, and urban scenarios proposed by the World Meteorological Organization [23]. Assuming no prior information regarding aerosol type, we considered equally all the 15 combinations produced by combining these three types of scenarios in quarters, e.g., maritime 25%, continental 25%, urban 50%, or maritime 25%, continental 75%, urban 0%. An e-folding height scale of 2 km was assumed for aerosols. Their concentrations were parameterized by the aerosol optical thickness at surface and at a reference wavelength of Remote Sens. 2021, 13, 1050 6 of 26 500 nm, τ a (500 nm), log-normally distributed with a mode of 0.06, determined from sunphotometric measurements made at the CEILAP-BA AERONET station in the 1999-2013 period [24]. The range of surface wind speeds (used by the SOS code to calculate surface roughness) was selected based on the wind measurements at the Aeroparque Meteorological Station-located just by the RdP river in the city of Buenos Aires-for the period 1976-2014. The values of 0, 2, 4, 8, and 16 m/s were chosen, geometrically distributed, given that the most probable wind speeds were within the range of 4 to 6 m/s, and winds greater than 14 m/s were recorded in much less than 1% of the measurements. In the simulations, the effect of whitecaps was not considered, given their negligible contribution in the highly reflective waters of interest in this study. Table 2. Atmospheric, surface, geometric, and optical parameters used as input to the sets of simulations using Centre National d'Études Spatiales/Successive Orders of Scattering (CNES-SOS). Rayleigh correction was performed over both sets of TOA simulations (Sets 0 and W) before computing the associated eigenvectors to Set 0 and applying the PCA scheme to Set W. To achieve this, the code was run twice for each subset to compute the TOA reflectance (ρ SOS TOA ) with two conditions: (i) using the input values specified by Table 2, and (ii) assuming a black marine environment (ρ w = 0) and no aerosol content (τ a = 0). The results from (ii) were subtracted from those from (i) to yield the RC signal, both in Sets 0 and W:

CNES-SOS Parameter
where Surf stands for the air-water interface. To prevent an overrepresentation of sun glint on the PCA eigenvectors, we discarded the subset of sun-observation geometries and wind speed combinations-in both Sets 0 and W-where the reflectance due to sun glint exceeded the threshold of 0.005. To perform this pre-selection, the sun glint reflectance, ρ g , was computed separately by generating a set of simulations with null water-leaving radiance reflectance and no atmosphere (ρ w = 0 and τ r = τ a = 0). The PCA-SWIR AC over the simulated RC reflectance ensemble (Set W) was applied using the four possible SWIR band combinations (see Table 1), i.e., PCA-SWIR12 when using SWIR bands 1 and 2, and the same for PCA-SWIR13, PCA-SWIR23, and PCA-SWIR123.

Other AC Algorithms and Match-up Procedure
We applied the PCA-SWIR scheme and other existing ACs on MODIS (onboard Aqua and Terra) and VIIRS (onboard Suomi-NPP and NOAA20 platforms) imagery over RdP for the dates when field measurements were collected within the period 2012-2020 (Section 3.1). Even though we tested all the four possible PCA-SWIR schemes over the imagery, in the current paper we only show the global intercomparison between existing ACs and PCA-SWIR13, since this is the only option that does not make use of the band SWIR2, which is no longer operational in MODIS/Aqua. Apart from PCA-SWIR13, the following ACs implemented in the SeaDAS/L2Gen ocean color Level2 processor [21] were also applied and considered in the match-up exercise (also summarized in Table 3): • NIR multi-scattering extrapolative approach (GW94-NIR): This AC was proposed by Gordon and Wang 1994 [1] and computes the aerosol signal accounting for multiple scattering effects and assuming black water in the NIR. It was the standard AC implemented for Sea-viewing Wide Field-of-View Sensor (SeaWIFS) imagery and is suitable for open waters; • Iterative scheme in the NIR (ITER-NIR): This AC is based on an iterative procedure to subtract the non-zero contribution of the water to the NIR signal based on a convergence strategy and on a semi-empirical model for the water reflectance in the NIR and VIS. It is described in [6,7]. This AC has been shown to underestimate water reflectance and even fail to produce any retrieval in RdP's maximum turbidity front [5,17], but it is considered here because it is the standard AC implemented by L2Gen in SeaDAS; • SWIR multi-scattering extrapolative approach (GW94-SWIR): This AC computes the aerosol signal accounting for multiple scattering effects and assuming black water in two specified bands in the SWIR. It is similar to the Gordon and Wang 1994 [1] approach, except it is adapted to SWIR bands [4,5,8]; • Rayleigh correction scheme (RC): This AC only removes the effect of Rayleigh scattering and molecular absorption, i.e., it does not correct the contribution of aerosols to the TOA reflectance. This AC was considered because it is reasonable to test the extent to which an aerosol correction is necessary over the RdP's highly reflective waters. Table 3. Specific inputs (aerosol correction type "aer_opt" and wavelengths of the corrector bands "aer_wave_short" and "aer_wave_long") of each of the studied L2Gen atmospheric corrections (ACs). The band wavelengths used for aerosol calculations are here represented with the band tags of Table 1, since they depend on each sensor. General inputs to L2Gen that were not changed among ACs are mentioned in the text. GW94-NIR, NIR multi-scattering extrapolative approach; ITER-NIR, standard NASA's iterative scheme in the NIR; GW94-SWIR, SWIR multi-scattering extrapolative approach; RC, Rayleigh correction.

Atmospheric Correction (AC) Aer_Opt (Description) Aer_Wave_Short Aer_Wave_Long
All the aforementioned ACs were run without applying the bidirectional reflectance distribution function (BRDF) correction (brdf_opt = 0), which accounts for the directionality of the upwelling radiance. This correction was not applied given that it generally degrades performance in turbid waters [25]. The gaseous absorption correction was set to account for ozone, carbon dioxide, nitrogen dioxide, and water vapor (gas_opt = 15). The land and clouds were masked (using a cloud threshold of 0.018 at the SWIR3 band) as well as pixels with high sun zenith angle (θ s > 60 • ) and view zenith angle (θ v > 70 • ). The sun glint mask was not applied (to eventually test the robustness of the PCA-SWIR schemes in sun glint conditions). Additionally, the HILT and STRAYLIGHT standard masks were switched off since they usually mask turbid water pixels erroneously. The spatial resolution was fixed at 1000 m.
For every VIIRS and MODIS scene that corresponded to a quasi-simultaneous field measurement (a maximum time difference of 30 min was considered), a window of 3 × 3 pixels was extracted, with its central pixel the closest to the reported coordinates of the field data. In the case of the coastal site at Buenos Aires located at the end of a 500 m long pier (black dot in Figure 3), the window was intentionally shifted away from land (1 pixel North, 1 pixel East in sensor geometry) to reduce the impact of land on the window. as pixels with high sun zenith angle (θs > 60°) and view zenith angle (θv > 70°). The sun glint mask was not applied (to eventually test the robustness of the PCA-SWIR schemes in sun glint conditions). Additionally, the HILT and STRAYLIGHT standard masks were switched off since they usually mask turbid water pixels erroneously. The spatial resolution was fixed at 1000 m.
For every VIIRS and MODIS scene that corresponded to a quasi-simultaneous field measurement (a maximum time difference of 30 min was considered), a window of 3 × 3 pixels was extracted, with its central pixel the closest to the reported coordinates of the field data. In the case of the coastal site at Buenos Aires located at the end of a 500 m long pier (black dot in Figure 3), the window was intentionally shifted away from land (1 pixel North, 1 pixel East in sensor geometry) to reduce the impact of land on the window. Once the window was selected, a filtering protocol was performed on a band-byband basis: 1.
If the number of AC-failure pixels (i.e., with NaN) was more than 4 out of 9 (>44.4%), the window was discarded (no match-up). Otherwise, the median and the standard deviation of the remaining pixels were taken as the reported value (x) and the absolute error (u x ) of the window, respectively; 2.
If the coefficient of variation (CV = u x /x) exceeded 20% at a given band, the station was discarded (no match-up).

Evaluation Metrics
To analyze the performance of the ACs, linear regressions between predicted (remotely estimated or simulated) and observed (in situ) reflectance were performed using a Theil-Sen linear regressor [26,27], which is less affected by outliers than least-squares regressors. From these regressions the usual linear regression statistics were obtained, i.e., slope, offset, and coefficient of determination R 2 . For the match-up exercise, the root mean square error (RMSE), the mean absolute difference (MAD), the mean difference (MD) and the mean absolute percentage difference (MAPD) were also computed as follows: where the subscript i accounts for each individual match-up, and N m is the total number of match-ups. Given that the selection criteria mentioned in Section 2.3 might yield different N m values for different ACs, the absence of effective retrieval in the computation of the absolute statistics (i.e., RMSE, MAD, and MAPD) was penalized in the following way: when the i-th station for a given AC was missing, its value was replaced by the worst value retrieved by the remaining ACs in the absolute sense (i.e., the one with highest |ρ w,i sat -ρ w,i field |). This was done to minimize biases caused by different N m among ACs. The total number of real match-ups, N m , will also be shown as a general indicator of performance of the ACs.

Study Area
The Río de la Plata (RdP) is a funnel-shaped estuary located at Eastern-Central South America ( Figure 3) and drains the second largest basin in South America, transporting between 80 and 160 million tons of sediments every year. Several measurements of suspended particulate matter (SPM) concentration were taken during field campaigns over the period 2012-2020, with values varying from 4.38 g/m 3 to 940.00 g/m 3 (median: 61.40 g/m 3 , interquartile range: 48.56 g/m 3 ). The RdP has a stationary maximum turbidity front associated with a saline front where the fresh waters of the estuary meet the South Atlantic salty waters [28,29]. The position of this front is mainly controlled by the bathymetry and coincides with the "La Barra del Indio" shoal in the north sector and with the 5 m isobath in the Samborombón Bay to the south. The RdP transports large amounts of sediments whose behavior is relevant for a number of applications in coastal areas. Several studies that make use of ocean color imagery have been conducted to determine the area of influence of the Rio de la Plata turbid plume over coastal waters of the Southwest Atlantic adjacent to the estuary [30][31][32][33]. The high-SPM and highly reflective waters of RdP and its extension (42 km to 220 km wide and 275 km long) make the estuary an excellent site to test ACs over moderate resolution ocean color radiometers such as MODIS and VIIRS.

Radiometric Measurements
A total of 55 in situ above-surface radiometric measurements collected during several field campaigns between 2012 and 2020 in the RdP (Figure 2a) were used both as input to radiative transfer simulations to theoretically test the PCA-SWIR AC and to perform match-ups with satellite data. From these 55 measurements, made over waters with SPM concentration ranging between 6 g/m 3 and 522 g/m 3 and nephelometric turbidity between 9 FNU and 708 FNU (measured using a portable HACH 2100P ISO), 28 were acquired using a set of three RAMSES/TriOS hyperspectral spectroradiometers (spectral range 350 nm to 950 nm, spectral resolution 3.3 nm, spectral discretization 2.5 nm), two measuring upwelling and downwelling radiances just above surface (L 0+ u and L 0+ d , respectively) at a relative azimuth angle from the sun of either 135 • or 90 • and reciprocal zenithal angles of θ = ±40 • , respectively, and one pointing toward zenith to measure downwelling irradiance, E 0+ d . The remaining 27 measurements were acquired with an Analytic Spectral Device (ASD FieldSpec FR) spectroradiometer (spectral range 350 nm to 2500 nm, spectral resolution 3 nm between 350 nm and 900 nm and of 10 nm to 12 nm in the range 900 nm to 2500 nm, with a spectral discretization of 1 nm) that consisted of a single radiometer connected to an optical fiber of 8 • aperture that was sequentially pointed to water and sky to determine L 0+ d and L 0+ u using the same sun-viewing geometry as with TriOS and pointed to a quasi-Lambertian Spectralon plaque to determine E 0+ d by multiplying the nadir radiance leaving the plaque by π. The reflectance of the plaque (>0.99 in the range 400-1500 nm) was assumed equal to 1. These measurements were combined to calculate the water reflectance ρ w according to where ρ M (W) is the surface radiance reflectance factor taken from Mobley 1999 [34], which depends on the wind speed, W. Since many of these measurements were taken at coastal sites or close to the coast, the relation between the surface roughness and the wind speed used to parametrize the surface radiance reflectance in Mobley 1999 no longer holds. In those cases, we fixed ρ M at 0.0256, i.e., at null wind speed. The protocols used with both radiometers followed the generic "Above-water The results presented in this and in several of the following sections vary slightly among the considered sensors, so they will not be shown for all of them. A set of PCA eigenvectors is displayed in Figure 4, calculated from SOS simulations (Set 0) in different bands/sensors for the model that used VIIRS 1241 nm, 1602 nm, and 2257 nm SWIR bands (PCA-SWIR123). Note that the eigenvectors have dimension N + 1, with N equal to the number of correction bands used in the corresponding model.
Similar to what was shown by Gross et al., 2007 [10], the first PCA eigenvector has a smooth spectral shape, showing that the magnitude-spectrally close to white-of the signal is the first cause of variability between different spectra. This eigenvector is mainly representing the different concentrations of aerosols and the residual component of the water surface Fresnel reflection, strongly dependent on the geometry of observationillumination and the surface wind. The second eigenvector has a blue signature, which is mainly due to multiple Rayleigh-aerosol scattering, ρ ra , highly variable depending on aerosol type, and more effective in the blue due to a larger Rayleigh optical thickness in that range. The third and fourth eigenvectors are associated with very small (<5% in the blue and <2% in the NIR) fractions of explained variance and represent high-order departures from the spectral shape of the main principal components. Even though they are not easily interpretable in basic physical terms, their marked spectral features suggest they are associated to the only physical component of our system whose spectral shape varies, i.e., aerosols (either through differential granulometry or chemical composition). As expected, the first two eigenvectors are associated with more than 95% of the variance, as is shown in Table 4. As an example, Table 5 shows the components of the eigenvectors of the PCA-SWIR13 AC scheme for sensors MODIS/Aqua and VIIRS/Suomi-NPP.  Similar to what was shown by Gross et al., 2007 [10], the first PCA eigenvector has a smooth spectral shape, showing that the magnitude-spectrally close to white-of the signal is the first cause of variability between different spectra. This eigenvector is mainly representing the different concentrations of aerosols and the residual component of the water surface Fresnel reflection, strongly dependent on the geometry of observation-illumination and the surface wind. The second eigenvector has a blue signature, which is mainly due to multiple Rayleigh-aerosol scattering, ρra, highly variable depending on aerosol type, and more effective in the blue due to a larger Rayleigh optical thickness in that range. The third and fourth eigenvectors are associated with very small (<5% in the blue and <2% in the NIR) fractions of explained variance and represent high-order departures from the spectral shape of the main principal components. Even though they are not easily interpretable in basic physical terms, their marked spectral features suggest they are associated to the only physical component of our system whose spectral shape varies, i.e., aerosols (either through differential granulometry or chemical composition). As expected, the first two eigenvectors are associated with more than 95% of the variance, as is shown in Table 4. As an example, Table 5 shows the components of the eigenvectors of the PCA-SWIR13 AC scheme for sensors MODIS/Aqua and VIIRS/Suomi-NPP.  The error amplification factors, i.e., the conditional numbers (Equation (6)) of the inversion matrices (left-hand side of Equation (5)) of each scheme, are reported for each band in Table 4. Large numbers are associated with matrices closer to non-invertibility and also with weaker error bounds after inversion (i.e., potentially larger error amplification). This means that the PCA-SWIR123 scheme, given the high conditional numbers in the blue (389.3 vs. 4.8), is likely to produce worse retrievals when compared to the two-SWIR-band schemes, even if the total variance explained by the total used N components is higher in the PCA-SWIR123 case (99.43 vs. 95.42). The aerosol reflectance retrieved by the PCA-SWIR12 AC applied to the SABIA-Mar BLUE and NIR1 bands in the conditions used to generate the PCA eigenvectors (Set 0), ρ PCA a , is displayed in Figure 5. This exercise allows a first evaluation of the aerosol reflectance estimation prior to the diffuse transmittance correction. In concomitance with what has been discussed in the previous sections, the performance is markedly worse in the BLUE than in the NIR due to a higher spectral distance to the SWIR spectral region. It must be noticed that the PCA-SWIR12 scheme was implemented-with no intermediate steps involving the NIR bands-to estimate ρ PCA a (412 nm). In the case of VIIRS bands (Figures 6 and 7), a better performance is observed with PCA-SWIR123 in the NIR when compared to two-band models (slope closer to 1 and lower MAD). On the contrary, we observe worse retrievals in the blue, i.e., for ρ a (489 nm), using the three-band model. This is also evident from an abrupt increase in the conditional number of the inversion matrix, as the band to be corrected is spectrally further away from the SWIR bands. The larger conditional number in the blue than in the NIR occurs for all models (both with N = 2 and N = 3, although markedly larger in the latter case (see Table 4), indicating that the inversion matrix of the system of Equation (5) is closer to non-invertibility, which in turn is associated with a low correlation between the bands in the blue and those in the SWIR. It may be seen as counter-intuitive that adding information (i.e., adding a band) can degrade the performance of the algorithm, but it must be emphasized that this only occurs in the blue bands where the correlation to the SWIR is lower than at higher wavelengths. In the case of the model with three corrective bands, this lack of correlation and the addition of a component leads to an overrepresentation of the variability of the reflectance in the blue, implying larger propagation errors. BLUE and NIR1 bands in the conditions used to generate the PCA eigenvectors (Set 0), ρ a PCA , is displayed in Figure 5. This exercise allows a first evaluation of the aerosol reflectance estimation prior to the diffuse transmittance correction. In concomitance with what has been discussed in the previous sections, the performance is markedly worse in the BLUE than in the NIR due to a higher spectral distance to the SWIR spectral region. It must be noticed that the PCA-SWIR12 scheme was implemented-with no intermediate steps involving the NIR bands-to estimate ρ a PCA (412 nm). In the case of VIIRS bands (Figures 6 and 7), a better performance is observed with PCA-SWIR123 in the NIR when compared to two-band models (slope closer to 1 and lower MAD). On the contrary, we observe worse retrievals in the blue, i.e., for ρa(489 nm), using the three-band model. This is also evident from an abrupt increase in the conditional number of the inversion matrix, as the band to be corrected is spectrally further away from the SWIR bands. The larger conditional number in the blue than in the NIR occurs for all models (both with N = 2 and N = 3, although markedly larger in the latter case (see Table  4), indicating that the inversion matrix of the system of Equation (5) is closer to non-invertibility, which in turn is associated with a low correlation between the bands in the blue and those in the SWIR. It may be seen as counter-intuitive that adding information (i.e., adding a band) can degrade the performance of the algorithm, but it must be emphasized that this only occurs in the blue bands where the correlation to the SWIR is lower than at higher wavelengths. In the case of the model with three corrective bands, this lack of correlation and the addition of a component leads to an overrepresentation of the variability of the reflectance in the blue, implying larger propagation errors.

NIR Water Reflectance Retrieval (SOS Set W)
The retrieved (PCA-SWIR) vs. observed (ASD) water reflectance values in the SABIA-Mar and VIIRS NIR bands are shown in Figures 8 and 9 using all the PCA-SWIR possible combinations of SWIR bands, i.e., one in the case of SABIA-Mar (lacking the SWIR3 band) and four in the case of VIIRS. The results obtained for each of the 22 field reflectance measurements used in the computation of Set W are displayed as the mean values of retrieved ρw PCA (red dots) and the interquartile ranges (IQRs) are represented with red vertical bars. In general, a good correspondence is observed between the estimated and measured reflectance values, with low biases and slopes and correlation coefficients close to 1. Similarly to what was found for the estimates of the aerosol reflectance using the simulations of Set 0, Figure 9 shows that the three-band model presents lower IQRs, i.e., less sensitivity to atmospheric variability, and slopes closer to 1 compared to the two-band schemes.

NIR Water Reflectance Retrieval (SOS Set W)
The retrieved (PCA-SWIR) vs. observed (ASD) water reflectance values in the SABIA-Mar and VIIRS NIR bands are shown in Figures 8 and 9 using all the PCA-SWIR possible combinations of SWIR bands, i.e., one in the case of SABIA-Mar (lacking the SWIR3 band) and four in the case of VIIRS. The results obtained for each of the 22 field reflectance measurements used in the computation of Set W are displayed as the mean values of retrieved ρ w PCA (red dots) and the interquartile ranges (IQRs) are represented with red vertical bars. In general, a good correspondence is observed between the estimated and measured reflectance values, with low biases and slopes and correlation coefficients close to 1. Similarly to what was found for the estimates of the aerosol reflectance using the simulations of Set 0, Figure 9 shows that the three-band model presents lower IQRs, i.e., less sensitivity to atmospheric variability, and slopes closer to 1 compared to the twoband schemes.   In order to study the impact of noise on the performance of the PCA-SWIR scheme, we estimated a noise-equivalent RC reflectance (NEρRC) from MODIS/Aqua imagery. We decided not to use the classic homogeneous area approach of Duggin et al., 1985 [41] given that it supposes a homogeneous "real" field (the field at zero noise) in the selected area, a condition that is never met in the Río de la Plata at the VIS/NIR bands. Instead, we imple-

Effect of Noise on MODIS/Aqua Imagery
In order to study the impact of noise on the performance of the PCA-SWIR scheme, we estimated a noise-equivalent RC reflectance (NEρ RC ) from MODIS/Aqua imagery. We decided not to use the classic homogeneous area approach of Duggin et al., 1985 [41] given that it supposes a homogeneous "real" field (the field at zero noise) in the selected area, a condition that is never met in the Río de la Plata at the VIS/NIR bands. Instead, we implemented a geostatistical method, which was first applied to satellite (AVIRIS) images by Curran and Dungan 1989 [42]. This method assumes a Gaussian noise model, a null autocorrelation of the noise field, and a null correlation between the noise and the "real" field. After having verified the geostatistical method on synthetic images of known noise, it was applied to a set of RC reflectance rasters (L2 product ρ s using SeaDAS v7.5) of 20 MODIS/Aqua images of the Río de la Plata over the period 2013-2019 over a small subregion between 35.26 • S and 35.52 • S and between 57.05 • W and 56.70 • W (windows of 40 × 40 pixels approximately, see Figure 3).
Taking the mean of the noise amplitudes in each of the selected subscenes described in the Methods section, the absolute value of the noise, σ, was obtained (Table 6) by estimating the nugget of the associated semi-variogram of each subscene, following [42]. These values were used to theoretically evaluate the impact of noise on the AC performance on MODIS/Aqua. It was assumed that the noise in the individual bands was not correlated spectrally. Table 6. Noise-equivalent RC reflectance for each of the MODIS/Aqua sensor bands in the NIR and SWIR, obtained by applying the geostatistical method to the RC reflectance of the mentioned set. The effect of noise in the MODIS/Aqua bands on the performance metrics of the PCA-SWIR algorithm over the set of radiative transfer simulations was found to be low, if taking into account the low signal-to-noise ratios of the SWIR bands (Table 7). Table 7. Comparison between performance metrics of simulated water reflectance ρ w (PCA) vs. ρ w (field) at MODIS/Aqua 857 nm band (Set W) with and without noise added to simulations. MAD, mean absolute difference; Int, intercept. It should be noted that, as already known, a normally distributed noise-of mean 0-is propagated through a linear system into a normally distributed noise-of mean 0-of the retrieved quantity, with the amplification of the noise regulated by the conditional number of the inversion matrix of the scheme. In our case, the output quantity is the aerosol-and eventually, the water-reflectance in the NIR bands. We consider this method to be more robust statistically under the effect of noise in comparison to the standard procedure of using ratios such as ε(λ 1 , λ 2 ), given that ratios of variables with normally distributed errors are prone to magnifying the effect of noise over the outputs (especially under low mean values of the denominator, i.e., under clear atmospheric conditions). This will be also discussed in the following section.

Match-ups
The global scatterplots of the match-up exercise are shown for the NIR1~750 nm ( Figure 10) and NIR2~860 nm ( Figure 11) bands of all the sensors, and the global performance metrics are shown in Table 8. We decided not to differentiate TriOS and ASD field reflectance measurements graphically since no systematic bias was observed between them. Inside each inset, the satellite-derived water reflectance (ρ w sat ) was compared with the fieldmeasured water reflectance (ρ field w ) for ITER-NIR, RC, GW94-SWIR13 and PCA-SWIR13 in the NIR bands) (see Table 3 for a description of these acronyms). It should be noticed that different total numbers of effective match-ups are obtained for each AC given the different AC failure and acceptance ratios after following the quality control criteria described in Section 2.3. The NIR-iterative (ITER-NIR, a) and standard NIR-black-water (GW94-NIR, b) ACs tend to underestimate water reflectance at all wavelengths and especially at high field-reflectance values, frequently failing to retrieve numerical values (seen as lower N values) or otherwise retrieving negative reflectance. reflectance measurements graphically since no systematic bias was observed between them. Inside each inset, the satellite-derived water reflectance (ρw sat ) was compared with the field-measured water reflectance (ρ w field ) for ITER-NIR, RC, GW94-SWIR13 and PCA-SWIR13 in the NIR bands) (see Table 3 for a description of these acronyms). It should be noticed that different total numbers of effective match-ups are obtained for each AC given the different AC failure and acceptance ratios after following the quality control criteria described in Section 2.3. The NIR-iterative (ITER-NIR, a) and standard NIR-black-water (GW94-NIR, b) ACs tend to underestimate water reflectance at all wavelengths and especially at high field-reflectance values, frequently failing to retrieve numerical values (seen as lower N values) or otherwise retrieving negative reflectance.    On the other hand, even if the results of the PCA-SWIR13 scheme in the blue and green show general overestimation (MD = 0.019 and MD = 0.011, respectively, see Table 8), they still show better agreement with field data in comparison with those of the GW94-SWIR13 approach (slopes closer to 1, less negative retrievals, and higher N-i.e., fewer AC failures).
In the NIR bands, the Rayleigh correction approach (RC, Figures 10b and 11b) tends to overestimate water reflectance (with intercepts of 0.012 and 0.014 and MD of 0.0103 and 0.0103 in bands NIR1 and NIR2, respectively) in comparison with the other schemes, all of which intend to remove the aerosol effect. The PCA-SWIR13 AC clearly shows an improved performance over the standard NIR procedures. When compared to GW94-SWIR13, a slightly better performance is attained in the NIR1 band-considering the obtained RMSE, MAD, MD, MAPD, slope, and R 2 , and a clear improvement is observed in the NIR2 band, with retrieved slope and R 2 close to 1 and RMSE, MAD, MAPD, and MD closer to 0. A subset of field water reflectance spectra together with the corresponding satellitederived reflectance spectra from sensors MODIS/Terra (MT) and VIIRS/Suomi-NPP (VS-NPP) is shown in Figure 12. These examples show general agreement in spectral magnitude and shape between field and PCA-SWIR13 spectra, with a tendency to achieve better performance in the green-red-NIR region of the spectrum when compared to ITER-NIR, GW-NIR, and GW-SWIR13 ACs. In the blue region, even though the PCA-SWIR13 scheme comparatively performs better than the other schemes, with fewer AC failures and negative retrievals, it tends to overestimate reflectance in the blue. This overestimation may be associated with a magnification of the error produced by higher condition numbers of inversion (Equation (6)) and the simple model of the diffuse transmittance factor expressed in Equation (7) (see Discussion section).

Spatial Analysis of Water and Aerosol Patterns
Another way of evaluating the quality of an AC is to compare the spatial patterns of the water and aerosol reflectance. If the algorithm is capable of decoupling both signals spectrally on a pixel-by-pixel basis, it is also expected that it will do so spatially, even though water properties/processes may partly affect local aerosol production. When comparing the performance of PCA-SWIR13 with other ACs (Table 3), we found lower spatial correlations and lower percentages of AC failures and negative water reflectance pixels, as well as less noisy aerosol retrievals in the VIS/NIR with the PCA-SWIR13 AC. This is exemplified in Figures 13 and 14, which display the water and aerosol reflectance maps produced for the MODIS/Aqua scene of 2014-05-16T17:35:00Z using the PCA-SWIR13, GW94-NIR, ITER-NIR, and GW94-SWIR13 ACs in the blue (443 nm) and NIR2 (859 nm) bands, respectively. In both bands, the PCA-SWIR13 AC yields the least cross-contamination of the water and aerosol reflectance, fewer AC failures (shown in white), negative reflectance pixels, and lower slopes on the ρ a vs. ρ w relation. Even though a lower slope is observed for GW94-SWIR13 than for PCA-SWIR13 at 443 nm, this may be associated with the much shorter range of retrieved water reflectance. In general, the performance is markedly improved when compared to standard NIR approaches (GW94-NIR and ITER-NIR). Note that the GW94-NIR AC is not shown in Figure 14 because this scheme retrieves trivially zero water reflectance in the NIR bands. When comparing PCA-SWIR13 with GW94-SWIR13, the slightly lower slope on ρ a vs. ρ w indicates that in the selected region (shown inside squared contours), where the maximum turbidity front of Río de la Plata is located [43], there is a slightly better water reflectance estimation. Additionally, a markedly lower percentage of negative reflectance pixels (0.00% vs. 7.11% in the BLUE band; 1.26% vs. 5.72% in the NIR2 band), lower percentage of AC failures (0.45% vs. 2.31% in the blue band; 0.45% vs. 0.71% in the NIR2 band) is observed. It is worth noting that PCA-SWIR13 aerosol reflectance maps are less noisy compared to the traditional extrapolative approach (compare insets (c) and (f) in both figures).    Comparison of spatial patterns of water and aerosol reflectance retrieved by PCA-SWIR13, GW94-SWIR13, GW94-NIR, and ITER-NIR ACs from MODIS/Aqua image of the Río de la Plata acquired on 2014-05-16T15:35:00Z. RGB composite (a), water (b,e,h,k), and aerosol (c,f,i,l) reflectance at 443 nm, retrieved using PCA-SWIR13 (b,c), GW94-SWIR13 (e,f), GW94-NIR (h,i), and ITER-NIR (k,l) schemes. Plots show aerosol vs. water reflectance in the Río de la Plata's turbidity front region marked with squares (d,g,j).

Discussion
The new AC presented here has multiple advantages. First of all, it relies on simple and plausible physical hypotheses which, given the simplicity of the algorithm, can be tracked and tested with ease. It also has the advantage that the generalization of the inversion procedure to the use of multiple sets of SWIR correction bands is straightforward, which is not the case of standard extrapolative procedures that rely on ratios over pairs of bands. Additionally, given the linear nature of the inversion scheme, this AC keeps the structure of the noise and also tends to bound the noise better than in the case of band ratios. This is a key point, since ratios of normally distributed random variables (i.e., assuming a Gaussian noise model), such as ε(λ 1 , λ 2 ), might produce drastic amplifications of noise when the mean value of the denominator (i.e., the aerosol signal at λ 2 ) approaches 0, i.e., in clear atmospheric conditions. This noise amplification combined with the error in the aerosol model produces frequent negative reflectances in spectrally distant bands, such as the blue. The PCA-SWIR AC has been able to show a considerable reduction of this undesired result. Although PCA-SWIR has shown general overestimation of the water reflectance in the blue, this might be easily circumvented in follow-up studies by enlarging or improving the radiative transfer simulation set (such as improving the quality of aerosol models or, e.g., assuming a spectrally dependent air-water relative refractive index, [44]) and revising the transmittance correction expression, whose error is amplified in the blue, given the Angstrom law specified for the aerosol optical thickness. The basic limitation in the blue, however, is that blue wavelengths are far from SWIR wavelengths, making the PCA representation of the atmospheric signal inaccurate. Reanalysis data (e.g., MERRA-2) could also be used to reduce uncertainties on aerosol transmittance.
The PCA-SWIR algorithm presents the main restriction that is only applicable to sensors containing at least two operational SWIR bands where the black water assumption is valid. This is an important constraint since many optical sensors lack these bands given their high marginal costs, and these bands tend to have lower associated signal-to-noise ratios when compared to the VIS and NIR bands. On the other hand, even though the performance of the AC was validated using several sources, it should be noted that the Río de la Plata field reflectance data used in the study are not evenly distributed among the complete range of expected values in the region; they are mostly representative of regimes with moderate suspended matter concentration. Therefore, the inclusion of a higher number of extremely high reflectance field measurements to the existing in situ database is highly desirable to assess the comparative performance of these ACs at extremely turbid regimes. Apart from these restrictions, the PCA-SWIR scheme(s) still has to be evaluated in optically complex waters under the effect of moderate to high sun glint. In follow-up studies it will be also necessary to perform intercomparisons between PCA-SWIR routines that make use of different SWIR band combinations and compare results over each specific band of each specific sensor. It will also be interesting to compare the PCA-SWIR scheme to a hybrid PCA and standard extrapolative scheme, in which water signal in the NIR is first estimated using the PCA-SWIR approach and water signal in the VIS is then computed using the standard approach (Gordon and Wang 1994) after subtracting the water signal in the NIR bands.

Summary and Conclusions
In this study an atmospheric correction algorithm (herein PCA-SWIR) is proposed to estimate the reflectance of Río de la Plata turbid waters in the NIR bands of moderateresolution sensors such as MODIS, VIIRS, and SABIA-Mar, which possess high-quality far SWIR bands. The procedure is based on the decomposition into principal components of the atmospheric VIS-NIR-SWIR signal, which is a different approach from the traditional "epsilon" approach to extrapolate from the SWIR bands to the NIR bands. The weight of the eigenvectors of the variance-covariance matrix of the set of simulations is determined from the signal at the SWIR bands for each pixel. These weights are then used to calculate the aerosol signal in the VIS-NIR. A simple diffuse transmittance model is applied to correct the water signal and finally yield the water reflectance.
The algorithm was theoretically tested using a set of simulated Rayleigh-corrected reflectances for a wide range of expected atmospheric conditions and field reflectances of Río de la Plata turbid waters. A geostatistical method, previously tested on synthetic images, was used to estimate the noise-equivalent Rayleigh-corrected reflectance of MODIS/Aqua images to theoretically assess the effect of noise on the water reflectance retrieval. These results show that even though SWIR bands tend to have a low signal-to-noise ratio in comparison to NIR and VIS bands, the high levels of noise have a relatively low impact on the retrievals in comparison to traditional extrapolative ACs based on ratios. When considering the effect of noise in the bands involved in the scheme, the relationships between estimated and observed water reflectance tend to show slopes slightly farther from 1, slightly higher biases, and slightly lower correlation coefficients. Nonetheless, the overall observed degradation effect on the performance is acceptably low. This is also observed in the retrieved aerosol reflectance maps over Río de la Plata, which tend to show smoother noise patterns in comparison with the traditional extrapolative SWIR approach.
Although this algorithm was originally conceived to estimate the water signal in the NIR in order to subsequently subtract it from the total signal and proceed with the usual extrapolation schemes, its design is not strictly limited to the NIR, so its performance in visible bands was also tested. As expected, given the decreasing correlations between the aerosol signal of the SWIR bands and the bands of interest, as the spectral distance increases, the water reflectance retrievals are not as accurate as in the NIR but are still better in comparison to existing standard AC approaches.
The theoretical results obtained in the NIR bands are optimal, with low biases and slopes and correlation coefficients close to 1. In particular, the PCA-SWIR AC that uses three SWIR bands theoretically performs better in the NIR than those that use two bands, given that the former manages to represent the atmospheric signal better than the latter. The PCA-SWIR AC is therefore useful for estimating turbidity and suspended particulate matter concentration from water reflectance in NIR spectral bands. The reverse occurs in the blue bands, where the combination of a low correlation between the atmospheric signals in the blue and SWIR and the use of three main components leads to an overrepresentation of the variability of the signal in the blue (a phenomenon that is quantitatively determined with high condition numbers of the inversion matrixes).
The match-up exercise showed frequent failures and negative reflectance retrieval using the standard approaches and a general overestimation using the just-Rayleighcorrection approach in the NIR. The SWIR approaches that use the SWIR1 (~1240 nm) and SWIR3 (~2130 nm) bands show better general performance in comparison to other choices of SWIR bands. In particular, the PCA-SWIR13 performs markedly better in the NIR2 (~860 nm) band compared to the traditional SWIR approach, and slightly better in the case of the NIR1 (~750 nm) band.
The proposed AC can be applied to any sensor containing at least two bands in the far-SWIR (i.e., where the black water assumption holds), including Landsat-8/OLI and Sentinel-2/MSI imagers and is also applicable to optically complex waters other than Río de la Plata.

Data Availability Statement:
The data presented in this study are not public but are potentially available on request from the corresponding author.