Spatiotemporal water quality analysis of Vrana Lake, Croatia

2.1 Study area and monitoring

Vrana Lake, located in Dalmatia near the eastern Adriatic coastline, is the largest natural lake in Croatia (Figure 1a), covering an area of approximately 30 km² [34], and is situated between 43°51′–43°57′N latitude and 15°30′–15°39′E longitude (WGS84 coordinate system). Its ecological characteristics and biodiversity have led to its protection as a nature park. The lake is characterized by shallow water that fluctuates seasonally, with higher levels in winter and spring (highest mean WL was 1.28 m above sea level (MASL) from 1990 to 2023) and lower levels in summer and autumn (lowest mean WL was 0.41 MASL from 1990 to 2023), based on data provided by the Croatian Meteorological and Hydrological Service (DHMZ) and visualized in Figure 2. The lake’s shallowness and geolocation make it particularly susceptible to wind influences that promote water column mixing [35]. In May 2023, heavy rainfall resulted in unusually high WL in subsequent months (Figure 2), highlighting the importance of considering meteorological and WL data in understanding water quality [36].

Figure 1

An overview of (a) Vrana Lake position and in situ stations, (b) preparation for measuring vertical profiles, and (c) and (d) in situ monitoring.

Figure 2

Comparison of AT, PP, and WL data for the research period of 2023/24 against historical averages.

Several factors contribute to the WL dynamics of Vrana Lake, including its connections to the Adriatic Sea, freshwater tributaries, brackish submerged groundwater discharge in the lake, underground karst fissures, and evaporation. During periods when WL drops below sea level, the lake’s salinity increases significantly due to the seawater influx through the artificial Prosika canal (Figure 1a). Additionally, the lake is affected by underground karst fissures along its western edge and the Jugovir brackish submerged groundwater discharge in the southern part, especially during strong south winds and high tides [37]. On the other hand, freshwater from surrounding karst fields and springs enters the lake through the Kotarka channel and the Lateral channel in the northern part (Figure 1a). Rainwater runoff from the northern and eastern hills is another significant source of water [32].

2.2 Data collection

In this study, we developed a comprehensive monitoring grid of 20 in situ stations for monthly in situ measurements of seven water quality parameters, including electrical conductivity (EC), turbidity, salinity, water temperature (WT), dissolved oxygen (DO), oxygen saturation (SO), and Chl-a, from July 2023 to June 2024 [38]. Furthermore, we analysed correlations between these parameters, WL, and six meteorological factors, such as AT, WS, PP, SH, relative humidity (RH), and atmospheric pressure (AP), in a 12-month period to understand their influence on lake ecosystem health. Additionally, we built upon the work of Rubinić and Katalinić [39] by analysing 34 years of WL and meteorological data (1990–2023), including AT and PP, to evaluate the more recent impact of climate change on the lake. Finally, we evaluated 15 GIS spatial interpolation methods to determine the best method for mapping water quality parameters, using cross-validation, root mean square error (RMSE), and mean error (ME). The flowchart is shown in Figure 3.

Figure 3

Flowchart illustrating the study’s methodology.

The lake’s ecological health is further impacted by anthropogenic influence from its tributaries [40] and influx of seawater via the Prosika canal, particularly during low WLs. To address these ecological concerns, we developed a new comprehensive monitoring grid designed to provide sufficient data for remote sensing and machine learning applications in subsequent research. Previously, monitoring efforts were limited to one station by the Water Institute Josip Juraj Strossmayer and three stations by the Public Institution Vransko Jezero Nature Park. We incorporated these existing stations and added new ones to ensure comprehensive spatial coverage. The 20 monitoring stations were strategically positioned (Figure 1a), considering lake’s characteristics, existing monitoring stations, the time required to complete the measurements by a vessel, and the hydrological model, including the lake’s bathymetry, tributaries, and connection to the sea [8].

We delineated the lake boundary using PlanetScope satellite imagery [41] and the Normalized Difference Water Index [42], focusing on the month of October 2023, when the WL was at its lowest during our year-long research period (Figure 2).

