Estimation of Water Quality Parameters in Oligotrophic Coastal Waters Using Uncrewed-Aerial-Vehicle-Obtained Hyperspectral Data

Abstract

Water quality monitoring in coastal areas and estuaries poses significant challenges due to the intricate interplay of hydrodynamic, chemical, and biological processes, regardless of the chosen monitoring methods. In this study, we analyzed the applicability of different monitoring sources using in situ data, uncrewed-aerial-vehicle (UAV)-mounted hyperspectral sensing, and Sentinel-2-based multispectral imagery. In the first part of the study, we evaluated the applicability of existing empirical algorithms for water quality (WQ) parameter retrieval using hyperspectral, simulated multispectral, and satellite multispectral datasets and in situ measurements. In particular, we focused on three optically active WQ parameters: chlorophyll a ($Chl,a$), turbidity (TUR), and colored dissolved organic matter (CDOM) in oligotrophic coastal waters. We observed that most existing algorithms performed poorly when applied to different reflectance datasets, similar to previous findings in small and optically complex water bodies. Hence, we proposed a novel set of locally based empirical algorithms tailored for determining water quality parameters, which constituted the second part of our study. The newly developed regression-based algorithms utilized all possible combinations of spectral bands derived from UAV-generated hyperspectral data and exhibited coefficients of determination exceeding 0.9 for the three considered WQ parameters. The presented two-part approach was demonstrated in the semi-enclosed area of Kaštela Bay and the Jadro River estuary in the Central Eastern Adriatic Sea. This study introduces a promising and efficient screening method for UAV-based water quality monitoring in coastal areas worldwide. Such an approach may support decision-making processes related to coastal management and ultimately contribute to the conservation of coastal water ecosystems.

1. Introduction

Coastal environments, as transitional ecosystems between land and sea, are complex areas that ensure sufficient resources for different living organisms. Moreover, their relevance is reflected in enabling social and economic benefits to humans [1]. However, the continuous expansion of urban settlements increases the pressure on these vulnerable systems [2], primarily by inputting anthropogenic pollutants into coastal waters. Various contaminants (pharmaceuticals, consumer and industrial products, microbiological pollutants) enter coastal waters through sewage and runoff from agricultural lands. Those discharges potentially cause excessive nitrogen and phosphorus concentrations. Such nutrient intrusion could eventually lead to eutrophication, which depends on several factors (water residence time, tidal exchange, stratification) as reported by Pinckney et al. [3]. Consequently, that can disturb the ecosystem’s natural balance and result in the growth of harmful algae, increased turbidity, and oxygen depletion [4,5]. In addition to algal blooms, the discharge of untreated domestic wastewater introduces fecal matter into rivers and estuaries, which promotes the growth of various bacterial and viral pathogens [6]. Thus, the regular monitoring of coastal areas is crucial to ensure the marine ecosystem’s quality and sustainability.

Historically, water body monitoring programs have heavily relied on labor-intensive and time-consuming in situ measurements [7]. Nevertheless, over the past two decades, the integration of remote sensing data has enhanced these monitoring systems by overcoming limitations in spatial and temporal resolution [8]. While traditional in situ methods provide detailed, localized information about water quality parameters, the advent of UAV-based remote sensing has introduced a complementary approach. This approach allows for broader spatial coverage and captures dynamic variations in water quality (WQ) over larger areas [9]. A review by [10], which covers various implementations of remote sensing in WQ assessment, identifies airborne sensors as particularly suitable for monitoring smaller water bodies, such as rivers and estuaries.

The most commonly developed remote sensing (RS) algorithms for WQ assessment include different band combinations of the visible (VIS) and infrared (IR) electromagnetic spectrum, depending on the area and WQ parameters being studied. Historically, algorithms were developed for ocean application, and they needed to be adjusted when used in coastal areas due to the different optical properties. Thus, the algorithms developed for the coastal waters are usually site-specific and regionally dependent.

Some studies used hyperspectral (HS) sensors, such as portable spectrometers (e.g., [11,12,13]) or boat-mounted sensors [14], and recently, even sensors fixed at the point of interest above the river surface [15]. The usage of uncrewed aerial vehicles (UAVs) with HS sensors, providing hundreds of narrow contiguous bands, has emerged as an important complement for monitoring micro-areas like bays and estuaries. The UAV-mounted HS sensors provide data with higher spatial, temporal, and spectral resolution than data obtained from satellites or point spectrometers [16] and can be used with low cost, high safety, and flexibility in flight planning [17]. With the ability to timely detect different small-scale features, the UAV technology compensates for most shortcomings of traditional RS technology, and their usage in WQ assessment has been growing at an accelerating pace [9,18,19,20,21,22,23].

