Deriving Nutrient Concentrations from Sentinel-3 OLCI Data in North-Eastern Baltic Sea

Abstract

Nutrients are important elements in marine ecosystems and water quality, and have a major role in the eutrophication of water bodies. Monitoring nutrient loads is especially important for the Baltic Sea, which is especially sensitive to the eutrophication. Using optical remote sensing data in mapping total nitrogen (TN) and total phosphorus (TP) is challenging because these substances do not have a direct influence on the water optics that remote sensing sensors can detect. On the other hand, it would be very rewarding. In this study, more than 25,000 Sentinel-3 Ocean and Land Colour Instrument (OLCI) data algorithms were tested in order to detect the TN and TP concentrations in the Estonian marine waters between 2016–2021. The TN estimations were well derived for Estonian marine waters (R² = 0.73, RMSE = 4.87 µmolN L⁻¹, MAPE = 14%, n = 708), while the TP estimations were weaker (R² = 0.38, RMSE = 0.23 µmolP L⁻¹, MAPE = 24%, n = 730). The Estonian marine waters were divided into six geographic regions in order to study the effect of regional water quality on the TN and TP retrievals. The nutrient concentrations were derived in every region when spring and summer periods were treated separately. In this study, the detection of both nutrients was more successful in more closed areas with P deficiency, while in open sea areas it was more challenging. This study shows that it is possible to estimate nutrients, especially TN, from remote sensing data. Consequently, remote sensing could provide a reliable support to the conventional monitoring by covering large marine areas with high temporal and spatial resolution data.

1. Introduction

Baltic Sea is the world’s largest inland brackish water sea and is very well studied. Eutrophication has been evident in the Baltic Sea for many decades, due to past high and still excessive loads of total nitrogen (TN) and total phosphorus (TP) [1]. The combination of a large catchment area with a high rate of human activities and a small body of water with limited exchange with the Atlantic Ocean through the narrow and shallow Skagerrak makes the Baltic Sea very sensitive to nutrient enrichment and eutrophication. Therefore, in the Baltic Sea, the large input of nutrients like phosphorus and nitrogen is a major environmental concern [2]. The conventional in situ water quality monitoring has failed to characterize nutrient dynamics because of the limitations in spatial sampling and poor availability of reliable data for nutrient loads [3,4].

Nutrient concentrations are very important because they cause eutrophication—an increase in aquatic biomass. Phytoplankton biomass is usually characterised by chlorophyll-a (Chl-a) concentration. Therefore, the Chl-a is the primary indicator of the waterbody ecological state [4,5]. In the Baltic Sea, the open basins are mainly nitrogen limited, especially in spring, when the spring phytoplankton bloom is peaking. In summer, massive blooms of nitrogen fixing cyanobacteria occur in the Baltic Sea. They are driven by excess phosphorus, along with high temperatures. The cyanobacterial blooms contribute significant quantities of new nitrogen to the pelagic ecosystem, hence reducing phosphorus loads is very important to reduce those blooms [4]. To develop an effective nutrient management strategy, better understanding of the Baltic Sea ecosystem is required [1].

Although nutrient concentrations are important water quality indicators in coastal waters, and at the same time the Baltic Sea is often claimed as the most studied sea in the world, the quality of regional monitoring of the nutrient pollution entering the area remains relatively poor [4,6]. Remote sensing has the abilities that could be highly beneficial to marine monitoring—it provides high spatial and temporal resolution which is impossible to achieve with in situ measurements. This is beneficial even if remote sensing cannot provide similar accuracy as the time consuming and expensive laboratory methods. Consequently, combining both conventional in situ sampling and remote sensing methods should be the most optimal way to study the marine environment, provided the remote sensing methods can provide sufficient accuracy.

The European Union’s Earth observation programme, the Copernicus program [7], has launched a mission Sentinel-3, which is a constellation of two satellites (A and B, launched in 2016 and 2018, respectively) [8]. Sentinel-3 has a medium resolution (300 m) Ocean and Land Colour Instrument (OLCI) onboard for marine and land research. It has 21 spectral bands and provides global coverage (at the equator) every two days. OLCI was built for marine monitoring and has well placed spectral bands for that purpose (ESA Sentinel Online).

Remote sensing has been widely used for water quality monitoring [9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24], because of its large advantage in temporal and spatial coverage compared with in situ methods. So far, most studies have focused on water quality variables, such as Chl-a, turbidity or transparency (measured with Secchi disk depth, SD), total suspended solids (TSS), and coloured dissolved organic matter (CDOM), which all are optically active substances. Remote sensing presents a challenge in estimating nutrients like TN and TP concentration in water, because they have no direct optical properties and spectral characteristics [25]. However, nutrients like TN and TP can be highly correlated with optically active variables that can be estimated by remote sensing [22,26]. When nutrients concentrations are in good correlation with some optical properties of water, the optical parameter can be used as a proxy for TN and/or TP. The previous work testing the retrieval of the optical parameters in the Baltic Sea based on OLCI data with Case 2 Regional CoastColour (C2RCC) atmospheric correction have given good results [27]. However, some studies have found only the band ratios based on the reflectances to work well, and not the Level-2 products themselves [11,17].

Remote sensing techniques have been used to estimate different forms (total, inorganic, etc.) of nitrogen and phosphorus in freshwater bodies [18,26,28,29,30,31,32,33,34,35,36,37,38]. Although less studies have been carried out in sea water [3,5,39,40], the results have been encouraging with errors as low as 11% for TN and 13% for TP [3]. While most studies in marine waters are focusing on detecting TN or/and other forms of nitrogen, there are only a few attempting to detect TP [40]. The chosen methods are either the neural network [3] and machine learning approaches [5,40], or multiple stepwise [39] regression models.

The main objective of current work was to find the best algorithms for estimating the concentrations of TN and TP in Estonian marine waters using Sentinel-3 OLCI data. In addition, the temporal and spatial variability of the best algorithms were studied to estimate nutrients with the highest possible accuracy. Accurate algorithms would give the possibility to monitor the nutrients over large spatial scales and with high temporally frequency, which is not possible with conventional methods. This is the first time where remote sensing with regionally tuned algorithms for nutrient estimations have been used in the north-eastern Baltic Sea.

2. Materials and Methods

2.1. In Situ Dataset and Study Area