DHMZ provided historical monthly data on average AT and total PP from 1990 to 2023 for the Biograd na Moru station (43°56′44.5″N, 15°26′52.8″E, WGS84), along with average WL data from the Prosika station (43°51′0.3″N, 15°37′45.3″E, WGS84) (Figure 1a, Table S1). Additionally, they supplied daily meteorological data from the Biograd na Moru and Zadar stations, as well as WL data from the Prosika station, for our research period from July 2023 to June 2024 (Table 1 and Table S2). We analysed meteorological and WL parameters for each measurement day, except for PP, which we assessed based on the amount from 3 days prior. We adopted this methodology because we aimed to conduct in situ measurements on sunny, non-cloudy days, and we plan to integrate satellite data for specific measurement dates in our future research.

Researchers from Ruđer Bošković Institute collected in situ data on seven physicochemical and biological parameters in the lake (Table 1). Using the YSI EXO2 (YSI, USA) multiparameter probe, the team collected vertical profiles on a monthly basis from July 2023 to June 2024 (Figure 1c and d). These measurements were conducted between 8:00 and 13:00 at 20 designated monitoring stations (Table S3). Salinity is determined automatically from the sonde conductivity and temperature readings according to algorithms found in Standard Methods for the Examination of Water and Wastewater [43]. The use of the practical salinity scale results in values that are unitless, since the measurements are carried out in reference to the conductivity of standard water at 15°C [44]. Optical sensors in probes are increasingly important for real-time water quality assessment in environmental applications [45]. Due to poor weather conditions, we had to postpone our fieldwork scheduled for November 2023 to December 4th, and we conducted the December survey on the 19th, while all other measurements were conducted in their respective months.

2.3 Data analysis

3.1 Temporal trends of meteorological factors

The Jarque–Bera test confirmed that all historical monthly meteorological data exhibited a normal distribution annually, with p-values consistently exceeding the significance level of 0.05. Consequently, we used mean values to derive key annual metrics.

The WL values were historically low in 1990, 2008, and 2012, with measurements of 0.28, 0.32, and 0.31 m, respectively (Table S1, Figure 4). These low WL values corresponded with low PP levels in 1989, 2007, and 2011, measured at 555.50, 516.30, and 435.90 mm, respectively (Table S1, Figure 4). The analysis of WL and PP historical data over 34 years revealed an increasing WL trend of 3.12 mm per year and an increasing PP trend of 2.28 mm per year (Table S1, Figure 4). Average AT increased by 1.54°C during the same period, with a yearly increase trend of 0.05°C since 1990 (Table S1, Figure 4).

Figure 4

Annual average and trends of AT, PP, and WL during 1990–2023.

During our research period, elevated WL prevented significant interchange between sea and lake water, with the lowest monthly average measured at 0.69 MASL (Figure 2). WL values were higher than historical averages due to heavy rainfall in May 2023 and subsequently dropped below historic averages in January and February 2024, which we attributed to lower monthly PP at the end of 2023 and the beginning of 2024. This downward trend continued into early spring 2024, leading to further decreases in WL during May and June. Throughout the research, average monthly AT remained consistently higher than historic averages (Figure 2).

3.2 Assessing vertical and temporal variability in water column data

We collected vertical profiles at the northern lake station (Figure 1a – point marked as “Vertical profile” at 15°30′40″E, and 43°55′41″N, WGS84) on March 15th and July 19th, 2024. On March 15th, we measured a water depth of 2.95 m with a WL of 1.41 MASL, while on July 19th, the depth decreased to 2.12 m with a WL of 0.58 MASL (Figure S1). We also observed the temporal sensitivity of parameters at the same monitoring site at 30 min intervals (Figure S2). Despite the difference in water depth and time frame of the two measurements, we found that the readings from both dates were comparable in relative units, indicating temporal sensitivity with change per hour and vertical variability with change per meter.

We observed that salinity and EC remained stable across all depths and time throughout our campaigns, with negligible differences of ≤0.01 in their respective units between the lowest and highest values for each parameter. We observed a slight increase in DO over time, with rates of 0.04 mg/L per hour in March and 0.12 mg/L per hour in July, while depth had little effect on it (decrease of 0.02 mg/L per meter in March and decrease of 0.07 mg/L per meter in July). WT and SO decreased with depth, showing a decrease of 0.44°C per meter for WT and 1.11% per meter for SO in March and a decrease of 0.17°C per meter for WT and 1.15% per meter for SO in July. However, both parameters increased over time due to changes in AT, with an increase of 0.31°C per hour for WT and 1.20% per hour for SO in March and increase of 0.21°C per hour for WT and 1.73% per hour for SO in July. We also recorded fluctuations in Chl-a and turbidity, with Chl-a increasing by 0.45 µg/L per meter in March and 0.35 µg/L per meter in July. Turbidity experienced an increase of 0.52 FNU per meter in March and 0.49 FNU per meter in July, with gradual decline over time (decrease of 0.10 FNU per hour in March and 0.03 FNU per hour in July).