Regardless of the type of water body, phytoplankton is one of the most important parameters for water quality (WQ) assessment, encompassing attributes such as bloom frequency, composition, and abundance. However, the practical challenges of assessing all these attributes have led to the common use of chlorophyll a ($Chl,a$) as a reliable proxy. For instance, a recent study conducted by [24] explored the feasibility of using a five-band multispectral sensor mounted on a UAV to estimate $Chl,a$ concentrations. Additionally, turbidity (TUR) and colored dissolved organic matter (CDOM) are employed to describe other WQ properties, such as the distribution of total suspended solids (indicative of water clarity) [25], and to detect the presence of surface runoff and point source pollution [26,27]. They also help assess water color [28]. Given their optical properties, $Chl,a$, TUR, and CDOM are commonly estimated using remote sensing technology, and these parameters will be utilized in the present study.

In this study, we investigate how to combine RS tools with UAV HS sensors and in situ data to develop new algorithms and produce maps for selected WQ parameters. The combination of these monitoring data sources has been approached from different aspects within the last two decades [29,30,31,32,33], and although there has been technological advancement, the solution still seems elusive. Hence, the problem formulation in this study builds on those foundations by further investigating how one can utilize UAV HS sensors, in situ monitoring, and satellite multispectral (SAT MS) imagery to obtain maps of WQ parameters. A similar setup, but with a boat-mounted HS point sensor instead of a UAV, was implemented on the river Elbe in Germany [14]. In this case, we focused on a small coastal water body (estuary), characterized by complex optical properties and rather a dynamic ecosystem behavior, which sets high requirements for monitoring procedures under different EU regulations [34,35,36].

The main goal of this study was to assess the effectiveness of existing algorithms and subsequently develop customized empirical algorithms that cater to the distinctive oligotrophic characteristics of our case study. Our particular focus was on utilizing the capabilities of UAV hyperspectral bands to accurately derive water quality parameters. To achieve this, we analyzed information from three distinct spectral data sources: hyperspectral data, simulated multispectral data, and satellite multispectral data. Through this evaluation process, we determined that hyperspectral data exhibited superior performance compared to other sources, thereby facilitating the development of site-specific regression-based algorithms for our study area.

2. Materials and Methods

The problem formulation in this study consists of a multiscale monitoring setup, followed by relevant data processing and algorithm analysis. Due to the different scales of data acquisition, we have obtained (1) a dataset from in situ measurements of $Chla$, TUR, and CDOM, (2) the UAV HS dataset, and (3) SAT MS data from Sentinel-2. For comparison purposes, we generated an additional (4) multispectral dataset from UAV HS data named UAV MS, explained in Section 2.3. Within algorithm analysis, we tested the most relevant algorithms from the literature using the above-mentioned datasets. Finally, the high-resolution HS data in this case study enabled the search for more appropriate site-specific algorithms with acceptable accuracy. The implemented methodology is illustrated in the flowchart (Figure 1) and will be further described in the following sections.

2.1. Study Area

Kaštela Bay is a semi-enclosed basin on the eastern Adriatic coast near Split, the second largest city in Croatia (Figure 2). Most Croatian coastal areas, including Kaštela Bay, rely heavily on beach-based tourism for economic growth and development [37]. The climate in this area is of Mediterranean type with hot and dry summers [38]. The average annual air temperature for Kaštela is 13.7 °C, while the long-term average annual precipitation is 1029 mm [39]. The average evaporation rate in the Adriatic is $2.8·{10}^{-8}\mathrm{m}/\mathrm{s}$, as reported by [40]. Hot summers contribute to higher evaporation rates in this region compared to the northern Adriatic.

With a maximum depth of 56 m and an average depth of 23 m, the bay is described as shallow and has a total surface of 61 km². Agricultural areas are located on the bay’s northern coast, while municipal and industrial wastewater is discharged into the eastern part of the bay. Although there have been positive results in recent years, Kaštela Bay has been well known as one of the most polluted areas in the eastern Adriatic coast due to the urban expansion lacking adequate wastewater collection and treatment [41]. Several studies were conducted in the bay, observing the negative effect of uncontrolled wastewater discharge and nutrient loadings [42,43,44].

The eastern part of the bay is dominated by the estuary of the karst river Jadro with an average annual discharge of 9.5 m³/s. The river Jadro has a 4.2 km long watercourse that flows through the small town of Solin near Split before it flows into the Kaštela Bay. The entire city of Split and the nearby region are supplied with water from the Jadro spring. As the Jadro River is the main freshwater inflow and nutrient supplier for Kaštela Bay [42], we chose this area for our testing site (Figure 2).

2.2. Data Collection