In the current study, the measurements from the Estonian National Monitoring Program were used. In total, 2103 TN and 2078 TP measurements were carried out during 2016–2021 in Estonian marine waters (approx. 36,500 km² (~10%) area of the Baltic Sea). Water samples for both nutrients were collected from 1 m depth. Unfortunately, the optical parameters, like the Chl-a, SD, TSS and the absorption of CDOM at 400 nm (a_CDOM) were not always measured together with the nutrients. The database included 1387 Chl-a, 1162 SD, 49 TSS, and only 39 a_CDOM measurements.

To eliminate the possible effect from the coast or the seafloor on the remote sensing signal, 29 sampling stations (382 TN and 378 TP measurements) that were closer than 1 km from the shore were removed from the analysis. After this elimination, the database contained 63 sampling stations with 1721 TN and 1700 TP measurements in total. All the stations used in the study are shown in Figure 1.

In this study, the Estonian marine waters were divided into six geographical regions (R1–R6; Figure 1) based on the § 3²⁴(4) of the Estonian Water Act [41], because there are large differences in the characteristics of the regions. Division used in the current study is similar to the state classification of the coastal waters in the Estonian Water Act, with the difference being that we included the territorial waters and the exclusive economic zone (EEZ) of the six regions together with the internal waters. The main characteristics of the six regions are:

• R1: south-east area of the Gulf of Finland. It is oligohaline (2.5–6 ppt) open water with measured max Chl-a 29.5 mg m⁻³ and min SD 0.8 m during 2016–2021. The largest inlet in region R1 is the Narva River. Twelve sampling stations were in this region.
• R2: Pärnu Bay is located in the north-eastern part of the Gulf of Riga. The bay is semi-enclosed, oligohaline (4.0–5.5 ppt), with a large inlet of nutrient rich Pärnu river. During 2016–2021, the measured max Chl-a was 45 mg m⁻³ and min SD 0.4 m. It is the smallest region by area (411 km²), and the max depth in the mouth of the bay is 12 m. Four sampling stations were included in this region.
• R3: western part of Gulf of Finland, it has mesohaline (4.5–6.5 ppt) and deep water with measured max Chl-a 25.3 mg m⁻³ and min SD 2 m during 2016–2021. Sixteen sampling stations were in this region.
• R4: Baltic Proper area of the West Estonian archipelago, open sea, with mesohaline (6–7 ppt) water. Region R4 is a shallow area open to waves with measured max Chl-a 15.6 mg m⁻³ and min SD 3 m during 2016–2021. Seventeen sampling stations were in this region.
• R5: Väinameri Sea or the Sea of Straits (2200 km²); it has mesohaline (3–6.5 ppt) unstratified water, and it is a shallow, concealed area with measured max Chl-a 6.7 mg m⁻³ and min SD 1.2 m during 2016–2021. Seven sampling stations were in this region.
• R6: north half of Gulf of Riga with mesohaline (4–6 ppt), shallow, sheltered and seasonally stratified waters with measured max Chl-a 71.8 mg m⁻³ and min SD 0.5 m during 2016–2021. Seven sampling stations were in this region.

2.2. In Situ Parameters

The TN and TP concentrations were measured by the accredited Estonian Marine Institute’s laboratory with a continuous flow automated wet chemistry analyser Skalar SANplus (Skalar Analytic B.V.,DeBreda, The Netherlands), using the standard methods EN ISO 11905-1 and EN ISO 15681-2. Detection limits were 10 µg L⁻¹ for both TN and TP, and measurement uncertainty did not exceed 25%.

Temperature and salinity data were extracted from the CTD profiles, measured with the Sea and Sun Technology M90 or SAIV SD201 probes. The Chl-a concentration in a sample was determined spectrophotometrically (light absorption) using ISO 10260 standard and HELCOM guidelines [42]. Water transparency (SD) was measured with a standard Secchi disk.

Absorption by CDOM, a_CDOM, was measured spectrophotometrically with a PERKIN ELMER Lambda 35 UV/VIS spectrometer in the range 350–750 nm from a filtered (Millex 0.22 μm) water sample in a 10 cm cuvette against distilled water. The measurements were corrected for residual scattering according to Davies-Colley and Vant [43], and a final calculation for a_CDOM at 400 nm was conducted by Lindell et al. [44].

2.3. Sentinel-3 Dataset

The processing of the Copernicus Sentinel-3 OLCI Level-1 (300 × 300 m resolution) data was done using the Estonian national satellite data centre portal for Earth observation data processing [45]. This platform has Sentinel-3 OLCI Level-1 archive and atmospheric correction tools. For atmospheric correction, the C2RCC v.1.5 processor [46] with the multisensor pixel identification tool (IdePix) were used. In extracting the match-ups, the same day 1 × 1 pixels with the sampling stations date and location were used. The 1 × 1 pixel size was used because a single OLCI pixel is large, and many in situ locations were close to each other. Pixels with cloud or other quality flags (except the “quality_flags.sun_glint_risk” flag) were removed from further analyses.

In total, 741 same-day cloud-free match ups (including Sentinel-3A and -3B) during 2016–2021 were found for TP and 719 for TN. The most match ups were from May to July (

Table 1. The distribution of the match-ups for total phosphorus (TP) and nitrogen (TN) between study regions (R1–R6) and months (April–November 2016–2021). In the case that there was a different number of match ups for TN, it is shown in parentheses.

TP (TN)
	R1	R2	R3	R4	R5	R6	Total
April	20	10 (5)	33 (29)	18	4	16 (12)	101
May	22	33 (28)	49	25	5	32 (28)	166
June	31	12	56	21	13	23	156
July	37	18	44	33	20	34	186
August	25	14	23	18	5	8	93
September	4	10	3	2	0	9	28
October	4	0	0	3	0	2	9
November	0	0	1	1	0	0	2
Total	143	97 (87)	209 (205)	121	47	125 (117)	741 (719)

). There were no cloud-, ice- or snow-free match-ups in January, February, March and December.

2.4. TN and TP Retrieval Methods