3.3 Seasonal and spatial variations in water quality parameters and their correlation with meteorological data

During the research period, all meteorological and WL data, except PP, exhibited normal distribution based on the Jarque–Bera test (Table S2). We transformed the PP variable using a cube root to achieve a normality. All water quality parameters had p-values below the 0.05 significance threshold, indicating non-normal distribution. While various transformations yielded marginal improvements, only the Chl-a parameter gained normal distribution following a square root transformation. Due to the predominance of non-normal distributions, we conducted our analysis of monthly water quality data using median values. We employed the median absolute deviation (MAD) to assess variation, as it is less sensitive to outliers compared to standard deviation [30]. Table S3 presents a comprehensive overview of the data distribution, presenting mean and median values alongside standard deviation and MAD. Our calculations were based on a maximum of 20 monitoring stations in the lake per date, and the quantity of measured data stations per date varied (Table S3).

We observed that Chl-a varied from 0.15 μg/L in April 2023 to 2.39 μg/L in November 2023, with a mean value of 0.94 ± 0.39 μg/L and median value of 0.90 ± 0.24 μg/L (Table S3). The range of Chl-a levels was 2.24 μg/L, with the maximum and minimum values differing almost 2.5 times the median. In June 2024, we measured the lowest EC value of 2.66 dS/m, which peaked in November 2023 at 4.80 dS/m, showing a gradual decrease by the end of our research period (Table S3, Figure 5). The mean and median EC measurements were 3.96 ± 0.65 dS/m and 4.19 ± 0.49 dS/m, respectively, with a range differing by 2.14 dS/m. This range represents 45 and 51% of the highest value and median value, respectively. The mean and median in situ values for salinity were 2.10 ± 0.36 and 2.22 ± 0.28, respectively, with values ranging from 1.37 in June 2024 to 2.58 in November 2023 (Table S3, Figure 5). The salinity values showed a 1.21 variation, accounting for 47 and 55% of the highest and median values, respectively.

Figure 5

Median values of water quality parameters from the dataset (Chl-a, EC, turbidity, salinity, WT, DO, SO).

We measured the turbidity levels ranging from 0.73 FNU in July 2023 to 15.93 FNU in August 2023, with a mean value of 5.76 ± 3.96 FNU and median value of 3.85 ± 1.63 FNU. Over this 12-month period, we observed a range of turbidity values of 15.20 FNU, which accounts for 95% of the maximum value and 395% of the median value. In July 2023, we measured the highest WT (30.22°C), which declined to 4.47°C in January 2024, with a mean of 16.98 ± 7.87°C and a median of 17.16 ± 7.33°C. The WT data exhibited a fluctuation of 25.75°C, representing 85 and 150% of the maximum and median values, respectively.

We observed that DO levels ranged from 7.60 mg/L in August 2023 to 12.60 mg/L in January 2024, with a mean value of 10.33 ± 1.15 mg/L and a median of 10.35 ± 0.85 mg/L. During this period, DO values fluctuated within a range of 5 mg/L, which represented 40 and 48% of the maximum and median values, respectively. We observed that SO decreased to 93.46% in December 2023 after reaching a high of 143.28% in July 2023. A mean of SO is 106.70 ± 8.85% and a median is 103.64 ± 4.95%, while the SO values ranged from 49.82%, which is equivalent to 35% of the highest value and 48% of the median value.

Spearman’s correlation analysis revealed that the daily vertical movement of phytoplankton (Figure S1), influenced by factors like SH, WS, water column mixing, and PP, leads to changes in Chl-a concentrations (Figure 6a). Notably, we observed higher Chl-a levels in November compared to July (Figure 5), but we found no significant correlation between Chl-a and DO levels (Figure 6b), likely due to unpredictable fluctuations from factors such as the absence of macrophytes and phytoplankton accumulation. Our study revealed that EC, salinity, and turbidity significantly impact the Chl-a values across the lake (Figure 6b). We found that factors such as WT, total dissolved solids, PP, flooding, evaporation, and water flow have an impact on EC and salinity levels. Increased water volume and level can reduce EC and salinity (Figure 6a).