The monitoring setup applied in this study consisted of in situ measurements of WQ parameters, which were accompanied by simultaneous UAV HS operations retrieving reflectance data on a pixel-size scale (2 × 2 cm), and satellite reflectance data from Sentinel-2 overpass in the nearest time frame to the measurement campaign on a 10–20 m spatial scale. During June and November 2021 and March and April 2022, four measurement campaigns were conducted in the Jadro River estuary. All field cruises were conducted under clear weather conditions between 10:00 and 13:00 and during four seasons to cover annual temporal variations of WQ parameters. Average river flow during the week of the campaigns was 4.1 m³/s and 4.3 m³/s for June and November 2021, followed by 7.6 m³/s and 6.1 m³/s for March and April 2022. Wind speed was below the limit defined by the UAV in all four campaigns, i.e., less than 7 m/s, and no precipitation was present (details available in Supplementary Materials, Table S1).

2.2.1. In Situ Measurements

In situ measurements of WQ parameters were conducted using a boat (Nautika 600K, by Remia-plast Ltd., Zagreb, Croatia) with calibrated submersible C3 fluorometer (Turner Design), which is equipped with optical sensors for $Chla$, TUR, and CDOM. Implemented sampling interval was one second. As described by [45], the $Chla$ measurements were obtained using a blue mercury lamp, with peak emission at 460 nm and fluorescence collection at 680 nm. The turbidity sensor includes an IR lamp, with peak emission at 850 nm, and the collection of scattered light at 90°, while CDOM was measured using a UV LED (central wavelength: 365 nm), with peak emission at 350 nm and fluorescence collection at 450 nm.

2.2.2. UAV Operations

We employed a customized and robust UAV Aermatica BLY-O octocopter platform equipped with DJI Ronin-MX gimbal system to carry the hyperspectral imaging sensor. Nano-Hyperspec^® by Headwall Photonics is a push-broom sensor acquiring images line by line along the flight direction, with each line captured at a specific time corresponding to the UAV’s current position and attitude. Each line comprises 640 spatial pixels containing data acquired in 273 spectral bands. The sensor collects data in the VIS and NIR spectrum from 400 to 1000 nm, with a spectral resolution of 2.2 nm and a full-width half maximum (FWHM) of about 6 nm. Before image acquisition, the two flight polygons were set: one above the seawater and the second above the reference reflectance tarp. Due to the local flying zone restrictions, the flight altitude was 50 m during all cruises, with a total polygon length of 500 m and a width of 40 m.

2.2.3. Satellite Imagery

A high-performing Sentinel-2 satellite was chosen due to the product’s free access through the European Space Agency (ESA) data hub and a short revisit time of 5–10 days. Sentinel-2 mission consists of two satellites (Sentinel-2A and Sentinel-2B), both containing multispectral instrument (MSI) that can sample a total of 13 spectral bands (B1–B12, and B8a) with different spatial resolutions and central wavelengths [46]. In this study, we used Sentinel Level-2A atmospherically corrected surface reflectance data processed with the Sen2Cor processor [47].

2.3. Data Processing

Multiple processing stages are performed on the collected raw reflectance data to obtain the three different reflectance datasets (UAV HS, UAV MS, and SAT MS), which, together with in situ measurements, are used for algorithm analysis in Section 2.4.

2.3.1. In Situ Data

Multi-parameter probe (C3 fluorometer), as described in Section 2.2.1, recorded numerical data for $Chla$ (μg/L), TUR (NTU) and CDOM (ppb) along the water column up to the depth of assessed visibility. A total of 32 sampling locations were considered for the present analysis due to their overlap with the UAV overflight (a more detailed presentation can be found in Supplementary Material, Figure S1).

2.3.2. Hyperspectral Data

The HS imaging sensor (Nano-Hyperspec^®) records data in unitless digital numbers (DN), which are then processed using Headwall’s SpectralView application (within software package Hyperspec^® III 3.1.3) to obtain the reflectance dataset. After obtaining the reflectance values, we selected target areas of a 2 m diameter circle (defined by the on-boat GPS accuracy level) overlapping the in situ points and extracted HS data for each location using Global Mapper 22.0 [48]. All further analysis and visualizations were performed in Matlab [49]. We removed the overexposed pixels from all selected pixels and applied the Savitzky–Golay filter with a range of 15 bands and the 3rd-order polynomial to smooth the spectral reflectance [50]. Additionally, sensitivity analysis was conducted by testing the chosen existing algorithms on raw reflectances (without implementing the denoise procedure), with Savitzky–Golay filter, and with simple scatter correction (standard normal variate) [50], which can be found in Supplementary Materials (Tables S2–S4). The described preprocessing procedure was applied for each of the 32 points obtained during the monitoring period, as presented in Section 3.1.

2.3.3. Multispectral Data

MS data from Sentinel was acquired as S2-MSI granules for June and November of 2021, and for March and April 2022 from the Copernicus Open Access Hub [51] corresponding to the time frame of measurement campaigns. Then, we formed the satellite MS dataset by extracting the reflectance date from the granules matching the locations of sampling points. Since the MS instrument (MSI) mounted on Sentinel satellites has distinctively different properties than the HS sensor we used in our measurement campaigns, the SNV-processed reflectance for both satellite and UAV data are compared in Supplementary Materials (Figure S2).