Fifteen different formulas were used with the atmospherically corrected angular dependent water-leaving reflectance (C2RCC reflectances) and the C2RCC processor Level-2 (L2) products for the retrieval of TN and TP (

Table 2. The general formulas used in this study. B indicates the atmospherically corrected angular dependent water-leaving reflectance band, and index a, b, or c indicates different Ocean and Land Colour Instrument (OLCI) bands (15 bands in different options) or the Level-2 product.

Formula
1. Ba + Bb
2. Ba − Bb
3. Ba/Bb
4. Ba * Bb
5. Ba + Bb + Bc
6. Ba + Bb * Bc
7. (Ba + Bb) * Bc
8. (Ba − Bb) * Bc
9. (Ba + Bb)/Bc
10. Ba * Bb/Bc
11. (Ba − Bb)/(Ba + Bb)
12. (Ba/Bb) * (Ba/Bb)
13. Ba/Bb − Ba/Bc
14. Ba − (Bb + Bc)/2
15. Ba/(Bb + Bc)

). Only the general formulas used in this study are shown in the Table 2. Every formula was tested with different combinations by altering C2RCC reflectances on the 15 bands, and the 13 different C2RCC L2 products (

Table 3. Sentinel-3 OLCI spectral bands, their central wavelengths for Case 2 Regional CoastColour (C2RCC) reflectances, and C2RCC Level-2 (L2) products [46].

Band	Centre (nm)	L2 Product	L2 Product Description
1	400	iop_apig	Absorption coefficient of phytoplankton pigments at 443 nm (m−1)
2	412.5	iop_adet	Absorption coefficient of detritus at 443 nm (m−1)
3	442.5	iop_agelb	Absorption coefficient of coloured dissolved organic matter (CDOM) at 443 nm (m−1)
4	490	iop_bpart	Scattering coefficient of marine particles at 443 nm (m−1)
5	510	iop_bwit	Scattering coefficient of white particles at 443 nm (m−1)
6	560	iop_adg	Detritus + CDOM absorption at 443 nm (m−1)
7	620	iop_atot	Phytoplankton + detritus + CDOM absorption at 443 nm (m−1)
8	665	iop_btot	Total particle scattering at 443 nm (m−1)
9	673.75	kd489	Irradiance attenuation coefficient (Kd) at 489 nm (m−1)
10	681.25	kdmin	Mean Kd at the three bands with minimum Kd (m−1)
11	708.75	kd_z90max	Depth where 90% of the water-leaving irradiance comes from (m−1)
12	753.75	conc_tsm	TSS dry weight concentration (g m−3)
16	778.75	conc_chl	Chl-a concentration (µg L−1)
17	865
18	885

) were included only to the simplest formulas (1–4 in Table 2). In total, 25,013 different algorithms were made for testing.

To find the best algorithm for deriving TN and TP, a quadratic polynomial regression model was used similarly by Huang et al. [35] in Xiangxi Bay, China. Equation (1) was used to derive TN or TP concentrations.

y = ax² + bx + c

(1)

where y is the derived TN or TP, x is the algorithm of L2 product or band ratios, and a, b, and c are the polynomial regression coefficients.

The entire match up database was tested with over 25,000 algorithms, but also every region (R1–R6) separately was tested on six different temporal divisions:

• April to November;
• May to September;
• April to May;
• June to September;
• April to June;
• July to September.

The first period includes all data. The second period (May to September) is agreed in the monitoring program as the biologically active period. For example, the indicator for eutrophication in the monitoring program is mean Chl-a for the May–September period. The third and fourth periods are the traditional way to divide spring and summer between months (April to May and June to September, accordingly). Lastly, the fifth and sixth periods are the new way to define spring and summer seasons, to distinguish the spring bloom from the summer cyanobacteria bloom (April to June and July to September, accordingly).

In addition, three combined regions were formed: the Gulf of Finland (R1 + R3), west Estonian archipelago waters and the Baltic Proper (R4 + R5), and the Gulf of Riga (R2 + R6). All of them were used to test all the algorithms on the six temporal divisions.

2.5. Statistical Analysis

Mean values were used to describe in situ characteristics. The mean is the average value of the dataset.

The performance of algorithms was evaluated using the determination coefficient (R²). R² is used to analyze how well observed in situ values are predicted by the model based on the proportion of total variation of outcomes explained by the model; it is calculated using Equation (2).

$$R^2=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\widehat{y}\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(2)

where, $\widehat{y}$ is the predicted value, y is the observed value, $\overline{y}$ is the mean value of observed y values, and n is the number of observations.

To evaluate the error of the derived nutrients, root mean squared error (RMSE) and bias were used (both are in the same units as the quantity being estimated). RMSE is a frequently used measure of differences between values observed in situ and predicted by a model; it is calculated using Equation (3). Bias shows a systematic error, and it is calculated using Equation (4).

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(\widehat{y}-y_i\right)}^2}$$

(3)

$$\mathrm{bias}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left(\widehat{y}-y_i\right)$$

(4)

where, $\widehat{y}$ is the predicted value, y is the observed value, and n is the number of observations.

In addition, the performance of algorithms was evaluated using accuracy through the mean absolute percentage error (MAPE). MAPE measures the percentage error of the estimated values in relation to the actual values; it is calculated using Equation (5).

$$\mathrm{MAPE}=\frac{100\%}{n}{\displaystyle \sum }_{i=1}^n\left|\frac{y_i-{\widehat{y}}_i}{y_i}\right|$$

(5)

where, $\widehat{y}$ is the predicted value, y is the observed value, and n is the number of observations. Accuracy represents correctly predicted values (expressed as a percentage). Accuracy is calculated as $100-\mathrm{MAPE}$.

The 95% confidence interval value (95%CI) is used to describe the uncertainties associated with the quadratic polynomial regression model, and it is presented for the coefficients of the regression model. It is calculated using Equations (6) and (7), where se_fit is the unbiased estimate of the standard error of the fit, and t_crit is the critical value of the t distribution with degrees of freedom of the residuals with significance level 0.05/2.

$$\mathrm{CI}=\pm s{e}_{fit}\times t_{crit}$$

(6)

$${se}_{\mathit{fit}}=\sqrt{\frac{s{s}_{res}}{d{f}_{res}}\times X_0^T{(X^{T}X)}^{-1}{X}_0}$$

(7)

where, ss_res is the sum of squares for the residuals, df_res is the degrees of freedom for the residuals, X₀ is the column array and X is the design matrix of the observed values, and ^T notes the transposed matrix.

3. Results

3.1. Match-Up In Situ Database

The means of the in situ TN, TP, TN:TP, Chl-a and SD in different regions and seasons of the match-ups stations are shown in

Table 4. The number of in situ unique sampling stations in different Estonian coastal regions (R1–R6), the in situ mean values for total nitrogen (TN), total phosphorus (TP), TN:TP ratio, chlorophyll-a concentration (Chl-a), and Secchi disk depth (SD) in April to September, spring (April to June), and summer (July to September) of the match ups during 2016–2021. The number of measurements is shown in the parentheses.

	Unique Stations	TN µmolN L−1	TP µmolP L−1	TN:TP	Chl-a mg m−3	SD m