Figure 6

Spearman’s coefficient of determination between: (a) WL, meteorological (AT, WS, RH, AP, PP, SH), and water quality parameters (Chl-a, EC, turbidity, salinity, WT, DO, SO) and (b) in-between water quality parameters (Chl-a, EC, turbidity, salinity, WT, DO, SO).

We found that turbidity can lead to increased WT and decreased DO levels (Figure 6b). We recognized that factors such as PP and WS can further impact turbidity by increasing stream volume and resuspending settled sediments (Figure 6a). Turbidity will often spike annually due to spring rains (Figure 5) and snowmelt. Sunlight and AT primarily influence WT in Vrana Lake (Figure 6a). We measured higher WT during summer, autumn, and spring, and lower temperatures during colder months (Figure 5). Warmer waters have higher EC, but colder waters can hold more DO and have lower levels of SO (Figure 6b).

Aeration sources such as wind, AT (Figure 6a), photosynthetic activity, and oxygen consumption by aquatic organisms influence the levels of DO in waterbodies. DO levels can vary based on WT (Figure 6b), water pressure, and salinity, with fluctuations occurring due to microbial breakdown and limited air contact. DO levels are higher in winter and lower in summer (Figure 5). SO levels are around or slightly exceed 100% throughout the year (Table S3 and Figure 5).

3.3.1 Spatial distribution

According to equation (3), we found that the coarsest grid resolution is approximately 300 m when considering 20 unevenly distributed monitoring stations in an area of approximately 30 km². The comparison of 15 GIS spatial interpolation methods revealed that the Global Polynomial Interpolation – Second order method achieved the lowest RMSE for Chl-a, while SK – Optimized achieved the lowest RMSE for EC, salinity, DO, and SO (Table S4, Figure S3). EBK – Advanced had the lowest RMSE for turbidity, and SK – Default had the lowest RMSE for WT (Table S4). Generally, Kriging methods outperformed deterministic interpolation methods regarding RMSE values, with EBK – Advanced and all SK methods ranking in the top five for average rank across all parameters. Overall, SK – Optimized ranked as the best for lowest RMSE using equation (1) across all parameters (Table S5). We calculated ME values using equation (2), and both Inverse Distance Weighted methods and EBK – Default consistently scored 0, indicating the highest accuracy (Table S4). However, when we considered the ranking of all GIS spatial interpolation methods, we found that SK – Trend emerged as the most suitable method for modelling all parameters (Table S5). We used this method to create distribution maps demonstrating the variability of water quality parameters (Figure 7).

Figure 7

SK – Trend interpolation showing parameter distributions: (a) Chl-a, (b) EC, (c) turbidity, (d) salinity, (e) WT, (f) DO, and (g) SO.

We observed that Chl-a concentrations peak in the northwestern region of the lake, particularly near the Kotarka channel, which serves as the main freshwater tributary (Figure 7a). Conversely, we found that the eastern part of the lake, being the shallowest (Figure 1a), consistently maintained lower Chl-a levels throughout the year (Figure 7a). Turbidity was the highest in the northwestern part of the lake, near the Kotarka channel, while the southeastern part exhibited the lowest turbidity levels (Figure 7c).

There was a strong correlation between EC and salinity in the lake, with both parameters exhibiting similar distribution patterns. The southern region, connected to the sea through the Prosika canal, experienced seawater intrusion, especially during warmer months with lower WL. This resulted in higher EC and salinity in the south and lower EC levels in the north (Figure 7b and d).

We monitored the lake in a counterclockwise direction, starting at station 4 and ending at station 2 (Figure 1a), which took us approximately 3 h. During this process, WT and SO levels gradually increased, influenced by AT (Figure S2). This resulted in higher WT at later monitoring stations and in shallower areas on the northern and eastern sides of the lake (Figure 7e). DO levels were highest in the southern, deepest part of the lake near the Prosika canal and lowest in the northern area near the freshwater tributary Lateral channel (Figure 7f). Similarly, SO levels were highest in the deep western part of the lake and lowest in the shallow northeastern region (Figure 7g).

4.1 Impact of climate change on Vrana Lake