SNV-processed data enable one to mitigate the impact of measuring sensors and environmental conditions relevant to optical measurements. Apparent differences in the density of spectral data are furthermore visible in Figure 3, where we depicted Sentinel bands relevant to the domain measurable by the HS camera, namely from B2 up to B8 and B8a, respectively. Thin blue bars represent the water response spectra from this study but also indicate the spectral density of the HS dataset.

Given the accessibility of satellite MS data, most existing and widely used algorithms are developed based on such data. In this analysis of algorithms, we assess the suitability of HS and MS data obtained from UAV and satellites. To facilitate a meaningful comparison, we generated a simulated MS dataset (UAV MS), as illustrated in Figure 3. This UAV MS dataset was created by spectrally resampling HS data using Sentinel-2 response functions recently updated by ESA [52]. It is worth noting that this approach of generating a simulated MS dataset has been employed in prior research, specifically for Natura 2000 areas identification, as demonstrated by [53].

2.4. Algorithm Analysis

Following the approach depicted in Figure 1, in this section, we tested several existing algorithms and searched for new potentially improved regression models.

2.4.1. Existing Algorithms

We tested the performance of 18 selected algorithms developed for coastal areas to estimate the three WQ parameters ($Chla$, $TUR$, and $CDOM$) in the Jadro River estuary in Croatia using HS data from the UAV (

Table 1. The list of existing algorithms tested for chlorophyll $a$ ($Chla$), turbidity (TUR), and colored dissolved organic matter (CDOM).

Chl a
Algorithm	Relation	Corresponding Sentinel-2 bands	Range $\left[mg/m^3\right]$	Reference
${\theta }_{Chla}^1$	${(35.75·R_{708}/R_{665}-19.3)}^{1.124}$	B5, B4	2.27–81.7	[54]
${\theta }_{Chla}^2$	$R_{443}/R_{555}$	B1, B3	0.03–2.75	[55]
${\theta }_{Chla}^3$	$R_{490}/R_{555}$	B2, B3	0.003–7.06	[56]
${\theta }_{Chla}^4$	$48.4579·{\tilde{R}}_{681}/{\tilde{R}}_{660}-48.1727$	B4, B4	0.24–7.02	[13]
${\theta }_{Chla}^5$	$-26.4373·{\tilde{R}}_{498}/{\tilde{R}}_{518}+29.0687$	B2, B2	0.24–7.02	[13]
${\theta }_{Chla}^6$	$R_{B5}-R_{B4}/R_{B5}+R_{B4}$	B4, B5	0.9–28.17	[57]
${\theta }_{Chla}^7$	$146.5·{\tilde{R}}_{687}/{\tilde{R}}_{672}-141.9)$	B4, B4	0.8–19.8	[30]
${\theta }_{Chla}^8$	$R_{B5}/R_{B4}$	B5, B4	/	[58]
${\theta }_{Chla}^9$	$\left(1/R_{B4}-1/R_{B5}\right)·R_{B8a}$	B4, B5, B8a	/	[58]
TUR
Algorithm	Relation	Corresponding Sentinel-2 bands	Range $\left[NTU/FTU\right]$	Reference
${\theta }_{TUR}^1$	$228.1·{\tilde{R}}_{645}/\left(1-{\tilde{R}}_{645}/0.1641\right)$	B4	1.8–988	[59]
${\theta }_{TUR}^2$	$3078.9·{\tilde{R}}_{859}/\left(1-{\tilde{R}}_{859}/0.2112\right)$	B8a	1.8–988	[59]
${\theta }_{TUR}^3$	$-61.251·R_{660}/R_{790}+174.629$	B4,B7	0.21–3.46	[60]
${\theta }_{TUR}^4$	$6834.7·{\tilde{R}}_{821}/{\tilde{R}}_{763}-6632.3$	B8,B7	0–1678	[30]
CDOM
Algorithm	Relation	Corresponding Sentinel-2 bands	Range $\left[m^{-1}\right]$	Reference
${\theta }_{CDOM}^1$	$0.00129+0.6543·\left(R_{670}/R_{490}\right)$	B4,B2	/	[61]
${\theta }_{CDOM}^2$	$ln\left(\left[\left(R_{490}/R_{551}-0.4247\right)/2.453\right]/\left(-13.586\right)\right)$	B2,B3	/	[62]
${\theta }_{CDOM}^3$	$0.2987·{\left(R_{B2}/R_{B4}\right)}^{-1.369}$	B2,B4	0–0.8	[63]
${\theta }_{CDOM}^4$	$5.13·{\left(R_{565}/R_{660}\right)}^{-2.76}$	B3,B4	0.68–11.13	[64]
${\theta }_{CDOM}^5$	$0.133366/{\left(R_{443}/R_{510}\right)}^{2.025}$	B1,B2	0.07–1.1	[65]

).