R1	12	23.2 (139)	0.82 (139)	32.5 (139)	6.8 (105)	3.3 (113)
Spring		22.5 (73)	0.88 (73)	30.4 (73)	7.8 (59)	3.4 (60)
Summer		23.9 (66)	0.75 (66)	34.8 (66)	5.6 (46)	3.2 (53)
R2	4	37.6 (87)	0.61 (97)	65 (87)	8.9 (78)	1.2 (93)
Spring		45.1 (45)	0.57 (55)	80.9 (45)	10.7 (42)	1.3 (52)
Summer		29.6 (42)	0.66 (42)	47.9 (42)	6.7 (36)	1.1 (41)
R3	16	20.9 (204)	0.74 (208)	31.8 (204)	6.4 (205)	4.7 (170)
Spring		20.6 (135)	0.78 (138)	29.8 (134)	6.9 (136)	5.3 (115)
Summer		21.5 (70)	0.64 (70)	35.7 (70)	5.3 (69)	3.7 (55)
R4	17	18.6 (117)	0.56 (117)	35.7 (117)	4.5 (106)	5.8 (89)
Spring		17.6 (64)	0.59 (64)	32.7 (64)	5.0 (57)	6.6 (49)
Summer		19.9 (53)	0.53 (53)	39.3 (53)	4.0 (49)	4.7 (40)
R5	7	20.6 (47)	0.44 (47)	58.5 (47)	2.1 (26)	4.1 (41)
Spring		21.8 (22)	0.38 (22)	76.8 (22)	1.2 (13)	4.5 (21)
Summer		19.6 (25)	0.49 (25)	42.3 (25)	2.9 (12)	3.8 (20)
R6	7	25.9 (114)	0.62 (122)	45.6 (114)	7.2 (121)	2.4 (95)
Spring		26.8 (63)	0.63 (71)	46.0 (63)	9.0 (70)	2.4 (54)
Summer		24.9 (51)	0.60 (51)	45.2 (51)	4.8 (51)	2.4 (41)
R1–R6	63	23.8 (708)	0.67 (730)	40.7 (708)	6.4 (641)	3.7 (601)
Spring		24.3 (401)	0.70 (423)	41.2 (401)	7.4 (377)	4.0 (351)
Summer		23.3 (307)	0.63 (307)	39.9 (307)	5.1 (264)	3.1 (250)

. The means of the nutrients, and especially the TN:TP ratios, are not very different between the seasons, except in some regions like R2 and R5. The SD in region R2 is significantly lower than in other regions. The Chl-a in spring is the highest in regions R2 and R6 (Table 4).

The strongest relationships between in situ Chl-a and in situ TN were in the region R1 summer season, and in the rest of the regions no remarkable relationship was shown. The TP did also not have any relationship with Chl-a, except in the region R4 in spring (Figure 2A). The transparency and TP had the strongest relationship in only region R5 (spring and summer) and in the summer periods of regions R2 and R4 (

Table 5. Determination coefficients (R²) of the linear regression, count (n) and the p-value (p) between total nitrogen (TN, µmolN L⁻¹) or total phosphorus (TP, µmolP L⁻¹) and chlorophyll-a concentration (Chl-a, mg m⁻³) or Secchi disk transparency (SD, m) in different regions and seasons: spring (April to June) and summer (July to September) of the in situ match-up dataset. Cases where p-value is > 0.05 is marked with red colour.

		Chl-a						SD
		TN			TP			TN			TP
		R2	n	p	R2	n	p	R2	n	p	R2	n	p
All dataset	Spring	0.05	355	0.000	0.07	377	0.000	0.27	329	0.000	0.00	351	0.685
	Summer	0.06	264	0.000	0.12	264	0.000	0.20	250	0.000	0.06	250	0.000
R1	Spring	0.18	59	0.001	0.08	59	0.035	0.14	60	0.003	0.00	60	0.858
	Summer	0.30	46	0.000	0.11	46	0.021	0.05	53	0.006	0.02	53	0.367
R2	Spring	0.03	32	0.307	0.01	42	0.646	0.01	42	0.651	0.10	52	0.024
	Summer	0.07	36	0.108	0.07	36	0.124	0.15	41	0.013	0.24	41	0.001
R3	Spring	0.04	132	0.021	0.18	136	0.000	0.02	111	0.123	0.10	115	0.001
	Summer	0.04	69	0.094	0.01	69	0.375	0.10	55	0.020	0.00	55	0.700
R4	Spring	0.00	57	0.768	0.36	57	0.000	0.01	49	0.418	0.00	49	0.828
	Summer	0.03	49	0.233	0.11	49	0.019	0.03	40	0.269	0.26	40	0.001
R5	Spring	0.01	13	0.717	0.11	13	0.280	0.04	21	0.385	0.24	21	0.025
	Summer	0.30	13	0.054	0.00	13	0.999	0.09	20	0.208	0.38	20	0.004
R6	Spring	0.00	62	0.792	0.00	70	0.633	0.11	46	0.022	0.02	54	0.373
	Summer	0.04	51	0.164	0.07	51	0.053	0.01	41	0.520	0.16	41	0.009

). The relationship between in situ TN and SD seems to be exponential rather than linear. For example, in spring period the R² was 0.27 with linear regression, but 0.39 with exponential regression (Table 5; Figure 2B). All the linear regressions between the in situ TN and TP and the in situ optically active substances within the match-up dataset according to the regions and seasons are shown in Table 5.

The match-up in situ database showed strong dependence of TN with a_CDOM (R² = 0.78, n = 21, p-value = 1E−07), but only a weak relationship between TP and a_CDOM (R² = 0.19, n = 21, p-value = 0.05) (Figure 3). Within the match-up in situ dataset, the relationships with TSS were not significant (TP: R² = 0.13, n = 18, p-value = 0.15 and TN: R² = 0.09, n = 18, p-value = 0.22), though with a larger dataset some relationships might occur.

3.2. Total Nitrogen