Historical data indicate rising global temperatures and extreme weather events, including heavy precipitation. Our study confirms the growing trend of overall AT from the early 1990s till 2020s (2.11°C in 34 years) anticipated by IPCC [3], impacting ecosystems and affecting lake water quality. PP levels exhibit significant variability from average values, reflecting the occurrence of extreme weather events throughout the years, such as drought years and flooding years. Rising WL trend may be influenced by the Adriatic Sea’s mean sea level rise of +2.6 mm/year from 1993 to 2019, as reported by Pandžić et al. [50]. Given Vrana Lake’s connection to the sea through the Prosika canal, it is especially vulnerable to sea-level changes – a trend also observed in other coastal lakes, such as those in the Netherlands [51]. Our trend analysis, based on the work of Rubinić and Katalinić [39], which covers the period from 1961 to 2010, indicates a similar increase in AT and lake WL.

4.2 Vertical stratification and water quality dynamics

Vertical profiles revealed a very well-mixed water column with minimal or no stratification. This aligns with findings by Holgerson et al. [52], who reported that shallow lakes over 4 ha in surface tend to experience frequent mixing driven by wind and convection, preventing sustained stratification. In Vrana Lake, we measured slightly elevated Chl-a and turbidity at greater depths with lower light availability, a pattern also observed by Girdner et al. [53]. Similar to the observations in Lake Taihu [54] and Siombak Lake [55], where vertical WT differences were typically within 1°C and showed no significant variation in DO across the water column, we observed that Vrana Lake also exhibits very limited vertical WT gradients and minimal variations in DO, confirming its classification as a shallow, well-mixed system. Likewise, the stability of salinity throughout the water column mirrors observations in shallow saline Lake Shunet in Russia [56], where salinity remains stable to depths of 5 m. Given the shallow nature of Vrana Lake, we find that vertical variations are not significant and can be disregarded, allowing us to use median values from vertical profiles as reference measures for each monitoring station.

4.3 Seasonal and spatial variations in water quality and meteorological influences

The spatial distribution of water quality parameters is important for understanding changes in water quality across a lake. Our findings indicate a slight increase in Chl-a concentrations outside the vegetation period, particularly during autumn and winter. This is consistent with a study by Kong et al. [57], in which it was found that phytoplankton growth in winter is particularly sensitive to light and temperature changes, making it vulnerable to the impacts of climate change.

According to the legislation in The Official Gazette [58], our 12-month study of Vrana Lake shows that it is mesotrophic (classified as “very good” based on eutrophication indicators for mean annual Chl-a values). This differs from the 2021 measurements from Hrvatske vode, which classified the lake as mesotrophic/eutrophic (classified as “good”) [40]. We acknowledge that the reliability of Chl-a measurements in our study is uncertain due to the use of the multiparameter probe without spectrophotometric analysis of water samples. A study by Zolfaghari et al. [45] shows that the connections between sonde and laboratory measurements of Chl-a depend on the site and the methods used in the laboratory.

Salinity concentrations in Vrana Lake rise due to evaporation and the entry of seawater, particularly during the summer (Figure 5), as demonstrated by recent findings in another coastal lake in Croatia [59]. While the Lateral channel remains fresh, heavy rainfall can introduce slightly saline water into the Kotarka channel due to the spring from Vrana polje in the north [35]. Seasonal changes, such as reduced seawater influx during colder periods and potential salinity changes in the Kotarka channel, further affect water quality.

Turbidity levels in Vrana Lake are influenced by factors like nutrient and water flow from tributaries, with warmer months leading to increased turbidity due to lower WLs and active vegetation periods. Contributing sources include soil erosion, seasonal variations, local geology, and algal blooms, which can negatively impact water quality and aquatic life [60]. Additionally, higher turbidity level in the northwestern part of the lake may be linked to the anthropogenic influence, as reported by Silva et al. [61], including nearby road network and industries.

WT is affected by AT and SH, showing higher temperatures in shallow parts of the lake, as found by Anamunda and Lamtane [62]. WT variations have effect on the DO and SO levels, biological activities, and other parameters, as reported by Wang et al. [63] and Khouni et al. [19]. Areas with higher Chl-a concentrations generally have elevated DO levels, which inversely correlate with WT and AT, as found by Saturday et al. [64]. In shallow lakes like Vrana Lake, higher DO levels are typically associated with increased SO concentrations, and SO levels typically remain near 100% or slightly above, sustained by photosynthesis, aeration, and shallow mixing [65], as found by Allesson et al. [66].