Selected algorithms in Table 1 are denoted with $\theta$, where the superscript indicates the algorithm number and the subscript identifies the WQ parameter. The performance of all algorithms was tested on UAV HS and UAV MS datasets, denoted with tilda and overbar on $\theta$, respectively. Similarly, we indicated reflectances obtained from the HS sensor by $\tilde{R}$, while the corresponding simulated MS dataset was denoted by $\overline{R}$. Reflectances obtained by the Sentinel-2, i.e., SAT MS dataset are presented by R.

We considered the following: the similarity between the ranges of WQ parameters in the original and in our study, their applicability to the data from UAV HS and Sentinel MS, or their good performance in different study areas based on a literature review. Although the algorithms considered in this analysis were developed for different satellites (MERIS, MODIS, SeaWifs, Landsat) applied in various coastal and inland waters, their performance has already been tested with the Sentinel-2 dataset [30,58,66]. In the literature reviewed, these algorithms were previously used for ranges of WQ parameters closest to ours but still with a wider range, e.g., chlorophyll from 1 to 5 μg/L in Sado Estuary [66] and turbidity between 2 and 6 NTU in Lake Shinji [30]. In addition, we included some algorithms that are not applicable to the MS dataset due to their higher spectral resolution and compatibility with our site (namely ${\theta }_{Chla}^4$, ${\theta }_{Chla}^5$, and ${\theta }_{Chla}^7$). In this study, we focused on simple algorithms from the blue to the NIR region to cover the relevant spectral range (single-band, band ratio, and three-band algorithms).

2.4.2. Proposing New Empirical Algorithms

For the purpose of proposing new empirical algorithms, which are presented in Section 3.3, we tested several regression-based functions against the ratio of two HS bands introducing the independent variable $\rho ={\tilde{R}}_{Mband}/{\tilde{R}}_{Nband}$. Reflectance band ratios are widely used in the WQ algorithms retrieval since they often eliminate background noise and rough water surface interference [67]. Since our datasets are spread over the four seasons, these applied retrieval regression functions bear the characteristics of empirical and site-specific algorithms.

The screening procedure for optimal band ratio consisted of all possible combinations of 273 bands available in the UAV HS dataset filtered through four regression models following logarithmic, rational, power, and exponential functions. Details on choosing the best-performing regression model are presented in Section 3.3.

2.4.3. Accuracy Assessment

To evaluate the comparison between the existing algorithms and the newly proposed ones with in situ data, we have carried out the statistical evaluation based on the coefficient of determination ($R^2$), the root mean square error (RMSE), absolute percentage difference (APD), and relative percentage difference (RPD). Tables with statistical evaluation for all tested algorithms are available in Supplementary Materials (Tables S5–S7). The performance analysis was obtained using MATLAB curve fitting toolbox with least absolute residuals (LAR) method [49,68], where robust nonlinear fit was implemented as derived by [69].

3. Results And Discussion

3.1. In Situ and Hyperspectral Data

Field observations indicate a seasonal pattern of WQ parameters. $Chla$ shows the highest median values during the autumn period (Figure 4), as also previously observed in Kaštela Bay [70]. The lowest CDOM and TUR median concentrations were observed during the dry period in June, while the highest variations were observed in March, during the rainy season, and the largest flow of the Jadro River (Supplementary Materials Table S1).

Spatial and temporal variability of $Chla$, TUR, and CDOM during four sampling campaigns are presented in Figure 5. Different shapes follow a different period of the year, while the intensity of the color is proportional to the quantity of the corresponding parameter.

We utilized reflectances that had been processed using the Savitzky–Golay filter, as shown in corresponding spectral profiles (Figure 6), and are used in the algorithm analysis.

Figure 6 shows the measured reflectance spectra in a wavelength range from 400 nm to 1000 nm, with a shape characteristic of coastal waters [13,71]. The observed valley in the spectra, up to 500 nm, is mainly due to the absorption of CDOM and $Chla$. In particular, a valley near the 420–440 nm is primarily caused by the $Chla$ absorption of blue light. A distinct reflection peak was observed near 550 nm in all spectral profiles, mainly due to phytoplankton’s spectral behavior [72,73]. In addition, the absorption peaks at about 670 nm and the reflection peaks near 690 nm, regarded as fluorescence peaks of $Chla$, are not so pronounced, which could be due to lower TUR and $Chla$ values [13,74].

3.2. Performance of Existing Algorithms

All selected algorithms from Table 1 were tested on the UAV HS and UAV MS datasets. Each algorithm is evaluated at the sampling location using corresponding reflectances from the UAS HS and UAV MS datasets and compared with the in situ measured concentration. A simple linear regression is used to examine match-ups, between the algorithm-derived concentration and measured concentrations for considered WQ parameters, assessed by statistical indicators. Finally, algorithms with the best fit on the UAV HS and UAV MS datasets were additionally applied to the SAT MS dataset, and full comparisons are given in Figure 7, Figure 8 and Figure 9.