After testing over 25,000 band/product ratio combinations on the 719 TN match-ups, the highest R² for each region (R1–R6) and each temporal division is shown in Table 4 with the combinations of regions and with no regions results. The highest R² was usually in the shortest period: April to May. Although the highest R² was in April to May, the combination of April to June + July to September periods was chosen in the current study. There are two reasons for that: (1) the algorithms worked much better in all regions in the July to September period than in the June to September period; (2) too low count of match-ups in the April to May period in R5 (n = 9) made the polynomial regression results unreliable (

Table 6. The number of match-ups (n) and the determination coefficients (R²) of the quadratic polynomial regression model for the best algorithm out of 25,013 for different regions (R1–R6) or region combinations for deriving total nitrogen (TN).

	April–Nov		May–Sept		April–May		June–Sept		April–June		July–Sept
	n	R2	n	R2	n	R2	n	R2	n	R2	n	R2
No regions	719	0.49	620	0.50	245	0.64	463	0.31	401	0.53	307	0.44
R1	143	0.17	119	0.14	42	0.46	97	0.17	73	0.31	66	0.21
R2	87	0.66	82	0.68	33	0.55	54	0.65	45	0.66	42	0.71
R3	205	0.14	175	0.09	78	0.18	126	0.14	134	0.15	70	0.19
R4	121	0.10	99	0.14	43	0.32	74	0.15	64	0.13	53	0.22
R5	47	0.16	43	0.15	9	0.68	38	0.22	22	0.19	25	0.51
R6	116	0.08	102	0.11	40	0.20	74	0.12	63	0.12	51	0.17
R1 + R3	348	0.16	294	0.12	120	0.30	223	0.13	207	0.22	136	0.15
R2 + R6	203	0.59	184	0.62	73	0.60	128	0.46	108	0.61	93	0.61
R4 + R5	168	0.09	142	0.08	52	0.20	112	0.09	86	0.18	78	0.16

).

Generally, the option with no regions applied was performing rather well in all temporal divisions (R² ranged between 0.31–0.64). In comparison of different regions, the algorithms were performing very well in R2. In the R3 and R6 regions, the performance was the weakest along all the temporal divisions (Table 6). The combination of regions did not improve the results and it was not used in the further study. Not using any regions might be a simple option for fast results, but using regions improves the MAPE of the deriving TN for the entire vegetation season (April to September) from 18.2% to 14.2% (Figure 4). The overall results in the different regions not applied and all regions combined versions for the entire vegetation period were relatively good: R² from 0.52 to 0.73, with RMSE from 6.5 to 4.9 µmolN L⁻¹, respectively (Figure 2).

There were significant differences in deriving TN in different regions. The strongest relationship with TN was in region R2 (Pärnu Bay) (R² = 0.74), while the weakest relationship was in regions R3 and R6 (R² = 0.17 and 0.15, respectively) (Figure 5). The RMSE was between 3.11 to 8.85 µmolN L⁻¹ (regions R4 to R2, respectively). Although the region R2 had the strongest relationship towards TN, the accuracy was the lowest (81.4%), because the absolute values of TN were at least twice as high compared with other regions. In region R2, the TN max is 105.6 µmolN L⁻¹, while in other regions, the TN max ranges from 29.8 µmolN L⁻¹ to 48.6 µmolN L⁻¹.

The statistics of the spring and summer seasons separately are shown in

Table 7. The best formulas (B in the formulas indicates the C2RCC reflectance band, and the number indicates the band number), the determination coefficient (R²), p-value, mean and the root mean square error (RMSE) (both in µmolN L⁻¹), the mean absolute percentage error (MAPE), and bias (µmolN L⁻¹) for each region in spring (April to June) and summer (July to September) for the derived total nitrogen (TN).

Region	Season	Formula	R2	p-Value	Mean	RMSE	MAPE	Bias
No regions	Spring	B16 * B17/B5	0.53	<0.001	24.3	7.6	19.9	0.0
	Summer	B18/B4 − B18/B5	0.44	<0.001	23.2	4.5	16.0	0.0
R1	Spring	(B8 − B10) * B17	0.31	<0.001	23.2	3.9	14.3	0.0
	Summer	(B8 + B18)/B10	0.21	<0.001	23.9	3.6	11.3	0.0
R2	Spring	(B4 + B17)/B5	0.66	<0.001	45.1	11.2	21.7	0.0
	Summer	B16/B6 − B16/B7	0.71	<0.001	29.6	5.2	15.3	0.0
R3	Spring	B7 * kdmin	0.15	<0.001	20.6	3.4	12.2	0.0
	Summer	iop_btot/conc_tsm	0.19	<0.001	21.5	3.0	11.7	0.0
R4	Spring	iop_apig/iop_agelb	0.13	<0.004	17.6	2.8	14.1	0.0
	Summer	B16/B6 − B16/B7	0.22	<0.001	19.9	3.4	12.4	0.0
R5	Spring	B17/B11 − B17/B12	0.18	0.05	21.8	6.7	17.9	0.0
	Summer	kd_z90max/B5	0.51	<0.001	19.6	3.3	15.5	0.0
R6	Spring	(B2 − B10) * B9	0.12	0.006	26.8	6.7	18.1	0.0
	Summer	conc_chl/iop_btot	0.17	0.003	25.0	3.7	13.7	0.0
R1–R6 combined	Spring	-	0.75	<0.001	24.3	5.6	15.2	0.0
	Summer	-	0.62	<0.001	23.3	3.7	12.9	0.0

with the formulas used. All the results are significant, except in R5 spring, where the p-value is 0.05. In addition, the coefficients of the quadratic polynomial regression model (a, b and c with the 95% confidence interval values (95%CI)) used in deriving TN are shown in

Table 8. The coefficients of the quadratic polynomial regression (a, b and c, used in Equation (1)) to derive TN with the 95% confidence interval values (95%CI).

Region	Season	a ± 95%CI	b ± 95%CI	c ± 95%CI
No regions	Spring	−7404652 ± 1171571	37505 ± 4133	18.7 ± 0.9
	Summer	40004 ± 1810	114.4 ± 110.6	21.3 ± 0.86
R1	Spring	1.9E+14 ± 9.1E+13	−4.4E+07 ± 16436166	20.1 ± 1.3
	Summer	134.6 ± 142.8	−259.8 ± 314.4	146.8 ± 173.3
R2	Spring	843.3 ± 490.1	−1452.6 ± 926.8	648.6 ± 431.9
	Summer	6723.3 ± 2037.9	504.6 ± 106.2	33.4 ± 2.6
R3	Spring	−8217.4 ± 89529.9	682.0 ± 811.6	19.1 ± 1.2
	Summer	15.1 ± 29.7	−8.0 ± 38.9	22.2 ± 9.9
