Jukuri, open repository of the Natural Resources Institute Finland (Luke) All material supplied via Jukuri is protected by copyright and other intellectual property rights. Duplication or sale, in electronic or print form, of any part of the repository collections is prohibited. Making electronic or print copies of the material is permitted only for your own personal use or for educational purposes. For other purposes, this article may be used in accordance with the publisher’s terms. There may be differences between this version and the publisher’s version. You are advised to cite the publisher’s version. This is an electronic reprint of the original article. This reprint may differ from the original in pagination and typographic detail. Author(s): Niko Kulha, Leena Ruha, Sakari Väkevä, Sampsa Koponen, Markku Viitasalo, Elina A. Virtanen Title: Satellite bathymetry estimation in the optically complex northern Baltic Sea Year: 2024 Version: Published version Copyright: The Author(s) 2024 Rights: CC BY 4.0 Rights url: https://creativecommons.org/licenses/by/4.0/ Please cite the original version: Kulha, N., Ruha, L., Väkevä, S., Koponen, S., Viitasalo, M., & Virtanen, E. A. (2024). Satellite bathymetry estimation in the optically complex northern Baltic Sea. Estuarine, Coastal and Shelf Science, 298, 108634. https://doi.org/10.1016/j.ecss.2024.108634. Estuarine, Coastal and Shelf Science 298 (2024) 108634 Available online 10 January 2024 0272-7714/© 2024 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Satellite bathymetry estimation in the optically complex northern Baltic Sea Niko Kulha a,*, Leena Ruha b,c, Sakari Vakeva a, Sampsa Koponen a, Markku Viitasalo a, Elina A. Virtanen a a Finnish Environment Institute, Latokartanonkaari 11, 00790, Helsinki, Finland b Research Unit of Mathematical Sciences, University of Oulu, P.O. Box 8000, 90014, Oulu, Finland c Natural Resources Institute Finland, Latokartanonkaari 9, 00790, Helsinki, Finland A R T I C L E I N F O Keywords: Satellite bathymetry estimation Optically complex waters Baltic Sea Sentinel-2 A B S T R A C T Satellite bathymetry estimation (SBE) is a cost-effective method for producing depth information for shallow marine and aquatic areas. However, most SBE methods have been developed and are refined in clear water conditions and their feasibility in optically complex waters is poorly understood. This hinders the use of SBE in optically complex waters. To test the feasibility of SBE in optically complex waters, we used multi-temporal Sentinel-2 (S-2) imagery and in-situ measured depth information from dive computers to fit bathymetric models for 14 study areas in the northern Baltic Sea. We compared different methods of using multi-temporal S2-scenes, atmospheric correction, and statistical modeling, as well as different band combinations to define the optimal model for each study area. We tested the dependence of the depth estimation on 13 environmental variables with potential influence on the accuracy of the SBE. The optimal bathymetric models differed among the study sites, with the differences manifested in different requirements for atmospheric correction, band combination, statistical model, and the use of multiple S2-scenes. The coefficients of determination (R2) of the optimal models varied between 0.42 and 0.79, with mean (X) R2 equaling to 0.65. The coefficients of determination were positively related to the proportion of sea within the 5 km search kernel and to the exposure of the study site. At the predicted depth of 2 m, the differences between observed and predicted depths varied between 0.11 and 0.60 m (X ˆ 0.35 m). At the predicted depth of 1 m, the differences varied between 0.04 and 0.34 m (X ˆ 0.21 m). SBE provides accurate depth estimates in optically complex waters, even when in-situ calibration data are limited. However, the maximum detectable depth obtained with SBE in optically complex waters is significantly lower than in clear water conditions. The SBE process can be fine-tuned to the specific characteristics of an area, improving its performance in optically complex conditions. In particular, the identification of the optimal band combination is critical for successful SBE. The depth estimates produced in this study can be used in a variety of applications that require accurate information on the seafloor bathymetry. 1. Introduction In the marine environment, accurate bathymetric information is needed for many applications, such as maritime safety, environmental monitoring, maritime spatial planning, and ecological research (Gao, 2009). However, mapping shallow water bathymetry is virtually impossible using traditional ship-based methods such as sonar or echo sounding (Tronvig, 2005) and expensive using airborne solutions such as bathymetric LiDAR (Guo et al., 2022). Shallow photic areas are often the most ecologically interesting, as they support the majority of submerged vegetation and a diverse biota (Kotta et al., 2018; Krause-- Jensen et al., 2009; Virtanen et al., 2018). This is particularly the case when light is rapidly attenuated in the water column (Luhtala et al., 2013), which physiologically limits the occurrence of submerged aquatic vegetation to near-surface depths (Lappalainen et al., 2019; Luhtala et al., 2016). Optical remote sensing imagery can be used to map bathymetry in shallow water areas (Chybicki, 2017; Geyman and Maloof, 2019). To date, observations from instruments such as Quickbird (Lyons et al., 2011), Worldview (Doxani et al., 2012), Landsat (Knudby et al., 2016), * Corresponding author. Current address: Natural Resources Institute Finland (Luke), Latokartanonkaari 9, 00790, Helsinki, Finland. E-mail address: niko.kulha@gmail.com (N. Kulha). Contents lists available at ScienceDirect Estuarine, Coastal and Shelf Science journal homepage: www.elsevier.com/locate/ecss https://doi.org/10.1016/j.ecss.2024.108634 Received 22 June 2023; Received in revised form 27 September 2023; Accepted 9 January 2024 Estuarine, Coastal and Shelf Science 298 (2024) 108634 2 and Sentinel (Caballero and Stumpf, 2023) have been successfully used for satellite bathymetry estimation (SBE). The advantage of Landsat and Sentinel data is their current free, open, and long-term availability, which enable cost-effective bathymetry estimation over large areas. However, compared to the Landsat data, the Sentinel instruments have shorter revisit times and higher spatial resolution in the visible bands, providing good premises for SBE (Hedley et al., 2018). Optical sensing of bathymetry is based on the notion that the spectral curve shape and the amount of radiant energy scattered by the water column and reflected from the bottom is a function of water depth. In empirical modeling, the relationship between sensed water radiance and in-situ measured depth values is established empirically (Gao, 2009). Among the empirical SBE methods that require in-situ data for calibra- tion, the linear band model (Lyzenga, 1978) and the log ratio model (Stumpf et al., 2003) are among the most widely used. While both ap- proaches estimate bathymetry based on the differences in the attenua- tion of different wavelengths of visible light in water, the linear band model is more sensitive to variations in bottom albedo than the log ratio model (Doxani et al., 2012; Geyman and Maloof, 2019). As a result, the log ratio model can outperform the linear band model under conditions of varying bottom albedo and provide accurate depth estimates even outside the range of in-situ data (Geyman and Maloof, 2019). The performance of the log ratio model depends on, e.g., the optical properties of the water column (and the processes that affect these properties), the bands used, and atmospheric the conditions (Gao, 2009). Several methods have been developed to address these chal- lenges. For example, the use of multiple satellite scenes (referred to as multi-scene compositing) can reduce the adverse effects of transient changes in the atmosphere (thin clouds or haze) or in the optical prop- erties of the water column, such as turbidity or algal blooms (Caballero and Stumpf, 2020a, 2023). Furthermore, the best combination of bands depends on the optical properties of the water column at a given location (Casal et al., 2019). For example, blue wavelengths are typically feasible in clear water, whereas longer wavelengths such as yellow and red may perform better in turbid conditions (Vahtmae and Kutser, 2016). Cor- recting the scenes used in SBE for atmospheric distortions may produce more accurate depth estimations especially in optically complex waters (Caballero and Stumpf, 2020b). Most SBE methods have been developed and fine-tuned in clear oceanic waters (Hedley et al., 2018; Lyzenga, 1978; Stumpf et al., 2003). Therefore, these methods may not be directly applicable in optically complex (Case 2) waters (Casal et al., 2019; Geyman and Maloof, 2019; Lafon et al., 2002), such as in the Baltic Sea (Kutser et al., 2018). Optically, the entire Baltic Sea can be considered a Case 2 sea, with high concentrations of colored dissolved organic matter (CDOM) (Darecki and Stramski, 2004) and suspended particulate matter, with the con- centrations peaking near the estuaries and during spring and fall, when river discharge is highest (Asmala et al., 2012; Kyryliuk and Kratzer, 2019; Neumann et al., 2021). Similarly, the of chlorophyll a concen- trations show spatiotemporal variability, with the highest concentra- tions generally occurring in May and July in the northern Baltic Sea (Attila et al., 2018; Suominen, 2018). Consequently, the optical prop- erties of the Baltic Sea, change throughout the ice-free season (Stock, 2015), even diurnally (Luhtala et al., 2013), and vary between its sub-basins (Kratzer et al., 2003; Stock, 2015). Due to the spatiotemporal variability of optical properties, oceanic remote sensing methods need to be calibrated to the site-specific conditions of the Baltic Sea (Darecki et al., 2003; Kratzer et al., 2003). In this study, we used Sentinel-2 (S2) imagery and in-situ depth ob- servations in a log ratio model to identify the optimal methods for generating SBEs in the optically complex waters by fitting models at 14 study sites distributed throughout the northern Baltic Sea. Specifically, our objectives were to determine (1) the optimal method for SBE is in the northern Baltic Sea, (2) how well the S2-based bathymetric models perform in these optically complex waters, and (3) what factors affect the model performance. 2. Materials and methods 2.1. Study area We selected study sites along the Finnish coast, where the shoreline types and profiles vary. Sharply sloping rocky shores are common in the south, while flat sandy shores are common in the north. We selected the sites with specific criteria (in order of importance): (1) the sites are distributed along the whole coastline, (2) a maximum of four sites are in the same sub-basin as defined by the Baltic Marine Environment Pro- tection Commission (HELCOM) (Fig. 1), (3) no major estuaries, ports, or other features that would cause variations in the optical properties of the area are located within or near the site, (4) the sites have been field- sampled, and (5) at least five cloud-free Sentinel-2 – scenes, captured in or after 2016 and depicting the site during clear water conditions, are available. Adhering to these criteria, we selected 14 study sites located throughout the Finnish coast (Figs. S1–S14). 2.2. Sentinel-2 data acquisition and processing For each study site, we downloaded five cloud-free Level-1C (L1C) S2-images taken during clear water conditions between 2016 and 2020 using the Sentinel Scientific Data Hub (see Table S1 for image details). Correcting the images for atmospheric distortions may produce more accurate depth estimations especially in optically complex waters (Ca- ballero and Stumpf, 2020b). Therefore, we used two atmospheric Fig. 1. The locations of the study areas. The red lines represent the HELCOM sub-basin boundaries. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) N. Kulha et al. Estuarine, Coastal and Shelf Science 298 (2024) 108634 3 correction procedures with demonstrated feasibility in SBE (Caballero and Stumpf, 2020b; Casal et al., 2019): ACOLITE Exponential model (EXP) and ACOLITE Dark Spectrum Fitting Model (DSF) alongside the L1C images. To test the performance of different band combinations in SBE, we produced depth estimates using three different band combinations: blue (central wavelength, CW, 490 nm)/red (CW ˆ 665 nm), blue/green (CW ˆ 560 nm), and green/red at all sites. We formed all possible combinations of the L1C, EXP and DSF-scenes and band combinations and implemented the log ratio model of Stumpf et al. (2003) as EDˆm1  log…Bm† log…Bn† m0 eq. 1 where Bm and Bn are water-leaving reflectance values at specific bands, and m1, m0 are constants scale and offset, respectively, that transform the model results to actual depth values. ED is estimated depth that is depicted as a raster map with 10 m spatial resolution. We denote these estimated depths ED-maps. Each study site had five L1C scenes and five scenes corrected with EXP and DSF procedures, respectively. For each site, we overlaid the ED-maps produced with scenes using same atmo- spheric correction and band combination (e.g., five L1C-maps produced with bands blue/red, five DSF-maps produced with bands blue/green). Henceforth we denote these ED-stacks. Each study area had nine different ED-stacks. We calculated a normalized difference water index (NDWI) and tried to use it to distinguish land areas from water. Due to the poor perfor- mance of NDWI near reed beds and in areas where other aquatic vege- tation grew above the surface, we also calculated a normalized difference vegetation index (NDVI) for each L1C-scene. Due to the su- perior performance of NDVI compared to NDWI, we used NDVI to mask land areas from each ED-stack by removing cells with NDVI 0.1 in any L1C-scene at a given site from all stacked ED maps for that site. 2.3. In-situ depth data used to calibrate the satellite bathymetry estimates In-situ-measurements are needed to calibrate the remotely sensed depth estimates, and they can also be used to quantify the error in the remotely sensed records (Kulha et al., 2018). We used diver-recorded depth information collected within the Finnish Inventory Programme for Underwater Marine Diversity (VELMU) to calibrate the S2-based depth estimates and to quantify the error in depth estimation. The VELMU programme collects information on species and habitats by scuba-diving and drop-video surveys (for details on the sampling see Virtanen et al., 2018). The dives are typically made in approximately 100 m long transects, perpendicular to shoreline, where sampling of species and substrates (% cover) is done in horizontal 10 m (along the seabed profile) or vertical 1 m (change in depth) intervals, from 4 m2 inspection squares. Transects end either at the deepest recorded vege- tation, or when the end of 100 m line is reached. The diver records the depth of each sampled site underwater, using a dive computer. We used the availability of dive transect data to define the study areas within the selected study sites. In sites where dive transects were scarce, a larger study area was selected to ensure sufficient in-situ observations for calibration of the satellite bathymetry estimates (Figs. S1–S14). We included in the analyses all sampling points that had a recorded depth3m within the study area. In total, we used 3433 in-situ depth observations, the number of observations per study site ranging from 96 in Pajukari to 584 in Trato (median number of observations per study site ˆ 279) (Table 1). 2.4. Optimal models for spatially explicit depth estimation In each study site, we identified the optimal method for SBE by building SBE-models using different elements. The elements included S2-scenes processed with different AC procedures, as well as different band combinations, statistical models, and multi-scene compositing methods. We overlaid the in-situ measurements on the ED-stacks and extracted ED-stack values to the dive points. Due to the nested structure of our data (samples within dive transects, multiple samples from a single S2-scene pixel), we used mixed-effect models with dive transect and S2-pixel as random effects. When the relationship between the in- situ measurements and the ED-stack values was linear, we implemented a linear mixed model (lmm) and, as our response variable is strictly positive, generalized linear mixed model with gamma-distributed penalized quasi-likelihood generalized mixed model (glmmPQL). If the relationship was nonlinear, we used a generalized additive mixed model (gamm). After accounting for the nested structure in our data, the model residuals were not spatially autocorrelated. The water level fluctuations evident between the images taken at different times may cause random error in the bathymetry estimation. However, since the Baltic Sea is non- tidal, it can be assumed that the water level differences are small be- tween the images and the influence of water level fluctuation is negli- gible. For the same reason, tidal corrections were not performed. Although the use of multiple satellite scenes (i.e., multi-scene compositing) may enhance SBE, the optimal compositing method is not evident. Hence, we first performed multi-scene compositing using all five scenes in one ED-stack. We denote these compositions as All. Then, in order to reduce the number of explanatory variables, we performed the multi-scene compositing with three additional methods, using (1) the site-specific mean of the five scenes in one ED-stack, (2) a prediction from a linear model where the values of one scene from an ED-stack were predicted using values from the other four scenes of the stack, and (3) a site-specific maximum pixel values among the five scenes in a particular ED-stack (Caballero and Stumpf, 2020b). We denote these compositions as Mean, Linear, and Maximum, respectively. We assessed the performance of each SBE-model version to identify the best combination of AC procedure, band combination, statistical model, and multi-scene compositing method individually in each study site. We examined model performances by comparing the coefficients of determination (R2) and root mean squared errors (RMSE) from each calibration model version using a leave-one-subject-out cross-validation (Hastie et al., 2009). Here, we dropped one dive transect at a time from the dataset and used the remaining data to fit the models. Then, we used the re-fitted models to predict depth values for the sampling points within the removed transect. We repeated the procedure for each dive transect and every calibration model version at each studied site to identify the model with the highest R2 (smallest RMSE) for each studied area. We denote these optimal models. We fitted the models in R 4.2.2 (R Core Team, 2019), using the package lme4 for lmm-models (Bates et al., 2015), mgcv for gamm-models (Wood, 2017), and MASS to for the Table 1 Description of the in-situ depth observations used to calibrate the satellite ba- thymetry estimates. Study area Study area size (km2) Number of in-situ depth observations Number of dive transects Average number of depth observations per transect Average in- situ measured depth (m) Bergon 889 373 55 7 1.16 Hamskeri 36 278 34 8 1.11 Kallahti 6 203 23 9 1.31 Kampholmen 200 478 51 9 0.91 Lepanen 53 105 16 7 1.26 Molpe 240 245 29 9 0.82 Mussalo 139 280 25 11 1.01 Pajukari 254 96 12 9 1.72 Putsaari 140 325 44 7 1.26 Storskar 1069 361 52 7 1.38 Trato 1296 584 82 7 1.16 Ulkokrunni 106 164 18 9 1.08 Vimpasaari 32 138 17 8 0.74 Åva 817 289 120 9 1.32 N. Kulha et al. Estuarine, Coastal and Shelf Science 298 (2024) 108634 4 glmmPQL-models (Venables and Ripley, 2002). 2.5. Spatial variation in the quality of satellite bathymetry estimation To visualize and compare the depth estimation quality and accuracy at different study sites, we used the optimal models to predict depth values and examined the difference between predicted and measured depths in points where depths of 2 m or 1 m were predicted. Moreover, we used visual examination and a linear model to test whether the co- efficients of determinations of the optimal models depended on the (1) size of the study area, (2) mean in-situ-measured depth in the study area, (3) average Secchi-depth in the study areas, (4) wave exposure of the study area, defined as the potential for wave action at a particular location and quantified following Isæus (2004), (5) the proportion of the study area that is sea, (6) proportion of sand, (7) bedrock, (8) hard reflective bottom substrate, (9) silt, (10) clay, (11) mud, (12) other soft bottom substrate on the study area, and (13) number of different bottom substrate types in the sites where in-situ-data were measured. We calculated the size of the study area using a convex hull drawn around the in-situ-sampling sites in the area. The depth, Secchi, and bottom substrate information were in-situ-observations obtained from the VELMU database, and the proportion of sea from Virtanen et al. (2018). 3. Results 3.1. Optimal models for satellite bathymetry estimation in the northern Baltic Sea The optimal models for satellite bathymetry estimation differed be- tween the study sites (Fig. 2, Table S2). In nine study areas, the models with highest R2 (smallest RMSE) had the scenes with no atmospheric correction as predictors (i.e. L1C-scenes; Fig. 2A). Scenes corrected with the ACOLITE Dark Spectrum Fitting Model (DSF) produced a model with the highest R2 (smallest RMSE) in three study areas, and scenes cor- rected with the ACOLITE Exponential model (EXP) in two study areas (Fig. 2A). The logarithmic ratio of blue to red returned the highest R2 in nine study areas, followed by blue/green (three study areas) and green/red (two study areas) (Fig. 2B). Gamm performed best in nine study areas, lmm in four and glmmPQL in one study area (Fig. 2C). In six study areas the optimal model had the mean of the five scenes as a predictor, while in five study areas the optimal model had all five S2-scenes as predictors (Fig. 2D). In two study areas a linear model between the scenes produced the highest R2, and in one study area a composite of maximum pixel values produced the highest R2 (Fig. 2D). 3.2. The accuracy of the optimal models in depth estimation The accuracy of the depth estimation, measured as the difference between predicted and measured depth at the predicted depths of 2 m and 1 m, varied between the study areas (Fig. 3A and B). At the measured depth of 2 m, the difference ranged from 0.11 m at Hamskeri in the Bothnian Sea to 0.60 m at Bergon in the Åland Sea (Fig. 3A). At the measured depth of 1 m, the difference varied between 0.09 m at Storskar in the Quark and 0.34 m at Bergon (Fig. 3B). In each study site except of Hamskeri, Mussalo (Åland Sea), and Trato (Åland Sea), the difference was smaller at the depth of 1 m than in the at the depth of 2 m. Similarly, the R2 and RMSE varied between the study areas, the R2 ranging from 0.42 to 0.79 and the RMSE from 0.07 m to 0.38 m (Fig. 3C and D). Based on visual estimation, the optimal models produced reasonable bathy- metric maps of the study areas (Fig. 4). The prediction accuracy of the optimal bathymetric models varied between the study regions (Fig. 5, Table S2). In certain regions, such as at Åva (Åland Sea) (Fig. 5A), Lepanen (Bothnian Bay) (Fig. 5F) and Putsaari (Bothnian Sea) (Fig. 5J), the underestimation in the predicted depth values increased when the measured depth values were larger than 2 m. We examined the coefficients of determinations (R2) of the optimal models to investigate what variables are associated with the model fits in different study areas. The R2s were positively related to the propor- tion of sea in the study area and to the exposure of the study area at the significance level of 0.05 (Fig. 6). A linear model that had the co- efficients of determination as a dependent variable and proportion of sea as a predictor had an R2 of 0.26 (RMSE  0.10, p ˆ 0.03), and a linear model that had the coefficients of determination as a dependent variable and the exposure of the study area as a predictor had an R2 of 0.47 (RMSE  0.09, p < 0.01). The R2s of the optimal models varied inde- pendent of study area size, mean in-situ measured depth, Secchi-depth, proportion of sand, bedrock, hard reflective bottom substrate, silt, clay, mud, soft bottom substrate, and the number of bottom substrate types in the field-sampling sites (Fig. S15). 4. Discussion The parameters of the optimal models for satellite bathymetry esti- mation (SBE) differed between the 14 study areas, even within the same Fig. 2. The atmospheric correction (A), band combination (B), statistical model (C), and multi-scene compositing method (D) used in the identified optimal models with which the satellite bathymetry estimates were produced. N. Kulha et al. Estuarine, Coastal and Shelf Science 298 (2024) 108634 5 sub-basin. The optimal parameters for SBE-models depend on the characteristics of the area where the depth estimation is performed (Casal et al., 2019; Gao, 2009). The fact that we did not identify a single optimal model for SBE in the northern Baltic Sea illustrates the vari- ability in optical properties within the northern Baltic Sea as well as spatial variability in atmospheric conditions that may influence the quality of SBE (Caballero and Stumpf, 2020b). The optical properties of the Baltic Sea may vary both between (Kratzer et al., 2003; Stock, 2015) and within (Luhtala et al., 2013; Suominen, 2018) its sub-basins, due to – for example – differences in the concentration of CDOM (Darecki and Stramski, 2004), chlorophyll a (Pitarch et al., 2016) and suspended particulate matter (Fleming-Lehtinen and Laamanen, 2012). These fac- tors vary due to, for example, differences in river runoff (Asmala et al., 2012) and nutrient availability (Vahtera et al., 2007), and as a function of distance to coast (Luhtala et al., 2013). The fact that we did not identify a single optimal SBE models suggests that various approaches should be tested when conducting SBE in optically complex waters. The underestimation in the predicted depth values started to increase when the measured depth values exceeded 2 m. This suggests that the reflectance signal saturates rapidly in the water column and indicates that the maximum depth from which reliable SBEs may be obtained is considerably less in the optically complex northern Baltic Sea than in clear water conditions where depths beyond ten (Caballero and Stumpf, 2019) or even twenty (Caballero and Stumpf, 2023; Hedley et al., 2018) meters or more have been accurately estimated with optical remote sensing. Due to the rapid saturation of the reflectance signal, the us- ability of SBE is particularly high on optically complex coasts with a gentle depth profile, such as the sandy coasts of the northern Baltic Sea. However, in areas where the depth profile is steep, such as the rocky shores common in the southern parts of the Finnish maritime area, less depth estimates may be obtained. Waters with different turbidity levels scatter incoming radiation differently (Lafon et al., 2002), highly turbid conditions lowering the detectable water depth as light rapidly attenuates in the water column Fig. 3. The difference between predicted and measured depth at the predicted depth of 2 m (A) and 1 m (B), the coefficients of determination (R2) (C), and root mean squared errors (D) of the optimal models. Fig. 4. Part of the bathymetric map generated with SBE for the Ulkokrunni study area (approximately Nº65.38, Eº24.83). Panel A shows the location of the area on a Sentinel-2 image. Panel B shows the estimated depths and dive transects, indicated by yellow circles, on a false-color aerial image derived from the National Land Survey of Finland. Panel C shows the false-color aerial image without the estimated depths, with light blue colors indicating shallow water areas. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) N. Kulha et al. Estuarine, Coastal and Shelf Science 298 (2024) 108634 6 before reaching the seafloor (Gao, 2009; Luhtala et al., 2013). Accord- ingly, in optically complex waters, such in the Dublin Bay, Ireland, the SBE is limited to <10 m depths, the best fit (highest R2) between the SBE-estimates and in-situ depth observations occurring between the depths of 2–6 m (Casal et al., 2019). In Singapore Strait, an approach similar to ours yielded depth estimates between the depths of 0–4 m (Bramante et al., 2013) that are comparable to our results. In the coast of Estonia, Baltic Sea, imagery from a space-borne remote sensor Hyperion has successfully been used for depth estimation down to 3 m depth (Vahtmae and Kutser, 2007), while airborne hyperspectral data pro- duced reliable depth estimates in depths shallower than 4 m (Vahtmae and Kutser, 2016). The model using the logarithmic ratio of blue to red (b/r) often had the highest R2 (smallest RMSE), followed by the models using the ratios of blue to green (b/g) and green to red (g/r), respectively. In each study area, the five best performing models typically used the same band ratio, indicating that identifying the correct band ratio is critical to successful SBE in optically complex waters. Due to the high penetration in the blue spectral range and low penetration in the red spectral range into the water column (Gao, 2009), the b/r-model has shown good feasibility especially in SBE of shallow waters with moderate turbidity (Caballero et al., 2019; Caballero and Stumpf, 2019). In turbid conditions, the spectral peak should shift to wavelengths longer than blue, such as green (Casal et al., 2019). This has also been observed in the turbid waters of the Baltic Sea (Vahtmae and Kutser, 2007). However, using airborne hyperspectral images, Vahtmae and Kutser (2016) showed the best correlation between the sensed radiance of water and in-situ measured depth using the band ratio of 515/606 nm that is close to the 490/560 nm available in Sentinel-2, respectively. Although the coastal aerosol band (B01) in Sentinel-2 may improve SBE in turbid conditions (Ca- ballero and Stumpf, 2019), its coarse spatial resolution (60 m), potential sensitivity for organic matter in water column and sensitivity to atmo- spherics effects difficult its use in the turbid and spatially heterogeneous coast of the northern Baltic Sea (Vahtmae and Kutser, 2007). Atmospheric correction enhanced SBE in only five study areas. The purpose of atmospheric correction procedures is to reduce the influence Fig. 5. The relationship between the predicted and observed depths in the study areas. The density distributions show the relative number of observations in different depth classes. Fig. 6. The optimal model coefficients of determination as a function of the proportion of the study area that is sea within the 5 km search kernel (A) and the wave exposure, a unitless measure of potential for wave activity, of the study area (B). The shaded gray areas show the 95% confidence interval of the predictions. N. Kulha et al. Estuarine, Coastal and Shelf Science 298 (2024) 108634 7 of atmospheric distortions and scattering in retrieving water-leaving reflectance. Despite its shown feasibility in SBE of optically complex waters (Caballero and Stumpf, 2020b; Casal et al., 2019), atmospheric corrections may also increase the error noise in derived radiance in dark and optically complex waters (Werdell et al., 2010), such as in the northern Baltic Sea. Hence, the introduced noise may only increase the relative error in SBE, decreasing the accuracy of SBE. Due to complex structure of the coast and variable depth profiles the atmospherically corrected scenes could not be processed with spatial filtering to reduce noise (Caballero and Stumpf, 2020b), as this would have excessively removed information from SBE. The depth estimation did not always result in linear relationship between the in-situ measured and satellite-derived depth values. Consequently, a bathymetric model that used gamm produced the highest R2 (smallest RMSE) in nine out of 14 study areas. The lack of linear relationship was plausibly due to the different performance of empirical modeling methods for SBE in optically complex environments compared to clear water conditions (Casal et al., 2019). Due to the op- tical complexity, oceanic remote sensing methods typically need to be calibrated to the site-specific conditions in the Baltic Sea (Darecki et al., 2003; Kratzer et al., 2003). As the mixed models accounted for the nested structure and removed spatial residual autocorrelation in the in-situ data, we did not need to use additional method to remove spatial autocorrelation, such as kriging with external drift (Casal et al., 2020). The coefficients of determination of the optimal models varied within the study areas, depending on the used multi-scene compositing method. Multi-scene compositing has been shown to reduce the adverse impact of momentary variation in the optical properties of water column due to, e.g., turbidity or algal blooms (Caballero and Stumpf, 2020a; Sagawa et al., 2019). However, our observation suggests that the method of multi-scene compositing is not trivial and multiple methods to use the multi-temporal scenes should be explored. Further, due to complex structure of the coast in our study areas and the consequent spatial variability in the optical properties of water column, we further argue that multi-scene compositing outperforms the use of a single optimal scene (Evagorou et al., 2019). Of the variables tested, only the proportion of the study area that is sea within 5 km search kernel (hereafter sea share) and the mean exposure of the study area showed a relationship between the co- efficients of determination of the optimal models, indicating that the model fits varied independent of, e.g., bottom albedo. In the Baltic Sea, the concentration of CDOM is high especially close to estuaries (Asmala et al., 2012) and the concentration of other suspended matter decreases as a function of distance to coast (Luhtala et al., 2013), influencing the optical properties of water column. Hence, the observed positive re- lationships between increasing sea share and exposure and increasing coefficients of determination probably reflect differences in water quality, the areas with high sea share and exposure having more clear waters than their counterparts. In the northern Baltic Sea, the depth zone from which we were able to obtain reliable SBEs represents an ecologically important zone, dominated by submerged macrophyte communities (Lappalainen et al., 2019; Luhtala et al., 2016). These shallow areas host a variety of species with narrow distributional ranges, and are of high ecological value in the northern Baltic Sea (Virtanen et al., 2018). Although our SBE was limited to shallow water areas, we argue that improving depth infor- mation in these particular areas is extremely valuable, e.g. for the con- struction of species distribution models, commonly used for instance in conservation planning. The accuracy of such models depends on the quality of environmental predictors (Araújo et al., 2019), including bathymetry in marine applications. Until now, the species distribution models used in the Finnish seas have been based on depth data used in coastal navigation charts. These data are rather inaccurate in the shallow areas where no echo soundings have been made. Our new depth estimates provide data with the potential to significantly increase the accuracy of species niche predictions. However, the spatial resolution of Sentinel-2, 10 m, makes it difficult to estimate depth at the shoreline. In future, SBE would benefit of the use of imagery from instruments with fine spatial resolution. 5. Conclusions Our comprehensive analyses show that optical sensing of bathymetry provides valuable depth information in the shallow coastal regions of the northern Baltic Sea. However, the maximum detectable depth ob- tained with optical remote sensing of bathymetry is significantly lower in these optically complex waters compared to clear water conditions. In optically complex waters, the optimal procedure for producing remotely sensed depth estimates should be adapted to the specific characteristics of the estimated site. In particular, identifying the optimal band com- bination is critical for successful satellite bathymetry estimation. Our study provides an example of how to define an optimal method for satellite bathymetry estimation in optically complex waters. Funding information This research was funded by the Finnish Ministry of the Environ- ment’s VELMU project and the LIFE-IP BIODIVERSEA project LIFE20 IPE/FI/000020. CRediT authorship contribution statement Niko Kulha: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. Leena Ruha: Writing – re- view & editing, Validation, Supervision, Methodology, Formal analysis. Sakari Vakeva: Writing – review & editing, Validation, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. Sampsa Koponen: Writing – review & editing, Validation, Supervision, Methodology, Conceptualization. Markku Viitasalo: Writing – review & editing, Validation, Supervision, Project administration, Funding acquisition, Conceptualization. Elina A. Virtanen: Writing – review & editing, Visualization, Validation, Supervision, Resources, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability The data that has been used is confidential. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.ecss.2024.108634. References Araújo, M.B., Anderson, R.P., Marcia Barbosa, A., Beale, C.M., Dormann, C.F., Early, R., Garcia, R.A., Guisan, A., Maiorano, L., Naimi, B., O’Hara, R.B., Zimmermann, N.E., Rahbek, C., 2019. Standards for distribution models in biodiversity assessments. Sci. Adv. 5, eaat4858 https://doi.org/10.1126/sciadv.aat4858. Asmala, E., Stedmon, C.A., Thomas, D.N., 2012. Linking CDOM spectral absorption to dissolved organic carbon concentrations and loadings in boreal estuaries. Estuar. Coast Shelf Sci. 111, 107–117. https://doi.org/10.1016/j.ecss.2012.06.015. Attila, J., Kauppila, P., Kallio, K.Y., Alasalmi, H., Keto, V., Bruun, E., Koponen, S., 2018. Applicability of Earth Observation chlorophyll-a data in assessment of water status via MERIS — with implications for the use of OLCI sensors. Remote Sens. Environ. 212, 273–287. https://doi.org/10.1016/j.rse.2018.02.043. N. Kulha et al. Estuarine, Coastal and Shelf Science 298 (2024) 108634 8 Bates, D., Machler, M., Bolker, B., Walker, S., 2015. Fitting linear mixed-effects models using lme4. J. Stat. Software 67 (1), 1–48. https://doi.org/10.18637/jss.v067.i01. Bramante, J.F., Raju, D.K., Sin, T.M., 2013. Multispectral derivation of bathymetry in Singapore’s shallow, turbid waters. Int. J. Rem. Sens. 34 (6), 2070–2088. https:// doi.org/10.1080/01431161.2012.734934. Caballero, I., Stumpf, R., 2020a. Towards routine mapping of shallow bathymetry in environments with variable turbidity: contribution of sentinel-2A/B satellites mission. Rem. Sens. 12 (3), 451. https://doi.org/10.3390/rs12030451. Caballero, I., Stumpf, R., Meredith, A., 2019. Preliminary assessment of turbidity and chlorophyll impact on bathymetry derived from Sentinel-2A and Sentinel-3A satellites in south Florida. Rem. Sens. 11 (6), 645. https://doi.org/10.3390/ rs11060645. Caballero, I., Stumpf, R.P., 2019. Retrieval of nearshore bathymetry from Sentinel-2A and 2B satellites in South Florida coastal waters. Estuar. Coast Shelf Sci. 226, 106277 https://doi.org/10.1016/j.ecss.2019.106277. Caballero, I., Stumpf, R.P., 2020b. Atmospheric correction for satellite-derived bathymetry in the Caribbean waters: from a single image to multi-temporal approaches using Sentinel-2A/B. Opt Express 28 (8), 11742. https://doi.org/ 10.1364/OE.390316. Caballero, I., Stumpf, R.P., 2023. Confronting turbidity, the major challenge for satellite- derived coastal bathymetry. Sci. Total Environ. 870, 161898 https://doi.org/ 10.1016/j.scitotenv.2023.161898. Casal, G., Harris, P., Monteys, X., Hedley, J., Cahalane, C., McCarthy, T., 2020. Understanding satellite-derived bathymetry using Sentinel 2 imagery and spatial prediction models. GIScience Remote Sens. 57 (3), 271–286. https://doi.org/ 10.1080/15481603.2019.1685198. Casal, G., Monteys, X., Hedley, J., Harris, P., Cahalane, C., McCarthy, T., 2019. Assessment of empirical algorithms for bathymetry extraction using Sentinel-2 data. Int. J. Rem. Sens. 40 (8), 2855–2879. https://doi.org/10.1080/ 01431161.2018.1533660. Chybicki, A., 2017. Mapping South baltic near-shore bathymetry using Sentinel-2 observations. Pol. Marit. Res. 24 (3), 15–25. https://doi.org/10.1515/pomr-2017- 0086. Darecki, M., Stramski, D., 2004. An evaluation of MODIS and SeaWiFS bio-optical algorithms in the Baltic Sea. Remote Sens. Environ. 89 (3), 326–350. https://doi. org/10.1016/j.rse.2003.10.012. Darecki, M., Weeks, A., Sagan, S., Kowalczuk, P., Kaczmarek, S., 2003. Optical characteristics of two contrasting Case 2 waters and their influence on remote sensing algorithms. Continent. Shelf Res. 23 (3–4), 237–250. https://doi.org/ 10.1016/S0278-4343(02)00222-4. Doxani, G., Papadopoulou, M., Lafazani, P., Pikridas, C., Tsakiri-Strati, M., 2012. Shallow-water bathymetry over variable bottom types using multispectral WorldView-2 image. ISPRS - International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences XXXIX-B8, 159–164. https://doi.org/ 10.5194/isprsarchives-XXXIX-B8-159-2012. Evagorou, E., Mettas, C., Agapiou, A., Themistocleous, K., Hadjimitsis, D., 2019. Bathymetric maps from multi-temporal analysis of Sentinel-2 data: the case study of Limassol, Cyprus. Adv. Geosci. 45, 397–407. https://doi.org/10.5194/adgeo-45- 397-2019. Fleming-Lehtinen, V., Laamanen, M., 2012. Long-term changes in Secchi depth and the role of phytoplankton in explaining light attenuation in the Baltic Sea. Estuar. Coast Shelf Sci. 102–103, 1–10. https://doi.org/10.1016/j.ecss.2012.02.015. Gao, J., 2009. Bathymetric mapping by means of remote sensing: methods, accuracy and limitations. Prog. Phys. Geogr. Earth Environ. 33 (1), 103–116. https://doi.org/ 10.1177/0309133309105657. Geyman, E.C., Maloof, A.C., 2019. A simple method for extracting water depth from multispectral satellite imagery in regions of variable bottom type. Earth Space Sci. 6 (3), 527–537. https://doi.org/10.1029/2018EA000539. Guo, K., Li, Q., Wang, C., Mao, Q., Liu, Y., Zhu, J., Wu, A., 2022. Development of a single- wavelength airborne bathymetric LiDAR: system design and data processing. ISPRS J. Photogrammetry Remote Sens. 185, 62–84. https://doi.org/10.1016/j. isprsjprs.2022.01.011. Hastie, T., Tibshirani, R., Friedman, J., 2009. The Elements of Statistical Learning: Data Mining, Inference, and Prediction, second ed. Springer, New York, NY. https://doi. org/10.1007/978-0-387-21606-5. Hedley, J.D., Roelfsema, C., Brando, V., Giardino, C., Kutser, T., Phinn, S., Mumby, P.J., Barrilero, O., Laporte, J., Koetz, B., 2018. Coral reef applications of Sentinel-2: coverage, characteristics, bathymetry and benthic mapping with comparison to Landsat 8. Remote Sens. Environ. 216, 598–614. https://doi.org/10.1016/j. rse.2018.07.014. Isæus, M., 2004. Factors structuring Fucus communities at open and complex coastlines in the Baltic Sea. Doctoral dissertation, Botaniska institutionen. https://www.diva -portal.org/smash/record.jsf?pidˆdiva2%3A200596&dswidˆ9522. Knudby, A., Ahmad, S.K., Ilori, C., 2016. The potential for landsat-based bathymetry in Canada. Can. J. Rem. Sens. 42 (4), 367–378. https://doi.org/10.1080/ 07038992.2016.1177452. Kotta, J., Valdivia, N., Kutser, T., Toming, K., Ratsep, M., Orav-Kotta, H., 2018. Predicting the cover and richness of intertidal macroalgae in remote areas: a case study in the Antarctic Peninsula. Ecol. Evol. 8 (17), 9086–9094. https://doi.org/ 10.1002/ece3.4463. Kratzer, S., Hhkansson, B., Sahlin, C., 2003. Assessing Secchi and photic zone depth in the Baltic Sea from satellite data. Ambio 32 (8), 577–585. https://www.jstor.org/ stable/4315443. Krause-Jensen, D., Carstensen, J., Dahl, K., Back, S., Neuvonen, S., 2009. Testing relationships between macroalgal cover and Secchi depth in the Baltic Sea. Ecol. Indicat. 9 (6), 1284–1287. https://doi.org/10.1016/j.ecolind.2009.02.010. Kulha, N., Pasanen, L., Aakala, T., 2018. How to calibrate historical aerial photographs: a change analysis of naturally dynamic boreal forest landscapes. Forests 9 (10), 631. https://doi.org/10.3390/f9100631. Kutser, T., Paavel, B., Kaljurand, K., Ligi, M., Randla, M., 2018. Mapping shallow waters of the Baltic Sea with Sentinel-2 imagery. In: 2018 IEEE/OES Baltic International Symposium (BALTIC), pp. 1–6. https://doi.org/10.1109/BALTIC.2018.8634850. Kyryliuk, D., Kratzer, S., 2019. Summer distribution of total suspended matter across the Baltic Sea. Front. Mar. Sci. 5, 504. https://doi.org/10.3389/fmars.2018.00504. Lafon, V., Froidefond, J.M., Lahet, F., Castaing, P., 2002. SPOT shallow water bathymetry of a moderately turbid tidal inlet based on field measurements. Remote Sens. Environ. 81 (1), 136–148. https://doi.org/10.1016/S0034-4257(01)00340-6. Lappalainen, J., Virtanen, E.A., Kallio, K., Junttila, S., Viitasalo, M., 2019. Substrate limitation of a habitat-forming genus Fucus under different water clarity scenarios in the northern Baltic Sea. Estuar. Coast Shelf Sci. 218, 31–38. https://doi.org/ 10.1016/j.ecss.2018.11.010. Luhtala, H., Kulha, N., Tolvanen, H., Kalliola, R., 2016. The effect of underwater light availability dynamics on benthic macrophyte communities in a Baltic Sea archipelago coast. Hydrobiologia 776 (1), 277–291. https://doi.org/10.1007/ s10750-016-2759-x. Luhtala, H., Tolvanen, H., Kalliola, R., 2013. Annual spatio-temporal variation of the euphotic depth in the SW-Finnish archipelago, Baltic Sea. Oceanologia 55 (2), 359–373. https://doi.org/10.5697/oc.55-2.359. Lyons, M., Phinn, S., Roelfsema, C., 2011. Integrating Quickbird multi-spectral satellite and field data: mapping bathymetry, seagrass cover, seagrass species and change in Moreton Bay, Australia in 2004 and 2007. Rem. Sens. 3 (1), 42–64. https://doi.org/ 10.3390/rs3010042. Lyzenga, D.R., 1978. Passive remote sensing techniques for mapping water depth and bottom features. Appl. Opt. 17 (3), 379. https://doi.org/10.1364/AO.17.000379. Neumann, T., Koponen, S., Attila, J., Brockmann, C., Kallio, K., Kervinen, M., Mazeran, C., Müller, D., Philipson, P., Thulin, S., Vakeva, S., Ylostalo, P., 2021. Optical model for the Baltic Sea with an explicit CDOM state variable: a case study with Model ERGOM (version 1.2). Geosci. Model Dev. (GMD) 14 (8), 5049–5062. https://doi.org/10.5194/gmd-14-5049-2021. Pitarch, J., Volpe, G., Colella, S., Krasemann, H., Santoleri, R., 2016. Remote sensing of chlorophyll in the Baltic Sea at basin scale from 1997 to 2012 using merged multi- sensor data. Ocean Sci. 12 (2), 379–389. https://doi.org/10.5194/os-12-379-2016. Sagawa, T., Yamashita, Y., Okumura, T., Yamanokuchi, T., 2019. Satellite derived bathymetry using machine learning and multi-temporal satellite images. Rem. Sens. 11 (10), 1155. https://doi.org/10.3390/rs11101155. Stock, A., 2015. Satellite mapping of Baltic Sea Secchi depth with multiple regression models. Int. J. Appl. Earth Obs. Geoinf. 40, 55–64. https://doi.org/10.1016/j. jag.2015.04.002. Stumpf, R.P., Holderied, K., Sinclair, M., 2003. Determination of water depth with high- resolution satellite imagery over variable bottom types. Limnol. Oceanogr. 48 (1part2), 547–556. https://doi.org/10.4319/lo.2003.48.1_part_2.0547. Suominen, T., 2018. Applying MERIS time series and dynamic time warping for delineating areas with similar temporal behaviour in the northern Baltic Sea. Ecol. Indicat. 95, 794–804. https://doi.org/10.1016/j.ecolind.2018.08.023. Tronvig, K.A., 2005. Near-shore bathymetry. Hydro Int. 9, 24–25. Vahtera, E., Conley, D.J., Gustafsson, B.G., Kuosa, H., Pitkanen, H., Savchuk, O.P., Tamminen, T., Viitasalo, M., Voss, M., Wasmund, N., Wulff, F., 2007. Internal ecosystem feedbacks enhance Nitrogen-fixing cyanobacteria blooms and complicate management in the Baltic Sea. AMBIO A J. Hum. Environ. 36 (2), 186–194. https:// doi.org/10.1579/0044-7447(2007)36[186:IEFENC]2.0.CO;2. Vahtmae, E., Kutser, T., 2007. Mapping bottom type and water depth in shallow coastal waters with satellite remote sensing. J. Coast Res. 50, 6. Vahtmae, E., Kutser, T., 2016. Airborne mapping of shallow water bathymetry in the optically complex waters of the Baltic Sea. J. Appl. Remote Sens. 10 (2), 025012 https://doi.org/10.1117/1.JRS.10.025012. Venables, W.N., Ripley, B.D., 2002. Modern Applied Statistics with S, fourth ed. Springer https://www.stats.ox.ac.uk/pub/MASS4/. Virtanen, E.A., Viitasalo, M., Lappalainen, J., Moilanen, A., 2018. Evaluation, gap analysis, and potential expansion of the Finnish marine protected area network. Front. Mar. Sci. 5, 402. https://doi.org/10.3389/fmars.2018.00402. Werdell, P.J., Franz, B.A., Bailey, S.W., 2010. Evaluation of shortwave infrared atmospheric correction for ocean color remote sensing of Chesapeake Bay. Remote Sens. Environ. 114 (10), 2238–2247. https://doi.org/10.1016/j.rse.2010.04.027. Wood, S.N., 2017. Generalized Additive Models: an Introduction with R, second ed. Chapman and Hall/CRC. https://doi.org/10.1201/9781315370279. N. Kulha et al.