3.2.1. Chlorophyll a

A variety of algorithms have been developed to estimate $Chla$ using remote sensing techniques [54,57,75,76,77,78,79,80]. In this study, the evaluation of nine different algorithms for $Chla$ determination was tested on both UAV HS data and UAV MS data, except for a few algorithms that were originally developed for hyperspectral data (namely ${\theta }_{Chla}^4$, ${\theta }_{Chla}^5$ and ${\theta }_{Chla}^7$). Selected algorithms cover the entire spectral range from the blue to the NIR channel. The statistical evaluation for tested algorithms regarding $Chla$ is provided in the Supplementary Materials (Table S4).

The lowest $R^2$ was found for the band-ratio algorithm containing reflectances from blue bands ($R^2$ = 0.05 for ${\tilde{\theta }}_{Chla}^5$). In contrast, the highest value was found for the algorithm that included a red-to-red reflectance ratio ($R^2$ = 0.59 for ${\tilde{\theta }}_{Chla}^7)$ and was determined using the hyperspectral dataset. The second best algorithm also used the hyperspectral dataset ($R^2$ = 0.58 for ${\tilde{\theta }}_{Chla}^4)$ and at the same time was closest to our location. These results confirm that the suitability of both spectral and spatial resolution, as well as site characteristics, plays an important role in evaluating the performance of the algorithm, especially in coastal regions.

As expected, the $R^2$ of the other algorithms using red, red edge, and NIR bands showed similar moderate correlations on both the UAV HS and UAV MS datasets. These observations support the findings that blue–green band-ratio algorithms are more suitable for ocean waters [81,82,83], while in the coastal area, the complexity of spectral properties increase, and optical characteristics are not predominantly determined by phytoplankton biomass as in oceans. In this case, other water constituents, such as CDOM and suspended solids, play an essential role in determining the spectral shape of seawater, highlighting the unsuitability of algorithms developed for ocean waters [84,85]. Most algorithms for $Chla$ determination in coastal areas focus on the red and NIR channels [54,77,86,87], indicating the best optimal applicability, which happened to be the case in our study as well.

$Chla$ in Kaštela Bay was previously analyzed by the RS technology by testing existing algorithms and deriving new ones [13,88]. Both studies revealed better results by applying new algorithms indicating pronounced site-specific behavior of such enclosed water bodies, as reported by several other studies [71,89,90].

Kisevic [13] derived two algorithms for $Chla$ retrieval; one with bands from the blue spectral range (${\theta }_{Chla}^5$) and the other with bands from the red part of spectra (${\theta }_{Chla}^4$). Although the differences were found to be minor, the blue band-ratio algorithm revealed lower statistical errors and was chosen as a more suitable one. During testing of those algorithms, on the UAV HS dataset, we noticed that the correlation of the blue band-ratio algorithm was weak ($R^2$ = 0.05 for ${\tilde{\theta }}_{Chla}^5$), while the second algorithm showed moderate correlation ($R^2$ = 0.58 for ${\tilde{\theta }}_{Chla}^4$). The observed difference in results between these two studies might be due to the more extensive dataset used in our study, which is consistent with the concluding remarks in [13].

For all tested algorithms, the UAV HS data exhibits higher $R^2$ for reflectances from the red and NIR spectral range than from the blue and green spectral regions. In general, all tested algorithms showed a higher correlation using the UAV HS than the UAV MS dataset, which implies how higher spatial and spectral resolution from the UAV HS system enriched obtained information rendering better regression.

The best-fit algorithm applicable to all three datasets (UAV HS, UAV MS, and SAT MS) was ${\theta }_{Chla}^8$, and it is presented in Figure 7. The poor correlation ($R^2$ = 0.08) was observed when applied to satellite data, while moderate correlations were obtained using the UAV MS ($R^2$ = 0.45) and UAV HS ($R^2$ = 0.48) datasets. This indicates the potential importance of using UAV hyperspectral sensors in small coastal zones due to the higher spectral resolution and, consequently, more information provided.

3.2.2. Turbidity

The general recommendation about retrieving TUR in coastal areas is that the red and NIR regions are closely related to TUR levels [27,91,92,93,94]. Different models for TUR determination are proposed using a single-band or band-ratio algorithm [89,92,95]. This study tested four selected algorithms with different band combinations (single band, NIR/red, and NIR/red edge band ratios). Although all tested algorithms exhibit a rather weak $R^2$, the numerically highest correlation was observed using the NIR/red edge band ratio for algorithm ${\tilde{\theta }}_{TUR}^4$ as presented in Figure 8. The algorithm performance details for the TUR parameter are provided in the Supplementary Materials (Table S5).