R4	Spring	0.4 ± 0.3	−2.7 ± 2.2	20.7 ± 3.0
	Summer	−28849.5 ± 18696.6	−3911.1 ± 2652.5	−111.7 ± 93.6
R5	Spring	−101724 ± 125315	−65036 ± 80636	−10368 ± 12969
	Summer	0.0002 ± 0.0001	−0.07 ± 0.04	25.4 ± 4.7
R6	Spring	−9.2E+08 ± 667741628	−146626 ± 105740	26.5 ± 1.8
	Summer	4.5 ± 3.9	−10.2 ± 7.3	29.2 ± 3.0

. Only regions R2 and R4 in the summer periods have the same band ratio formula with the exact same bands. Almost half of the used formulas have the L2 product in them. In regions R3 and R6, the total particle scattering and either TSS or Chl-a concentrations gave the best results in summer period. In region R4, the best way to evaluate TN in spring is through the phytoplankton pigments and CDOM. The transparency products (K_d) are used in regions R3 and R5. Regions R1 and R2 algorithms used only reflectances in both periods. Only the reflectance bands at 400 and 442.5 nm were not used (Table 7). From the Table 6 and Table 7 (and Figure 4), it is evident that deriving TN is challenging in most of the regions (R1, R3, R4 and R6), and more in the spring period than in summer (MAPE is always higher in spring).

3.3. Total Phosphorus

There were more TP measurements than TN; the TP database contained 741 TP match-ups. The highest R² for each region (R1–R6) and each temporal division is shown in

Table 9. The number of match-ups (n) and the determination coefficients (R²) of the quadratic polynomial regression model for the best algorithm out of 25,013 for different regions (R1–R6) or region combinations for deriving total phosphorus (TP).

	April–Nov		May–Sept		April–May		June–Sept		April–June		July–Sept
	n	R2	n	R2	n	R2	n	R2	n	R2	n	R2
No regions	741	0.11	629	0.06	267	0.16	463	0.09	423	0.14	307	0.15
R1	143	0.09	119	0.04	42	0.26	97	0.06	73	0.17	66	0.20
R2	97	0.30	87	0.37	43	0.26	54	0.47	55	0.24	42	0.70
R3	209	0.18	175	0.05	82	0.25	126	0.09	138	0.23	70	0.33
R4	121	0.26	99	0.17	43	0.46	74	0.18	64	0.42	53	0.32
R5	47	0.34	43	0.45	9	0.78	38	0.43	22	0.34	25	0.52
R6	124	0.08	106	0.15	48	0.27	74	0.21	71	0.20	51	0.41
R1 + R3	352	0.11	294	0.04	124	0.21	223	0.05	211	0.20	136	0.12
R2 + R6	221	0.10	193	0.20	91	0.15	128	0.27	126	0.09	93	0.43
R4 + R5	168	0.28	142	0.19	52	0.46	112	0.20	86	0.45	78	0.25

. Additionally, similar to TN, combinations of different regions and no regions were used for algorithm testing. Generally, the relationships between derived algorithms and TP were weaker than with TN. The highest R² was usually found for the shortest period, April to May; however, at the same time, the July to September period gave very good results compared to other time divisions. For the same reasons as for TN, April to June + July to September periods were chosen for this study.

On the contrary to TN, the algorithms were performing very poorly in all temporal divisions when no regions were differenced (R² ranged 0.06–0.16). Using regions improves the accuracy of the estimation of TP for all the temporal divisions (Table 8). Similarly with TN estimations, algorithms performed better in regions R2 and R5; in the R1, R3 and R6 regions, the performance was the weakest along all the temporal divisions. The combination of geographically-close regions (R1 + R3, R2 + R6, and R4 + R5) did not improve the estimation of TP, and therefore was excluded from the further study (Table 9). The accuracy of the TP estimation improves from 69.9% to 76.2% (R² from 0.15 to 0.38) for the entire vegetation period (April to September) when regions are applied to the dataset (Figure 6).

The accuracy of estimating TP concentrations had significant differences between the different regions, being highest in region R2, at 84.7%, and lowest in region R1, at 67.6% (Figure 7). In all the cases, the remote sensing algorithms were underestimating TP concentrations. Still, the detection of TP was better in regions R2, R4, and R5 (R² was 0.60, 0.42 and 0.50, accordingly). The RMSE was between 0.11 and 0.36 µmolP L⁻¹.

Table 10. The best formulas (B in formulas indicates the C2RCC reflectance band, and the number indicates the band number), the determination coefficient (R²), p-value, mean and the root mean square error (RMSE) (both in µmolP L⁻¹), the mean absolute percentage error (MAPE), and bias (µmolP L⁻¹) for each region in spring (April to June) and summer (July to September) for the derived total phosphorus (TP).

Region	Season	Formula	R2	p-Value	Mean	RMSE	MAPE	Bias
No regions	Spring	(B2 + B5)/B4	0.14	<0.001	0.70	0.31	33.5	0.0
	Summer	conc_tsm/B4	0.15	<0.001	0.63	0.20	25.4	0.0
R1	Spring	(B8 − B10) * B1	0.17	<0.001	0.88	0.44	36.9	0.0
	Summer	kd489/iop_agelb	0.20	<0.001	0.75	0.23	27.5	0.0
R2	Spring	(B9/B3) * (B9/B3)	0.24	<0.001	0.57	0.10	15.1	0.0
	Summer	B17 * B18/B3	0.70	<0.001	0.66	0.12	15.5	0.0
R3	Spring	(B10/B9) * (B10/B9)	0.23	<0.001	0.78	0.28	25.9	0.0
	Summer	B5/(B4 + B10)	0.33	<0.001	0.64	0.15	18.9	0.0
R4	Spring	iop_atot/iop_agelb	0.42	<0.001	0.60	0.15	22.8	0.0
	Summer	B16/B4 − B16/B9	0.32	<0.001	0.53	0.11	15.6	0.0
R5	Spring	iop_bpart * conc_chl	0.34	0.004	0.38	0.13	43.8	0.0
	Summer	iop_apig * iop_adg	0.52	<0.001	0.49	0.11	20.6	0.0
R6	Spring	(B16/B5) * (B16/B5)	0.20	<0.001	0.63	0.19	23.2	0.0
	Summer	iop_adg * kdmin	0.41	<0.001	0.60	0.14	20.6	0.0