Previous studies concluded that the reflectances in the red part of the spectrum correlate well with low and medium turbidity values. At the same time, NIR bands were used to identify high turbidity waters [59,96]. Although poorly correlated, our results showed bias toward the red edge and NIR bands, similar to some other studies [27,97]. For instance, comparable in situ TUR values (<4 NTU) were observed in Sent et al. [66], delivering the highest correlation to the NIR band. Furthermore, the obtained results support the findings about the NIR region being equally sensitive and less influenced by bottom reflectance than the red region in shallow waters [91,94].

Even within the very low correlations, the tested algorithms showed relatively better performance when derived from the hyperspectral dataset. For all practical purposes, all these algorithms showed poor correlations for our case study, rendering their prediction capabilities questionable. However, these results indicate again the necessity for locally derived algorithms, which we further investigate in Section 3.3.

3.2.3. Colored Dissolved Organic Matter

Since the CDOM is optically visible and can be monitored by remote sensing, the retrieval models have been proposed for CDOM determination in the ocean and coastal environments [33,98]. Furthermore, several band-ratio algorithms have been proposed for CDOM absorption, and in this study, five algorithms were tested on our datasets. Similarly to the TUR algorithm testing, the selected algorithms for CDOM also showed poor correlation. According to the statistical evaluation (details available in Supplementary Materials, Table S6), the numerically best algorithm was ${\tilde{\theta }}_{CDOM}^2$, which uses the reflectances from the blue and green bands, a combination of bands commonly used for CDOM retrieval in coastal areas.

Similarly to the previously observed results, all algorithms derived from the UAV HS dataset performed better than those derived from the UAV MS dataset. The ${\tilde{\theta }}_{CDOM}^2$ algorithm delivered the highest statistical indicators for the UAV HS dataset ($R^2$ = 0.12) with rather weak performance for the UAV MS ($R^2$ = 0.06) and satellite MS datasets ($R^2$ = 0.03), as presented in Figure 9.

3.3. Proposed New Regression-Based Algorithms

The above analysis of selected algorithms has shown that their application has potential limitations in oligotrophic waters, which are characterized by a very narrow range of values for all of the considered WQ parameters. However, to take advantage of the use of UAV-based hyperspectral sensors, we tested and developed new site-specific models tailored to the oligotrophic coastal water body at hand. The necessity of using the UAV HS dataset as an intermediary between ground truthing and satellite data in such a complex and dynamic ecosystem is also addressed by some studies [9,96,99] and is supported by our study as well.

In this section, we propose four empirical (regression-based) models which are built upon the statistical relationships between the UAV HS dataset and the measured concentration of WQ parameters considered. The empirical models tested are logarithmic, rational, power, and exponential regression-based models denoted as ${\tilde{\mathrm{\Psi }}}_x^i$, where x stands for the WQ parameter and i denotes logarithmic, rational, power, and exponential regression models with $i=1,\dots ,4$, respectively.

$$\begin{array}{ccc}\hfill {\tilde{\mathrm{\Psi }}}_x^1& =& a+bln\left(\rho \right)\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill {\tilde{\mathrm{\Psi }}}_x^2& =& \left(a\rho +b\right)/\left(\rho +c\right)\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill {\tilde{\mathrm{\Psi }}}_x^3& =& a{\rho }^b+c\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill {\tilde{\mathrm{\Psi }}}_x^4& =& aexpb\rho \hfill \end{array}$$

(4)

where $a,b,$ and c are fitting parameters for each considered algorithm and $\rho ={\tilde{R}}_{Mband}/{\tilde{R}}_{Nband}$ is a band-ratio independent variable to be screened through all possible combinations of 273 bands available in the UAV HS dataset. Visual presentations of $R^2$ performance for various models and band combinations, as well as the implemented robust nonlinear fitting setup, can be found in the Supplementary Materials (Figures S3–S5 and Figures S6–S8, respectively).

3.3.1. Chlorophyll a

All tested regression models for $Chla$ retrieval show a very high correlation when using the ${\tilde{R}}_{673}/{\tilde{R}}_{664}$ band ratio (Table S8 in Supplementary Materials). Numerically, the best-fit algorithm was the exponential regression model ${\tilde{\mathrm{\Psi }}}_{Chla}^4$ (although the logarithm regression model showed almost the same performance). This algorithm used bands from the red part of the spectra and showed a high correlation with in situ $Chla$ values ($R^2$ = 0.958) (Figure 10).

The performance of different regression models (details can be found in the Supplementary Materials) is relatively similar and uses the same hyperspectral wavelength bands, which is somewhat expected for homogenous in situ data with low variability.

An adequately performing algorithm for Kaštela Bay by Kisevic [13] was also found in the red area, using a portable spectrophotometer, and the obtained bands are close to the ones we discovered. Conversely, most satellite-based algorithms used band ratio red/red edge, among other reasons, for the lack of distinctive spectral resolution.

3.3.2. Turbidity