R1–R6 combined	Spring	-	0.34	<0.001	0.70	0.27	26.4	0.0
	Summer	-	0.46	<0.001	0.63	0.16	20.1	0.0

and

Table 11. The coefficients of the quadratic polynomial regression (a, b and c, used in Equation (1)) to derive TP with the 95% confidence interval values (95%CI).

Region	Season	a ± 95%CI	b ± 95%CI	c ± 95%CI
No regions	Spring	8.7 ± 5.4	−28.2 ± 18.5	23.4 ± 15.8
	Summer	3.8E−10 ± 3.4E−09	4.5E−05 ± 2.9E−05	0.5 ± 0.04
R1	Spring	4.8E+11 ± 4.6E+11	−787,522 ± 429,007	0.7 ± 0.2
	Summer	0.004 ± 0.002	−0.08 ± 0.04	1.1 ± 0.2
R2	Spring	−0.0005 ± 0.0002	0.02 ± 0.009	0.5 ± 0.06
	Summer	−41,430 ± 46,114	342.2 ± 141.5	0.5 ± 0.06
R3	Spring	44.9 ± 32.1	−90.3 ± 67.2	46.0 ± 35.1
	Summer	23.7 ± 8.2	−35.9 ± 12.5	14.1 ± 4.7
R4	Spring	0.001 ± 0.004	0.02 ± 0.06	0.4 ± 0.2
	Summer	−373.6 ± 210.5	−92.4 ± 55.3	−5.1 ± 3.6
R5	Spring	−0.0009 ± 0.001	0.05 ± 0.05	0.2 ± 0.2
	Summer	210.2 ± 104.9	−17.0 ± 10.3	0.7 ± 0.2
R6	Spring	16.3 ± 7.9	−3.5 ± 1.8	0.7 ± 0.08
	Summer	−0.15 ± 0.06	0.6 ± 0.2	0.5 ± 0.06

show the best formulas and the coefficients used in TP deriving algorithms. Similarly with TN algorithms, almost half of the formulas have L2 products in them. In region R5, Chl-a and scattering of the marine particles seems to be influencing TP in spring, compared to CDOM together with TSS and pigments in summer. The best algorithms for regions R1 and R6 have relationships with transparency L2 products (K_d). The algorithms for regions R2 and R3 used only reflectance bands. The reflectance bands on 560–620 nm and 709–754 nm were not used in any of the best algorithms. In deriving TP, the spring season had lower R² and higher MAPE (with some exceptions (R2 and R4)) compared with the summer season (Table 10).

4. Discussion

The detection of nutrients from optical remote sensing data still remains a challenge because TN and TP have no spectral response in the visible and near-infrared regions [25]. Several studies have estimated the nitrogen and phosphorus based on their strong relationships with Chl-a, TSS, and other optically sensitive parameters in the water [5,26,34,35,47,48,49]. In addition, a previous study has shown that Chl-a and SD are able to explain 41% of the variance in TP for the Swedish rivers discharging into the Baltic Sea [30]. Unfortunately, in the current study, no strong relationships were observed between the in situ Chl-a or SD with the nutrients in the match-up dataset (Table 5). TN showed only a moderate relationship with Chl-a in the R1 region. TP only had a significant relationship in the R4 region during spring (Figure 2A). The overall (all dataset) relationship between TN and SD is stronger than the relationship with Chl-a (R² = 0.27 and 0.05, accordingly), while no significant relationships between TN and SD were in any of the regions separately (R² = 0.01–0.15) (Table 5).

The analysis of the relationships of the nutrients with the other two optically active substances (a_CDOM and TSS) are limited because of the small dataset of in situ a_CDOM and TSS. Nevertheless, the match-up in situ database showed strong dependence of TN with a_CDOM, but only a weak relationship between TP and a_CDOM (Figure 3). The region R2 plays an important role in this strong relationship between a_CDOM and TN because of the optically dark river Pärnu, which is bringing freshwater rich with dissolved organic matter into the Pärnu Bay (region R2). This strong relationship makes the detection of TN from remote sensing data possible in the Pärnu Bay area. In the areas (or regions) where a_CDOM was lower, the results of deriving TN were weaker. The weak relationship between in situ TP and a_CDOM increases significantly (R² = 0.48, n = 13, p-value = 0.008) when R6 data is removed, but the dataset is too limited to draw any significant conclusions. In marine systems generally, N has been identified as the growth limiting nutrient, whereas in estuaries, P may be limiting in the freshwater part and N in the marine part [50]. The optimal TN:TP ratio for phytoplankton growth is 16:1; this is the so-called Redfield ratio [51]. However, except in deep-oceans, it is more of an exception than a rule. A TN:TP ratio > 50 indicates severe P limitation in the environment, a ratio < 20 indicates N limitation, and anything in between could indicate either of the nutrient deficiencies [52]. In the Estonian marine waters, the TN:TP ratio is very rarely <20, but at the same time, it is also not often >50. Only regions R2 and R5 show P limitation with an average TN:TP ratio of 65 and 58.5, respectively (Table 4). Table 4 indicates that the TN:TP ratio is very high in R2 and R5 during the spring period; they are also the regions where nutrients have been estimated rather well. In Guildford and Hecky [52], the TP was controlling the phytoplankton growth only when its concentration was less than 0.5 µmolP L⁻¹, regardless of the concentration of the TN. Based on this, we can say that only the region R5 is P limited. This is an addition to confirming the P deficiency in the R5 region, besides the TN:TP ratio being > 50 in the spring period. As seen in Table 4, it is shown that the TN decreases and the TP increases from spring to summer in regions R2 and R5, while the other regions show the opposite trend (except a small decrease of TN in R6). This kind of TN:TP fluctuation is another sign of the P limitation in the R2 and R5 regions, where TN can be depleted fast during the growing season and P can be resuspended from the sediments. High TP loads can lead to N-fixing cyanobacteria blooms in the summer [53]. Well-derived TP in the region R4 might be related to the SD changes in the water. Although in the region R4, the nutrients and Chl-a were quite stable throughout the seasons, the variability in the SD were the largest of all the regions (Table 4). The low Chl-a values in the presence of high TN:TP ratios are a sign of other factors that are limiting the growth of phytoplankton (other minerals, light, etc.) [54].

Low TN:TP ratios may cause N-fixing cyanobacterial blooms. In the regions R1, R3 and R4, N-fixing cyanobacteria can occur in summer [55,56], when the main driver for the phytoplankton growth might be temperature, wind speed, or an upwelling effect [57]. Therefore, TN and TP detection is even more challenging when the optically active parameters, like Chl-a, TSS, CDOM, have no major effect or relationship with nutrients.

The overall detection of TN was successful (R² = 0.73, MAPE = 14.2%, RMSE = 4.87 µmolN L⁻¹, n = 708) in the case that TN was first derived for each region separately, and then all the results were combined (Figure 4: R1–R6 combined). The overall success relies in a large part to the good results in the region R2 (Pärnu Bay) (R² = 0.74, MAPE = 18.6%, RMSE = 8.87 µmolN L⁻¹, n = 87) (Figure 5). As discussed above, the good results of the TN estimations in Pärnu Bay were because of the strong dependence with a_CDOM in the region.

The means of the in situ and derived TN show slight difference only in theR1 region in spring and the R6 region in summer (3% and 0.4%, accordingly) (Table 4 and Table 7). This suggests that remote sensing can be used to make the estimations on the spatio-temporal means of TN, which are used for the determination of the ecological states of the coastal waters.

The overall detection of TP was less successful compared to the TN results (R² = 0.38, MAPE = 23.8%, RMSE = 0.23 µmolP L⁻¹, n = 730) (Figure 6: R1–R6 combined). However, the TP detection was much more successful in four regions compared to TN: region R3 (R² = 0.28), R4 (R² = 0.42), R5 (R² = 0.50), and R6 (R² = 0.27), and only in R1 (R² = 0.20) and R2 (R² = 0.60) were the results were somewhat weaker (Figure 5 and Figure 7). In addition, the means of in situ and derived TP in different regions and seasons were the same, except a minor difference of 1.7% in the R4 region spring season (Table 4 and Table 10). This also shows the suitability of remote sensing for the TP estimations, especially where spatiotemporal means are monitored.

The spatial variation of TN and TP can have discontinuity at the borders of regions with the method used in the current study. Therefore, the algorithms of neighbouring regions should be blended at the borderline. Algorithm blending has shown good results in chlorophyll retrieval studies by Moore et al. [58,59]. The overall work scheme would include a weight determined by the distance from the border of the region for each pixel in the border area, which is used to obtain a final blended TN or TP concentration. This approach would be suitable for progressive transitions between regions, and blend separately tuned algorithms without suffering from the discontinuity associated with hard-classification schemes.

Past studies showed better results on TP than TN whenever the nutrients were derived with remote sensing data from freshwater systems [29,32]. Moreover, most of the inland nutrient studies derived only TP [28,30,31,33,34,36,37,38]. Only He et al. [18] and Huang et al. [35] obtained better results with TN than with TP in a freshwater reservoir. The studies performed in marine systems were more successful in detecting TN values [3], or derived only TN values [5,35,39]. Our study confirms that, in marine waters, detecting TN is more successful than deriving TP.

In coastal regions, which are influenced by large river run-off, nutrients are usually conservatively mixed with the sea water and have quantitative relationships with salinity [60]. Thus, Wang et al. [3] suggest that the use of salinity and spectral data may improve retrieval accuracy. Our database showed only relationships between TN and salinity in the spring period. In all the other periods, R² was less than 0.1. TP did not show any relationship with salinity in any period. As salinity cannot be directly detected in coastal waters with remote sensing and our database shows weak relationships with the nutrients, an optically active substance as a proxy is needed, and, at least in Baltic Sea, use of salinity might not improve results.

None of the past studies have had as large a dataset as in current study with the range of temporal and spatial variations (741 match-ups from 63 sampling stations over 36,500 km² study area over six years). Yu et al. [39] based their study only on a single remote sensing image. Wang et al. [3] had very high R² values (for TN 0.98–0.99, and for TP 0.75–0.86), but based their study mostly on in situ measured reflectances and an artificial database, and only 18 match-ups. In Tampa Bay (1000 km²), which is half of the size of the region R5 in the current study, the TN estimations were based on 103 match-ups over a three-year period, and had R² values of 0.75 (calibration dataset) and 0.63 (validation dataset) [5]. Only Chang et al. [40] had 740 match-ups from 52 cloud-free days in total for Tampa Bay, where inverse modelling resulted with R² 0.53 and 0.58 based on the calibration and validation datasets, respectively. Taking into consideration the large spatial and temporal variability of the data that the current study is based on, the results show high accuracy for the detection of nutrients with remote sensing data, and TN in particular is reliable for marine waters. Although not applicable for the entire year, remotely sensed data with high spatial and temporal resolution could be used in nutrient monitoring, ecosystem modelling, or estimation of the ecological state of the coastal and offshore waters during the vegetation period. These algorithms developed in the current study can be used as guidance for nutrient mapping in other coastal waters. The ranges of optically active substances concentrations and TN:TP ratios should be used as an indication of whether the same algorithms can be used in different coastal waters. However, coastal waters around the world are extremely variable. Therefore, it is highly likely that the best approach to find out whether TN and TP can be mapped with remote sensing in a particular site is simply to repeat this study with the whole set of algorithms.

5. Conclusions

Nutrients play a major role in the eutrophication. Therefore, monitoring of the nutrients loads needs to be improved. The estimation of nutrients like TN and TP is very challenging with optical remote sensing data, because TN and TP are not directly related to the water colour that the remote sensing instruments are detecting. Nevertheless, we showed that remote sensing is very useful in detecting those non-optically active substances in the water. In this study, TN was well derived with R² = 0.73, RMSE = 4.87 µmolN L⁻¹, MAPE = 14%, and n = 708. The TP estimations in Estonian marine waters were not as good, with R² = 0.38, RMSE = 0.23 µmolP L⁻¹, MAPE = 24%, and n = 730. The estimations of both nutrients were the most successful in the region R2. This region represents the semi-closed, small, and turbid Pärnu Bay, which is a P limited water body with very high CDOM and a nutrient rich river inlet. The deficiency of P seemed to have a positive influence on the accuracy of the nutrient detection, which was also shown by the rather good results derived in the region R5 (Väinameri Sea). Based on the results of this study, we may say that deriving TN from optical remote sensing data is feasible in all Estonian coastal waters throughout the whole ice-free season. Mapping of TP concentrations in Estonian marine waters using remote sensing data should be taken more cautiously. Only the Pärnu Bay, the Baltic Proper area of the west Estonian archipelago, and the Väinameri Sea (regions R2, R4, and R5, accordingly) are the regions where TP can be mapped with sufficient accuracy.