The highest correlation between the reflectance band ratio and in situ TUR concentrations was found using the red edge part of spectra ${\tilde{R}}_{730}/{\tilde{R}}_{777}$ with logarithmic regression model ${\tilde{\mathrm{\Psi }}}_{TUR}^1$. The statistical summary for all tested algorithms is available in the Supplementary Materials (Table S9). In Figure 11, the model performance is presented and a high correlation ($R^2$ = 0.92) is in line with findings by [66,89,96], who also noticed the importance of using the red edge part in TUR retrieval.

3.3.3. Colored Dissolved Organic Matter

The best-fit algorithm for CDOM retrieval was the exponential model using the band ratio from the violet part of the spectra. In Figure 12, the regression model ${\tilde{\mathrm{\Psi }}}_{CDOM}^4$ is displayed with a high correlation $R^2=0.92$. The performance of all tested algorithms is available in the Supplementary Materials (Table S10). Most of the algorithms for CDOM retrieval used blue/red or green/red band ratios, although CDOM absorption decreases exponentially with increasing wavelength [31,64]. The reason for that is the problematic atmospheric correction of the blue sensor [31,64] and excessive absorption by CDOM and low natural water-leaving radiance at low wavelengths, reducing the usable signal [31,64,100]. In our case, a UAV does not need an atmospheric correction, and this could explain why the algorithms using shorter wavelengths performed well in our study. This hypothesis should be considered cautiously, and more studies of CDOM retrieval using UAVs or point spectrophotometers should be conducted. Also, for the CDOM retrieval, a benefit of incorporating UV spectral range [26,72] has been noticed, implicating the importance of using sensors with higher spectral range and resolution compared to satellites.

3.3.4. Mapping Application

Following the proposed new local-regression-based algorithm for $Chla$ (Figure 10), we obtained the WQ parameter map for the area that was recorded by the UAV. For presentation purposes, in Figure 13 we present $Chla$ map based on the hyperspectral data retrieved by UAV flights for all field campaigns covering four seasons. As expected, during the late spring and summer season, the area covered by UAV HS exhibits larger $Chla$ concentrations compared to the winter season. Interestingly, $Chla$ pattern in all seasons shows the trend probably resulting from the exchange flow pattern within the estuary (net seaward flow and its counterpart as a density-driven flow landward). The exchange flow dynamics in the Jadro River case are relatively fast and create different timing between the process of converting dissolved nutrients into the new biomass and the freshwater residence time [101]. The work performed represents a valuable application in the field of environmental monitoring, as it provides a high-resolution and efficient method for mapping and assessing water quality parameters, particularly chlorophyll-a concentration, across different seasons. These results can support long-term decision-making processes related to coastal management, ecological assessments, and the development of targeted remediation strategies, ultimately contributing to the preservation and restoration of coastal water ecosystems.

4. Conclusions and Future Remarks

Remote sensing has become a crucial tool for water quality monitoring and is evolving rapidly. This study sought to explore the efficiency of existing algorithm, commonly used for water quality parameter retrieval from readily available satellite imagery and contrast them with data from UAVs. Despite being less accessible, UAV-derived hyperspectral data offer superior spectral and spatial resolution. Our findings echo previous studies conducted in coastal waters, where existing algorithms demonstrate a subpar correlation due to the unique site-specific characteristics, underscoring the need for locally calibrated algorithms. This shortcoming is especially noticeable in our case study given its oligotrophic nature. While some algorithms performed poorly, the hyperspectral data generally outperformed both simulated and satellite-obtained multispectral data.

The lackluster performance of existing models motivated us to develop new, regression-based algorithms tailored to our study area, utilizing UAV hyperspectral data and in situ observations. The resulting algorithms—devised from an array of 273 spectral bands available with the utilized UAV—boasted high correlation with in situ data. This reiterates the necessity for site-specific algorithms and comparable monitoring scales. The efficacy of our monitoring approach underlines the potential of UAV-based water quality mapping as a swift and efficient screening tool in oligotrophic coastal waters.

Hyperspectral data bring with them the challenge of demanding data processing and redundancy issues. However, the relentless advancement in UAV technology and sensor versatility render this data source invaluable for complex monitoring needs, such as sudden algal blooms or spills. In cases where simple prediction models are already established, UAVs can be quickly deployed to update site-specific algorithms with a growing pool of data. Despite the necessity for numerous field campaigns to account for different water quality parameters and seasonal variations, the escalating availability of UAVs and corresponding sensors suggests this is a viable next step in monitoring development.

Our approach offers a more accurate, locally-focused alternative to the pursuit of universally applicable models. UAV-derived data bridge the spatial scale gap between in situ and satellite data and signify a promising area of future focus. It is important to acknowledge that this interplay between disparate spatial scales necessitates particular consideration in future monitoring endeavors. Our work accentuates the scientific value and contributions of high-resolution, locally calibrated algorithms for water quality assessment while also acknowledging the challenges and outlining future advancements in this dynamic field.